\section{Introduction}
Partially connected hybrid structures for millimeter wave (mmW) massive MIMO systems are drawing considerable attention from researchers due to the benefits of spatial multiplexing coupled with low complexity, low power consumption and high energy efficiency \cite{ref_01}, \cite{ref_02}. However, such systems can only accommodate a limited number of users due to restrictions on the number of RF chains \cite{ref_03}. In the environments where the number of users is large compared to the number of RF chains, user selection is essential. In the existing literature, many model driven approaches are proposed for user or antenna selection. The performance of optimization and model driven methods for user selection is evaluated in terms of either bit-error-rate (BER) or achievable data rate. Evolutionary algorithms for user selection are employed in \cite{ref_04} to find a suboptimal solution with a reduced computational complexity compared to user selection based on exhaustive search. A greedy algorithm for user selection with downlink beamforming is applied in \cite{ref_05} to further reduce the computational complexity, but at the cost of a reduction in achievable rate. User selection based on matching theory is proposed in \cite{ref_06}. A method based on an Adaptive Markov Chain Monte Carlo (AMCMC) algorithm for multi-cell multi-user massive MIMO down-link systems is proposed in \cite{ref_07}. The authors in \cite{ref_08} utilized a generalized power iteration precoding (GPIP) algorithm to generate a joint solution for user selection, power allocation, and downlink precoding. However, all these algorithms are sub optimal in nature. The optimal solution to the user selection problem can be obtained by employing an exhaustive search (ES) algorithm which enumerates all possible combinations of active users, but leads to an exponential growth of the complexity with the increase of the number of active users. Thus, the ES algorithm is not suitable for massive MIMO systems due to high computational complexity.

Recently, data driven deep learning (DL) methods have shown great potential in dealing with the problems of channel estimation \cite{ref_09}, \cite{ref_10} and RF precoder design \cite{ref_11}. Similarly, data driven methods were applied to the transmit antenna selection problem, to obtain the best antenna subset \cite{ref_11}, \cite{ref_12}. Based on the availability of appropriate datasets, classifiers such as k-nearest neighbors (KNN) and support vector machine (SVM) can be trained to map the available channel dataset into the selected antenna subset. The author in \cite{ref_12} has proved that the performance of KNN and SVM based classifiers exceeds that of conventional optimization driven methods. Compared to optimization driven methods, data driven models are computationally efficient as well, which makes them convenient to be implemented in real-life applications. Note that, even though data driven methods have been widely applied to address multiple problems \cite{ref_09}-\cite{ref_12}, no work has been done so far on their application to user selection. Therefore, inspired by the efficiency and great potential of data driven techniques, in this letter we apply DL based convolutional neural networks (CNN) to the user selection problem in a massive MIMO wireless communication system. Our motivation can be explained as follows:
\begin{itemize}
	\item For the user selection problem, it is well known that the computational complexity of the ES algorithm grows exponentially with the number of active users. Hence, this algorithm is not appropriate for implementation in a real time environment.
	\item Traditional user selection algorithms require many serial iterations which generally result in a high complexity. However, in a neural network (NN), the training phase is generally performed offline, and therefore the online deployment has very low complexity and is limited to a few matrix multiplications and additions.
	\item The availability of a large dataset and the use of a huge number of iterations in the training phase enable CNNs to understand the complicated features of wireless channels.
\end{itemize}
\subsection{Novelty and Contribution}
In this letter, we propose a DL-based user selection design approach and develop a CNN which is trained using the channel realizations as the dataset and learns how to optimize the user selection to maximize the sum rate. The contributions can be summarized as follows:
\begin{itemize}
	\item \textit{New Design Approach}: A CNN based technique is proposed to directly perform the user selection without the involvement of any iterative algorithm after the CNN has been trained. The dataset is obtained after performing exhaustive search on each channel realization to select a subset of the active users which maximizes the sum rate. It is pertinent to mention that both training and dataset generation are performed offline.
	\item \textit{Robustness to imperfect CSI}: The proposed algorithm ensures robustness to imperfect channel state information (CSI). Firstly, during offline training, the CNN learns how to optimize the user selection to approach the ideal sum rate using practical channel estimates. Secondly, during the online deployment, the CNN can itself adapt to imperfect CSI and achieve robust performance in the presence of channel estimation errors.
\end{itemize}
\section{System Model}
\begin{figure}[tbph]
\centering
\includegraphics[height=1.8in,width=3.4in]{sys.eps}
\caption{System Model for User Selection}
\label{fig:sysmod}
\end{figure}
We consider a massive MU-MIMO downlink system with a transmitter equipped with $N_{T}$ transmit antennas, $N_{s}$ data streams and $N_{RF}$ RF chains serving $N_{R}$ single antenna users. We assume a multi-user precoding case where the transmitter communicates with each user via only one stream. Hence, the number of data streams is equal to the number of users. Further, we assume that the maximum number of users that can be simultaneously served by the transmitter equals the number of RF chains. We assume that $N_{R}$ $>$ $N_{RF}$. Hence, we have to choose a subset of $N_{r}$ users so that $N_{RF}$ = $N_{r}$. Our proposed CNN method will select $N_{r}$ users out of the available $N_{R}$ users and afterwards we will perform hybrid precoding for the selected users. For the sake of simplicity, let $N_{s}$ = $N_{RF}$ = $N_{r}$ = $N$. On the downlink, the transmitter applies an $N \times N$ baseband precoder $\mathbf{F}_{BB} = [\mathbf{f}^{BB}_{1}, \mathbf{f}^{BB}_{2}, ..., \mathbf{f}^{BB}_{N}]$ followed by an $N_{T} \times N$ RF precoder $\mathbf{F}_{RF} = [\mathbf{f}^{RF}_{1}, \mathbf{f}^{RF}_{2}, ..., \mathbf{f}^{RF}_{N}]$. The transmitted signal is represented as
\begin{equation}
    \mathbf{x} = \mathbf{F}_{RF} \mathbf{F}_{BB} \mathbf{s}
\end{equation}

where $\mathbf{s} = [s_{1}, s_{2}, ..., s_{N}]^{T}$ is the $N \times 1$ transmitted signal vector. Furthermore, $\mathbf{F}_{RF}$ and $\mathbf{F}_{BB}$ are evaluated using the algorithm proposed in \cite{ref_03}. We have adopted a narrowband block fading channel \cite{ref_13}, where the \textit{n}th user observes the received signal as
\begin{equation}
    \mathbf{r}_{n} = \mathbf{h}_{n} \sum_{u = 1}^{N} \mathbf{F}_{RF} \mathbf{f}_{u}^{BB} \mathbf{s}_{u} + \mathbf{n}_{u}
\end{equation}

where $\mathbf{h}_{n}$ denotes the channel between the transmitter and the \textit{n}th user. Since the mmW propagation environment has a limited number of scatterers, in this paper we follow the geometric Saleh-Valenzuela model to represent the low rank mmW channel \cite{ref_13}:
\begin{equation}\label{eq:chan}
    \mathbf{h}_{n} = \sqrt{\left(\dfrac{N_{T} N_{R}}{\epsilon L_{n}}\right)} \sum_{l=0}^{L_{n}}\eta_{n,l}\textbf{a}_{n,R}(\phi_{n,R},\theta_{n,R})\textbf{a}^{H}_{T}(\phi_{n,T},\theta_{n,T})
\end{equation}

$L_{n}$ represents the number of effective channel paths linked to the limited number of scatterers for each user, $\epsilon$ is the path-loss, $\eta_{n,l}$ is the complex gain associated with the \textit{l}th path, $\textbf{a}_{n,R}$ and $\textbf{a}_{n,T}$ are the spatial signatures of the receiver and the transmitter, respectively, $\phi_{R}$ and $\theta_{R}$ represent the elevation and azimuth angle of arrival (AoA) and $\phi_{n,T}$ and $\theta_{n,T}$ represent the elevation and azimuth angle of departure (AoD) of the \textit{l}th path at the receiver and at the transmitter, respectively. For uniform planar arrays (UPAs), the spatial signatures are estimated as given in \cite{ref_01}.

The aim is to select $N_{r}$ users from a total of $N_{R}$ single antenna users so that the sum rate is maximized. The received signal to interference plus noise ratio (SINR) for the \textit{k}th user after precoding is given as
\begin{equation}
	\gamma_{k} = \frac{\left| \mathbf{b}_{k}\mathbf{F}_{RF} \mathbf{f}_{k}^{BB} \right|^2} {\sum_{\acute{k}\neq k} \mathbf{b}_{k}\mathbf{F}_{RF} \mathbf{f}_{\acute{k}}^{BB} + \sigma^2}
\end{equation}

where $\sigma$ is the noise power and $\mathbf{b}_{k}$ is the channel between the transmitter and the \textit{k}th selected user. Furthermore, let $b_{k}$ represent the \textit{k}th row of $\mathbf{B}$ $\in C^{N_{r} \times N_{T}}$ which consists of rows of channel matrix $\mathbf{H}$ that correspond to the channel of selected users obtained using our proposed CNN architecture. The sum rate of the system after user selection is expressed as
\begin{equation}
	R = \log \bigg(\mathbf{I}_{N_{r}} + \sum_{k=1}^{N_{r}} \gamma_{k}\bigg)
\end{equation}
\begin{figure*}[htbp]
\centering
\includegraphics[height=1.6in,width=6.0in]{cnn.eps}
\caption{Deep CNN Architecture}
\label{fig:DL}
\end{figure*}

\section{CNN Based User Selection}
In this section, we first introduce the dataset and the label generation which will serve as the input to the CNN model for training. Then, we will describe the CNN architecture in detail.
\subsection{Data Normalization}
A channel matrix realization is a data sample. Since we have assumed single antenna users, each row of $\mathbf{H}$ represents the channel for the respective user. Since the CNNs do not accept complex entries as input, complex entries ($h_{ij}$) of channel matrix $\mathbf{H}$ are divided into real and imaginary parts and concatenated into a single matrix.
\subsection{Label Generation}
In order to select $N_{r}$ users from the $N_{R}$ available users, we have $W$ = $N_{R} \choose N_{r}$ combinations. Every combination is represented as a class pattern, so we have a total of $W$ class labels. We have a one to one correspondence between each set of a selected users and a class label. For example, the first pattern of combinations of selected users is given label 1, the second pattern of combinations of selected users is given label 2, and so on.

In order to generate a dataset for the training of our CNN, an exhaustive search is performed on each channel realization. The search algorithm examines all possible $W$ combinations generated and selects the combination of $N_{r}$ users that maximizes the sum rate. Hence, each channel realization is labeled with a class and given as input for the training of the CNN.
\subsection{Deep CNN Architecture} \label{cnn}
We have adopted a LeNet architecture for our CNN as shown in Fig \ref{fig:DL}. There are 2 convolutional layers, 2 pooling layers, 1 fully connected layer and 1 soft-max output layer. The input of the CNN is the normalized channel matrix. The first convolutional layer filters the normalized input channel matrix with 16 kernels of size $3 \times 3$. Then the first pooling layer normalizes and pools the input into a $3 \times 3 \times 16$ output response. The max-pooling kernels have a size of $2 \times 2$ and a stride of 2. The second convolutional layer filters the response with 32 kernels of size $3 \times 3$. The second pooling layer converts the input response to a $2 \times 2 \times 32$ output response. The fifth layer is a dense fully-connected layer with 1024 fully-connected kernels of size $1 \times 1$, which is a dense layer accelerating the convergence. Then a dropout layer, which randomly resets the output of each hidden neuron to zero, is added behind the fully-connected layer in order to avoid overfitting. The output of our CNN is the soft-max layer, which produces a class label. There is a one-to-one correspondence between a class label and a user selection. For the convolutional and fully-connected layers the rectified linear unit (ReLU) is used as the nonlinear activation function. The cross entropy is employed as the loss function with gradient descent optimizer.
\section{Simulations and Results}
Python 3.6 with tensor flow is used for CNN training and prediction, whereas the comparison with model driven algorithms is performed in MATLAB. The dataset for the simulations is generated for 100,000 independent channel realizations. The number of transmit antennas $N_{T}$ is set to 144. An exhaustive search method is used to generate the class label for each channel realization. The number of epochs for the training of CNN is kept as 200. Batch size is kept as 100. Fig \ref{fig:res3} shows the result obtained by selecting 6 out of 10 active users and Fig \ref{fig:res4} shows the result obtained by selecting 3 out of 10 active users.
\begin{figure}[tbph]
\centering
\includegraphics[height=1.6in,width=2.4in]{3.eps}
\caption{Spectral Efficiency vs SNR for $N_{R}=10, N_{r}=6$}
\label{fig:res3}
\end{figure}

\begin{figure}[tbph]
\centering
\includegraphics[height=1.6in,width=2.4in]{5.eps}
\caption{Spectral Efficiency vs SNR for $N_{R}=10, N_{r}=3$}
\label{fig:res4}
\end{figure}
The results verify the performance of our CNN based algorithm over conventional iterative algorithms. It can be seen that the CNN based algorithm outperforms the sum rate achieved by the iterative binary particle swarm optimization (BPSO) and greedy algorithms. The greedy algorithm has the worst performance due the fact that it cannot select the “best” set of users from the available users. Although our proposed CNN approach only marginally outperforms the BPSO algorithm in terms of sum rate, it achieves a large gain in terms of reduction in computational complexity, as discussed in Section \ref{ca}. Even though the proposed CNN algorithm has a performance gap with respect to the exhaustive search algorithm, the latter is not suitable due to high computational complexity.
\subsubsection{Performance Analysis with Imperfect CSI}
Next we evaluate the case of uncertainties affecting the channel model. The estimated channel matrix with imperfect CSI $\mathbf{\widehat{H}}$ can be modeled as \cite{ref_02}
\begin{equation}
    \mathbf{\widehat{H}} = \xi\textbf{H} + \sqrt{1-\xi^2}\textbf{E}
\end{equation}
where $\textbf{H}$ is the actual channel matrix, $\xi \in$ [0,1] is the CSI accuracy parameter and $\textbf{E}$ is an error matrix with entries following an i.i.d. $\mathbb{CN}$ (0,1). Fig. \ref{fig:resd} shows the effect of imperfect CSI on selecting 6 out of 10 active users with $\xi$ = 0.9 and 0.7. It is evident that the proposed system is not overly sensitive to CSI accuracy. It can also be seen that the achievable rate of the proposed system is quite close to that of the perfect CSI scenario even when $\xi = 0.7$.
\begin{figure}[tbph]
\centering
\includegraphics[height=1.6in,width=2.4in]{1.eps}
\caption{Impact of imperfect CSI on the proposed system}
\label{fig:resd}
\end{figure}
\subsection{Complexity Analysis}\label{ca}
In this section, we perform a complexity analysis of the proposed CNN method and compare it with that of conventional model driven methods. We consider $N_{T}$ = 144 transmit antennas and  $N_{R}$ = 10 total available users. For our proposed CNN, the complexity in the offline training stage is normally not counted \cite{ref_14}, and only that in the online deployment stage is counted. With the parameters mentioned in Section \ref{cnn}, the total number of operations for CNN is around 6 million.

The exhaustive search algorithm needs to compute $N_{R} \choose N_{r}$ complex determinants to evaluate the sum rate and each determinant requires O($N_{T}^{2}N_{r}^{2}$) complex operations to compute the hybrid beamforming matrices ($\mathbf{F}_{RF}$ and $\mathbf{F}_{BB}$) \cite{ref_14}. So for $N_{T}$=144, $N_{R}$=10 and $N_{r}$=6, the total number of operations is around 156 million. The BPSO based algorithm requires $N_{pop} \times N_{iter}$ determinants, where $N_{pop}$ is the population size and $N_{iter}$ is the total number of iterations for the BPSO algorithm. Both variables are set to 10 during the simulations. So, the total number of operations is around 74 million for the BPSO algorithm. The greedy algorithm needs to compute O($N_{T}^{2}N_{r}^{2}$) complex operations for a total number of $N_{R}$ users. So, the total number of operations is around 7.5 million for the greedy algorithm.

It can be seen that the proposed CNN has competitive computational complexity when compared with traditional model-based algorithms. In addition, the main operation of any CNN based algorithms only involves large-scale matrix multiplications and additions, which can be effectively accelerated by using graphics processing units. Conversely, most traditional model-based algorithms involve serial iterations where the optimization of the next iteration depends on the result of the previous iteration, which is not suitable for parallel computing.
\section{Conclusions}
We have proposed a CNN based design approach for user selection involving large antenna arrays. It is shown that the proposed CNN based approach is more efficient in terms of sum rate than model driven methods and performs close to the optimal exhaustive search based user selection. In addition, the computational complexity of data driven methods is acceptable and calculations only involve matrix multiplications and additions in the online phase.


\section{Introduction} \label{sec:introduction}

Parity-time ($\mathcal{PT}$) symmetry \cite{Bender1998} was first introduced as a fundamentally non-Hermitian extension of quantum mechanics \cite{Bender2002} and is now often used to model open systems \cite{El-Ganainy2007, Guo2009, Ruter2010} with symmetrically balanced loss and gain. The primary feature of $\mathcal{PT}$-symmetric systems is that up to a critical strength of gain and loss, the non-Hermitian system Hamiltonian has all real eigenvalues, and it shares its eigenstates with the $\mathcal{PT}$ operator. Above this threshold, pairs of complex conjugate eigenvalues form, breaking the $\mathcal{PT}$ symmetry.

Static systems with $\mathcal{PT}$ symmetry have been studied for the past several decades \cite{Feng2017, El-Ganainy2018, Ozdemir2019}, and have proven to be experimentally achievable in a wide variety of settings including waveguides \cite{Guo2009, Ruter2010, Zeuner2015, Kremer2019}, optical resonators \cite{Hodaei2014, Hodaei2016}, electrical circuits \cite{Schindler2011}, mechanical systems \cite{Bender2013}, and atomic systems \cite{Zhang2016}. While these static systems provide ample grounds to explore many of the basic properties of $\mathcal{PT}$-symmetric systems, they have several drawbacks. Chiefly, for extended systems, gain is difficult to implement in a way which exactly balances the loss, which often has its origins in the natural coupling of the system to its environment. For this reason, it is often convenient to transform a $\mathcal{PT}$-symmetric system to a loss-only one, which can be more easily implemented. Experimentally, large arrays of optical waveguides have proven to be useful in the study of these ``lossy'' $\mathcal{PT}$-symmetric systems \cite{Zeuner2015, Kremer2019}.

Building on this theory, a lattice model was proposed by Rudner and Levitov \cite{Rudner2009} to study topology in a $\mathcal{PT}$-symmetric system. This model is a natural extension of the Hermitian Su-Schrieffer-Heeger (SSH) model \cite{Su1979}, a simple one-dimensional dimer lattice which exhibits a well-known topological transition and topologically protected mid-gap edge states \cite{Asboth2016}. The model of Rudner and Levitov also exhibits many of these qualities; however, the associated energies of the localized edge states are pure imaginary, so that they are dynamically unstable \cite{Hu2011}. Despite this, the topological transition in this model is well-defined \cite{Liang2013, Harter2018, Lieu2018} and has been experimentally observed \cite{Zeuner2015}.

In the present paper, we study a time-periodic (Floquet) system \cite{Shirley1965} which is fashioned after that of Rudner and Levitov \cite{Rudner2009}, but with a modulated gain and loss rate \cite{Joglekar2014, Lee2015, Li2019, Leon-Montiel2018}. A similar system has previously been shown to exhibit an entirely real spectrum in the high-driving-frequency regime \cite{Yuce2015}. According to the Floquet theory, the long-time dynamics of such a time-dependent system cannot be determined by examining a single static system; rather, we may introduce an effective static system, which is determined by the evolution over one period of the modulation \cite{Shirley1965, Kitagawa2010}. The introduction of the periodic time dependence allows states which are dynamically unstable in the static case to be stabilized by the periodic modulation while introducing a novel topology unique to Floquet systems \cite{Rechtsman2013, Rudner2013}.

We begin, in Sec.~\ref{sec:static}, with a quick overview of the topologically relevant, Hermitian SSH model in Sec.~\ref{sec:ssh}, after which we discuss, in Sec.~\ref{sec:pt-ssh}, the static model proposed by Rudner and Levitov \cite{Rudner2009} which is its extension. We also provide a quick discussion of the analysis of this model in momentum space in Sec.~\ref{sec:momentum-space}. In Sec.~\ref{sec:floquet-pt}, we introduce the Floquet driving of the $\mathcal{PT}$-SSH model. We begin with a discussion of the Floquet theory in Sec.~\ref{sec:floquet}, followed by an analytical treatment in Sec.~\ref{sec:analytical-floquet}, and an analysis of the edge states in Sec.~\ref{sec:edge-states}. Finally, in Sec.~\ref{sec:conclusion}, we conclude with a brief discussion of the merits of our model.

\section{The Static \texorpdfstring{$\PT$}{PT}-SSH Model} \label{sec:static}

We first briefly review the relevant results of the static $\mathcal{PT}$-SSH model. To understand this non-Hermitian model, we begin with the Hermitian SSH model \cite{Su1979} in Sec.~\ref{sec:ssh}, which is a simple dimer lattice with alternating couplings. Then, in Sec.~\ref{sec:pt-ssh}, we add to this a $\mathcal{PT}$-symmetric extension \cite{Rudner2009, Zeuner2015} which attempts to respect the lattice symmetry.

\subsection{The SSH Model} \label{sec:ssh}

The starting point for our study is the Su-Schrieffer-Hegger (SSH) model \cite{Su1979}, a one-dimensional dimer chain, the bulk of which is described by the Hermitian tight-binding Hamiltonian
\begin{equation}\label{eqn:bulk-ssh}
    H_{SSH} = \sum_m\left(v\ket{m,A}\bra{m,B} + w\ket{m,B}\bra{m+1,A}\right) + \mathrm{h.c.}
\end{equation}
with intra-dimer coupling $v$ and inter-dimer coupling $w$, both real-valued quantities. Each dimer cell is labeled $m$ with sub-lattice labels $A$ and $B$. We further impose open boundary conditions with $m$ ranging from $1$ to $M$ and the summation of Eqn.~(\ref{eqn:bulk-ssh}) cut off appropriately so that there are exactly $M$ complete dimers.

\begin{figure}
    \centering
    \includegraphics[width=0.9\columnwidth]{ssh-bands}
    \caption{We show, for reference, the energy bands of the static SSH model plotted against the changing lattice parameter $v/v_T$. In the smaller inset, we show the wavefunction amplitudes of the two edge states with the eigenvalues $E = 0$ and a bulk state with $E = -J$ (plotted against each spatial site on the horizontal axis) for the parameterization corresponding to $v/v_t = 0.25$ (which is indicated by the vertical dashed red line in the larger plot). The edge states are highly localized to the left and right of the graph (red and blue), while the amplitude profile for a bulk state is highly delocalized (dark-gray).}
    \label{fig:ssh-bands}
\end{figure}

We explore parameterizations of this model in terms of its intra-dimer coupling $v$ by defining the scale of relevant coupling frequencies $v_T \equiv v + w$. Specifically, along the parameter range $v/v_T \in [0, 1]$, the couplings range between two fully dimerized limits $(v, w) = (0, v_T)$ and $(v, w) = (v_T, 0)$.

In Fig.~\ref{fig:ssh-bands}, we see that the bulk energy bands are gapped by $\Delta E = 2|v - w|$, so that at exactly $v = w$, namely the uniform lattice, the gap closes. Due to the topological nature of the lattice, when $v < w$, or equivalently $v/v_T < 0.5$, the boundary condition of the lattice introduces a pair of protected edge states (see the inset of Fig.~\ref{fig:ssh-bands}), and the lattice exhibits a topologically nontrivial phase, as can be seen by the presence of the mid-gap zero-energy eigenvalue. However, when $v > w$, or equivalently, $v/v_T > 0.5$, these states are absorbed into the bulk.

At the dimerization limits $v/v_T \rightarrow 0$ and $v/v_T \rightarrow 1$, we can easily inspect the topological nature of the real-space Hamiltonian. In the case of $v = 0$, both edges of the system are completely decoupled from the bulk of the system, resulting in two zero-energy eigenvalues given by the states $\ket{1,A}$ and $\ket{M,B}$, which are exactly localized on either edge of the system. In general, moving away from this limit, the edge states remain localized, as is demonstrated in the inset of Fig.~\ref{fig:ssh-bands}, until the gap closes at exactly $v/v_T = 0.5$. In the other limit, when $v/v_T = 1$ ($w = 0$), the unit cells are completely decoupled from each other. Its eigenstates are given by the direct products of independent dimers, and hence, no zero-energy edge states exist. Away from this limit, the dimers couple with each other, but the edge states remain nonexistent until the gap closes at $v/v_T = 0.5$.

\subsection{The \texorpdfstring{$\PT$}{PT}-Extended SSH Model} \label{sec:pt-ssh}

\begin{figure}
    \centering
    \includegraphics[width=0.9\columnwidth]{pt-ssh-lattice}
    \caption{Depiction of the $\mathcal{PT}$-SSH lattice with intra-dimer coupling $v$ and inter-dimer coupling $w$. The gain site ($+i\gamma$) is indicated by the red color, and the loss site ($-i\gamma$) is indicated by the blue color. Each pair of sites A and B (red and blue) comprises a full unit cell (labeled by $m$). Inside the cell, the sites are coupled by $v$, and the cells themselves are coupled by the alternate parameter $w$. We define the coupling scale by $v_T = v + w$.}
    \label{fig:pt-ssh-lattice}
\end{figure}

A non-Hermitian extension of this model was proposed by Rudner and Levitov \cite{Rudner2009}. In this model, an additional $\mathcal{PT}$-symmetric component which is local to each dimer is added to the original SSH model,
\begin{equation}\label{eqn:pt-ssh}
    H_{\mathcal{PT}} = H_{SSH} + i\gamma\sum_m \left( \ket{m,A}\bra{m,A} - \ket{m,B}\bra{m,B} \right) \,,
\end{equation}
as depicted in Fig.~\ref{fig:pt-ssh-lattice}, where $\gamma$ is a real parameter that controls the gain and loss rate in the system.

\begin{figure*}
    \centering
    \includegraphics[width=0.9\textwidth]{pt-ssh-bands}
    \caption{Depiction of the complex energy spectrum of the $\mathcal{PT}$-extended SSH model with fixed gain/loss rate $\gamma/v_T = 0.25$. In (a) and (e), we show the real and imaginary parts, respectively, of the energy spectrum as a function of the configuration parameter $v/v_T$. The vertical dashed lines correspond to the six values of $v/v_T$ highlighted by (b)-(d) and (f)-(h). In panels (b)-(d), we show the complex energy spectrum of the Hamiltonian for configurations having $v/v_T < 0.5$; each complex energy corresponds to a complex coordinate which is colored according to its IPR (from blue at zero to red at one). Specifically $v/v_T = 0.25$ in (b), $v/v_T = 0.35$ in (c), and $v/v_T = 0.45$ in (d). In panels (f)-(h), we show the complex energy spectrum for configurations having $v/v_T > 0.5$, with $v/v_T = 0.55$ in (f), $v/v_T = 0.65$ in (g), and $v/v_T = 0.75$ in (h).}
    \label{fig:pt-ssh-bands}
\end{figure*}

In this case, the Hamiltonian is no longer Hermitian, and the energy eigenvalues may take on imaginary components. In Fig.~\ref{fig:pt-ssh-bands}, we show the energies in the complex plane as a function of $v/v_T$ \cite{Lieu2018}. In Fig.~\ref{fig:pt-ssh-bands} (b)-(d), we observe the presence of non-bulk states (red), although they are pure imaginary. Between (c) and (d), the bulk states pick up imaginary components and begin to extend towards the isolated states on the imaginary axis. In Fig.~\ref{fig:pt-ssh-bands} (f)-(h), we observe the recession of the imaginary bulk states from (f) to (g); in (g), a gap begins to open between the purely real bulk bands, but the edge states are removed.

In order to see whether the states are localized or extended, we can also calculate the inverse participation ratio (IPR) for each state; the IPR, defined by
\begin{equation}
    \mathrm{IPR}(\psi) = \sum_m |\psi_m|^4 \,,
\end{equation}
is a measure of the extent to which a given state is localized. For the states indicated by blue dots in Fig.~\ref{fig:pt-ssh-bands} (f)-(h), the IPR is small $\leq 0.02$, whereas for the states with dark gray dots in (d), the IPR is $\approx 0.20$, and for the edge states indicated by red/yellow states in (b) and (c), the IPR is large $\approx 0.55$ (yellow) in (c) and $\approx 0.80$ (red) in (b). It is immediately clear that the localized modes are the ones which lie in the imaginary plane and they are complex conjugates.

Moving away from these limits, to the full domain of $v$ and $w$, this model was found (under the appropriate imaginary guage transformation) to admit a topological invariant which is quantized to $1$ for $v/v_T < 0.5$, and $0$ for $v/v_T > 0.5$ \cite{Rudner2009, Liang2013}. This invariant consists of a global integration over space and time, and is not associated to the upper or lower bands as is the case in the Hermitian SSH model \cite{Rudner2009, Lieu2018}. While there has been much discussion \cite{Rudner2009, Liang2013, Lieu2018, Kowabata2019} as to the nature of the topology of this and other non-Hermitian models, we focus our effort on the presence (or absence) of dynamically stable edge states.

In light of this, the $\mathcal{PT}$-extended SSH model does not have a topological phase with dynamically stable edge states \cite{Hu2011}; when $v/v_T < 0.5$, the edge states are inherently robust against perturbations in $v$, and separately, when $\gamma < |v - w|$, where the energy spectrum is entirely real (possessing unbroken $\mathcal{PT}$ symmetry), the system is dynamically stable; however, these two features do not coexist in this extended model. In the case where there are localized edges states, we can see that the system always has a pair of imaginary eigenvalues which correspond to the topological states; and thus the topological states are inherently dynamically unstable.

\subsection{Momentum Space Analysis} \label{sec:momentum-space}

Considering just the bulk of the system, let us neglect the boundaries for the moment. For a given intra-dimer coupling $v$ and gain/loss strength $\gamma$, the bulk Hamiltonian is block-diagonal in momentum space, with each block corresponding to a value $k \in [-\pi, \pi)$ having the form
\begin{equation} \label{eqn:hptk}
    H_\mathcal{PT}(k) = (v + w\cos k)\sigma_x + \left(w\sin k\right)\sigma_y + i\gamma\sigma_z \,,
\end{equation}
where $\sigma_i$ are the Pauli matrices. Therefore, each block in momentum space effectively represents a two-site system with a complex, $k$-dependent coupling. Furthermore, in a rotated frame which changes with $k$, we may write
\begin{equation} \label{eqn:hptkrot}
    H'_\mathcal{PT}(k) = r(k)\sigma_x + i\gamma\sigma_z \,,
\end{equation}
where
\begin{equation}\label{eqn:rk}
    r(k) = \sqrt{v^2 + w^2 + 2vw\cos k} \,
\end{equation}
is real and takes values between $|v-w|$ and $v + w$ over the range of $k$. The rotation of the frame is given by
\begin{equation} \label{eqn:krotation}
H_\mathcal{PT}(k) = e^{-i\sigma_z\phi(k)/2}H'_\mathcal{PT}(k) e^{i\sigma_z\phi(k)/2} \,,
\end{equation}
where we have defined the angle
\begin{equation}
    \phi(k) = \tan^{-1} \frac{w\sin k}{v + w\cos k} \,.
\end{equation}
The eigenvalues of the matrix of Eq.~(\ref{eqn:hptkrot}) are $E(k) = \sqrt{r^2(k) - \gamma^2}$, which indicates that when $\gamma \leq |v - w|$, $E(k)$ is real over the entire range of $k$.

Thus, we may determine the time development under the real-space matrix $H_\mathcal{PT}$ in terms of the time development of $H'_\mathcal{PT}(k)$ and we focus our attention to the $k$-space evolution. The time-evolution operator for the $k$th rotated block reduces to
\begin{align}
    \nonumber G_\mathcal{PT}(k,t) &= e^{-iH'_\mathcal{PT}(k)t} \\
    &= \cos[E(k)t]\,I - i\sin[E(k)t]\frac{H'_\mathcal{PT}(k)}{E(k)} \,,
\end{align}
where $I$ is the $2 \times 2$ identity matrix. It is important to note that the quantities $\cos[E(k)t]$ and $\sin[E(k)t] / E(k)$ are always real even when $\gamma > r(k)$ because $E(k)^2$ is real.

\section{The Floquet-Driven \texorpdfstring{$\PT$}{PT}-SSH model} \label{sec:floquet-pt}

It has previously been discovered by Yuce \cite{Yuce2015} that periodic temporal (Floquet) modulation can provide a means to completely restore stability to the topologically nontrivial phase in the high-driving frequency regime. In this paper, we propose a concrete, exactly-solvable, two-step modulation which will allow us to gain insights into this process of stabilization and explore regions of stability below the high-driving frequency. Furthermore, utilizing the simplicity of the model, we seek to identify the underlying symmetry. This setup has also proved to be experimentally beneficial as it only requires control of the gain or loss in pulsed fashion \cite{Li2019}.

Consider the following two-step time sequence for $0 \leq t < T$, where $T$ is the period of modulation,
\begin{equation} \label{eqn:pt-pt-driving}
    H(t) = \begin{cases}
        H_\mathcal{PT} & \quad 0 \leq t < T/2 \,, \\
        \tilde{H}_\mathcal{PT} & \quad T/2 \leq t < T \,,
    \end{cases}
\end{equation}
with $H_\mathcal{PT}$ being the Hamiltonian (see Eq.~(\ref{eqn:pt-ssh})) of the static $\mathcal{PT}$-extended SSH lattice described in Sec.~\ref{sec:pt-ssh}, while $\tilde{H}_\mathcal{PT}$ is the Hamiltonian of $H_\mathcal{PT}$ with the sign of $\gamma$ inverted to $-\gamma$. This corresponds to a periodic modulation of the gain/loss strength which instantaneously interchanges the location of gain and loss, $\gamma$, at each step. In fact, we can say that $\tilde{H}_\mathcal{PT}$ is the time-reversed version of $H_\mathcal{PT}$, as in $\tilde{H}_\mathcal{PT} = \mathcal{T}H_\mathcal{PT}\mathcal{T} = H_\mathcal{PT}^*$.

\subsection{Floquet Analysis} \label{sec:floquet}

According to the Floquet Theorem, the time evolution operator of a periodically driven system can be written
\begin{equation}
    G(t) = P(t)e^{-iH_Ft} \,,
\end{equation}
where the matrix $P(t)$ is periodic with the same frequency as the driving, while $H_F$ is a constant matrix known as the effective Floquet Hamiltonian. Thus the long-term dynamics of the evolution are determined completely by $H_F$; further, at stroboscopic times $nT$, for any integer $n$, $P(nT)$ is the identity because $G(0) = P(0) = I$. Consequently, the system is dynamically unstable when at least one eigenvalue of $H_F(k)$ has an imaginary component, which results in exponential growth or decay in the long-term evolution of the system.

\begin{figure*}
    \centering
    \includegraphics[width=0.9\textwidth]{pt-phases}
    \caption{The $\mathcal{PT}$ phase diagrams for the Floquet $\mathcal{PT}$-symmetric systems discussed in Sec.~\ref{sec:floquet-pt}. In each figure, the magnitude of the maximum imaginary eigenvalue is plotted over a plane in parameter space, where the color indicates increasing magnitude (thus, the black regions indicate $\mathcal{PT}$-symmetric phases). As a visual aid, a solid line has been added which outlines the border between the $\mathcal{PT}$-symmetric and $\mathcal{PT}$-broken phases. (a) The phase diagram for the two-site system of Eq.~(\ref{eqn:hpt2}), over a range of gain/loss parameters $\gamma$ and driving frequencies $\omega$. (b) The phase diagram for the 40-site Floquet $\mathcal{PT}$-extended SSH driven between $\gamma/v_T = 0.5$ and $\gamma = 0$ over a range of driving frequencies and couplings $v/v_T$. (c) The phase diagram over the $\omega$-$\gamma$ plane as in (a), but for the 40-site Floquet $\mathcal{PT}$-extended SSH model driven by swapping the gain ($+i\gamma$) and loss ($-i\gamma$) sites and with fixed coupling $v/v_T = 0.2$. (d) The phase diagram for the same system of (c), but in the $v$-$\omega$ plane with gain/loss magnitudes fixed at $\gamma/v_T = 0.2$. The dashed red lines in (c) and (d) indicate a portion of the $\mathcal{PT}$-phase with equivalent parameter values. In each panel, the dashed white lines indicate locations of the first three resonances $\omega/v_T = 2/n$ ($\omega/J = 2/n$ in the two-site case) with $n = 1, 3, 5$. Note that there are additional resonances between these in the asymmetrically driven case of panel (b).}
    \label{fig:pt-phases}
\end{figure*}

For a two-site Hamiltonian
\begin{equation} \label{eqn:hpt2}
H_2 = J\sigma_x + i\gamma\sigma_z \,,
\end{equation}
where $J$ is a real coupling parameter, similar Floquet $\mathcal{PT}$-symetric models have previously been studied \cite{Joglekar2014, Lee2015} and found to exhibit a rich $\mathcal{PT}$-phase diagram featuring re-entrant $\mathcal{PT}$-symmetry breaking depending on the driving frequency $\omega \equiv 2\pi/T$. As is shown in Fig.~\ref{fig:pt-phases}(a), at specific driving frequencies, i.e. $\omega/J = 2, 2/3, 2/5, \ldots$, the $\mathcal{PT}$ symmetry of the system is broken for any nonzero value of $\gamma$ (i.e. the $\mathcal{PT}$-threshold is zero at these locations). The periodic modulation at these driving frequencies is in resonance with the natural coupling frequency of the two-site system, namely $2J$.

Utilizing the similarity between Eq.~(\ref{eqn:hptkrot}) and Eq.~(\ref{eqn:hpt2}), we can explore the $\mathcal{PT}$-phase diagram of the Floquet $\mathcal{PT}$-SSH model. In this case, as seen in Fig.~\ref{fig:pt-phases}(c), the resonance points occur for each value of $k$ at, i.e. $\omega/r(k) = 2, 2/3, 2/5, \ldots$, with $r(k)$ as defined in Eq.~(\ref{eqn:rk}), and there is an apparent broadening of the zero-$\mathcal{PT}$-threshold regime. In this case, the $\mathcal{PT}$-phase diagram can be thought of as the composition of the two-site case of Fig.~\ref{fig:pt-phases}(a) for multiple scaling values $r(k)$ determined by the allowed momenta $k$ replacing $J$ in Eq.~(\ref{eqn:hpt2}). The vertical dashed lines are placed at the resonance frequencies of the two-site case, demonstrating that they correspond to the two-site scaled system with the choice of $r(k=0) = v_T$.

We can also analyze this model over a range of coupling parameters $v$ with a fixed value of gain/loss parameter $\gamma$, as is done in both Fig.~\ref{fig:pt-phases}(b), (d), where we show the $\mathcal{PT}$ phase in a sweep along the lattice parameter $v$ with driving frequencies $\omega$ increasing from bottom to top. For comparison to the previous two-site experimental studies \cite{Li2019, Leon-Montiel2018}, in Fig.~\ref{fig:pt-phases}(b), we modulate the system between a $\mathcal{PT}$-symmetric one ($\gamma > 0$) and a Hermitian one ($\gamma = 0$), instead of the modulation of Eq.~(\ref{eqn:pt-pt-driving}). In this case, when $v/v_T < 0.5$, we see that the $\mathcal{PT}$ symmetry is always broken, as is expected by the presence of edge states with pure imaginary eigenvalues.

Similarly in Fig.~\ref{fig:pt-phases}(d), we modulate the Hamiltonian between opposing $\mathcal{PT}$-symmetric systems as in Eq.~(\ref{eqn:pt-pt-driving}). Note that since Fig.~\ref{fig:pt-phases}(c) and (d) have the same driving type, there is correspondence between the two graphs, which we have demonstrated by the vertical dashed line in (d) and the horizontal dashed line in (c). We observe that in this case, when $v/v_T < 0.5$, regions of unbroken $\mathcal{PT}$ symmetry may arise (corresponding to the black areas of the plot). Above the high driving-frequency threshold, we see that the spectrum is completely real for all coupling choices, and thus resembles the Hermitian SSH model, in accordance with the previous study \cite{Yuce2015}. We reveal that, even below this threshold, there are other regions of $\mathcal{PT}$-symmetry with $v/v_T < 0.5$, which arise just above the corresponding driving-frequency resonance points and which still retain localized edge states. These resonance points are highlighted by the dashed white lines.

\subsection{Analytical Result} \label{sec:analytical-floquet}

To better understand the topology and stability, we wish to analytically explore the effective Floquet Hamiltonian. Because the time dependence is a simple two-step function, the time evolution of the system over one period of driving is the product of the propagators for the static Hamiltonians associated with each step:
\begin{equation}\label{eqn:monodromy}
    G(T) = e^{-i\tilde{H}_\mathcal{PT}\tau}e^{-iH_\mathcal{PT}\tau} \,,
\end{equation}
where we introduced the shorthand $\tau \equiv T/2$.

Thus, we move to the momentum space, where each 2$\times$2 block, as in Eq.~(\ref{eqn:hptk}), of the Hamiltonian associated with momentum $k$ is independent. Therefore we may define a rotation which rotates each corresponding subspace around the $z$-axis such that we can identify the $k$-dependent rotated Hamiltonian,
\begin{equation}
    H_\pm(k) = r(k)\sigma_x \pm i\gamma\sigma_z = \vec{\lambda}_\pm(k) \cdot \vec{\sigma} \,,
\end{equation}
where $H_\mathcal{PT}(k) = H_+(k)$ and $\tilde{H}_\mathcal{PT}(k) = H_-(k)$ correspond to the two steps of the Floquet sequence, while $\vec{\lambda}_\pm$ are the vectors $[r(k), 0, \pm i\gamma]^T$ respectively (with $^T$ the vectorial transposition operation) and $\vec{\sigma} = [\sigma_x, \sigma_y, \sigma_z]^T$ is the Pauli vector such that $H_\pm = \vec{\lambda}_\pm(k)\cdot\vec{\sigma}$. We may write $\vec{\lambda}(k)_\pm = E(k)\hat{\lambda}_\pm(k)$, where $E(k) = \sqrt{r(k)^2 - \gamma^2}$ is a (possibly complex) energy eigenvalue of both $H_\pm(k)$, whereas $\hat{\lambda}_\pm(k)$ is a (possibly complex-valued) unit vector.

Dropping the argument $k$ from $E(k)$ and $\vec{\lambda}_\pm(k)$ for convenience, we may then obtain the time evolution operator for evolution of a rotated, $k$-space block, up to one driving period $T = 2\tau \equiv 2\pi/\omega$ (with driving frequency $\omega$) Eq.~(\ref{eqn:monodromy}):
\begin{align}
    \nonumber G(T) &= (\cos{E\tau} - i \sin{E\tau} \hat{\lambda}_+\cdot\vec{\sigma})(\cos{E\tau} - i \sin{E\tau} \hat{\lambda}_-\cdot\vec{\sigma}) \\
    &= \left(\cos^2{E\tau} - \sin^2{E\tau} \frac{r^2 + \gamma^2}{E^2}\right)I \nonumber \\
    &+ 2\frac{r}{E}\sin{E\tau} \left(\cos{E\tau} \hat{x} - i \frac{\gamma}{E}\sin{E\tau} \hat{y}\right)\cdot\vec{\sigma} \,,
\end{align}
where $\hat{x} = [1, 0, 0]^T$, $\hat{y} = [0, 1, 0]^T$ are unit vectors and we used the formula $(\vec{a}\cdot\vec{\sigma})(\vec{b}\cdot\vec{\sigma}) = (\vec{a}\cdot\vec{b})I + i(\vec{a}\times\vec{b})\cdot\vec{\sigma}$. From this form, we can easily identify the Floquet Hamiltonian $H_F$ and its eigenvalues (up to arbitrary terms $n\omega$ for integer $n$) $\pm\mathcal{E}$ in the form $\mathcal{E}\hat{H}_F\cdot\vec{\sigma}$, where $\mathcal{E}$ is found by solving:
\begin{equation}\label{eqn:floqeig}
    \cos(2\mathcal{E}\tau) = 1 - 2\sin^2(\mathcal{E}\tau) = 1 - 2\frac{r^2}{E^2}\sin^2{E\tau} \,,
\end{equation}
and $\hat{H}_F$ is found from
\begin{equation}\label{eqn:floqvec}
    \sin(2\mathcal{E}\tau)\hat{H}_F = 2\frac{r}{E}\sin{E\tau} \left(\cos{E\tau} \hat{x} - i \frac{\gamma}{E}\sin{E\tau} \hat{y}\right) \,,
\end{equation}
by utilizing the value of $\mathcal{E}$ found in Eq.~(\ref{eqn:floqeig}). Then the effective Floquet Hamiltonian is fully described by $H_F = \mathcal{E}\hat{H}_F$.

By the form of Eq.~(\ref{eqn:floqvec}), we see that $H_F$ can be written
\begin{equation} \label{eqn:hf}
    H_F = c\left(\cos{E\tau}\sigma_x - i\frac{\gamma}{E}\sin{E\tau}\sigma_y\right)\,,
\end{equation}
with $c \equiv (2\mathcal{E}/\sin{2\mathcal{E}\tau})(\sin{E\tau}/E)$. Restricting the real part of $\mathcal{E}$ to $[-\omega/2, \omega/2)$, we see that there are only two cases to handle: either $\mathcal{E}$ is real, or, when $|(r/E)\sin{E\tau}| > 1$, $\mathcal{E} = \omega/2 + i\eta$ where $\eta$ is real.

We see that when $\mathcal{E}$ is real, $c$ is real, and $H_F$ has the symmetries
\begin{align}
    \sigma_z H_F \sigma_z &= -H_F \,, \label{eqn:sublattice}\\
    \sigma_x H_F \sigma_x &= H_F^\dagger\,, \label{eqn:pseudohermitian}\\
    \sigma_y H_F \sigma_y &= -H_F^\dagger\, \label{eqn:chiral},
\end{align}
which correspond to sublattice symmetry for Eq.~(\ref{eqn:sublattice}), pseudo-Hermitian symmetry for Eq.~(\ref{eqn:pseudohermitian}), and chiral symmetry respectively for Eq.~(\ref{eqn:chiral}); and consequently, the Hamiltonian belongs to the non-Hermitian symmetry class AIII according to Ref.~\cite{Kowabata2019}.

When $\mathcal{E}$ is complex, so is $c$, and the latter two symmetries of Eqs.~(\ref{eqn:pseudohermitian}) and (\ref{eqn:chiral}) are removed. However, as previously mentioned, the solutions to Eq.~(\ref{eqn:floqeig}) also admit the solutions $\mathcal{E} + n\omega$ for any integer $n$. Rather than considering the eigenvalue pair $\pm\mathcal{E}$ which correspond to $\omega/2 + i\eta$ and $-\omega/2 - i\eta$, we may consider a solution which, for example, gives the complex conjugate pair $\omega/2 + i\eta$ and $\omega/2 - i\eta$. An alternate form of $H_F$ which takes into account this situation is
\begin{equation}
    H_{F,s} \equiv V\left[\mathcal{E}_s\sigma_z + \frac{\omega}{2} I\right]V^{-1} \,,
\end{equation}
where $\mathcal{E}_s = \mathcal{E} - \omega/2$ is the shifted eigenvalue and $V$ is the matrix of eigenvalues such that $H_F = \mathcal{E}V\sigma_zV^{-1}$. Importantly, when $\mathcal{E}$ is complex, $\mathcal{E}_s$ is pure imaginary. While this Hamiltonian is not traceless, as in Eq.~(\ref{eqn:hf}), it admits the pseudo-Hermitian symmetry of Eq.~(\ref{eqn:pseudohermitian}) even when the eigenvalues are complex. Thus, for all parameterizations of the model, the effective Hamiltonian can be seen to have this symmetry. Starting with a Floquet-modulated $\mathcal{PT}$-symmetric Hamiltonian, we have obtained (in the Floquet picture) a non-Hermitian Hamiltonian with pseudo-Hermitian symmetry \cite{Mostafazadeh2002}.

\subsection{Analysis of the Edge States} \label{sec:edge-states}

We also may analyze the fate of the edge states in the Floquet-driven $\mathcal{PT}$-SSH lattice.

\begin{figure*}
    \centering
    \includegraphics[width=0.9\textwidth]{floq-pt-ssh-bands}
    \caption{The complex energy spectrum of the Floquet $\mathcal{PT}$-symmetric system of Sec.~\ref{sec:floquet-pt} as a function of $v/v_T$, with driving frequency $\omega/v_T = 0.7$ and gain/loss rate $\gamma/v_T = 0.2$. The complex energies are plotted in the complex plane, with each marker colored by its IPR (blue at zero to red at one). In panels (b)-(d), the spectrum of the Hamiltonian is shown on the complex plane for specific configurations below the topological transition $v/v_T < 0.5$, (b) $v/v_T = 0.1$, (c) $v/v_T = 0.2$ and (d) $v/v_T = 0.4$. (f)-(h) the complex energy spectrum above the topological transition $v/v_T > 0.5$, (f) $v/v_T = 0.6$, (g) $v/v_T = 0.8$, and (h) $v/v_T = 0.9$. (a) and (e) the real and imaginary parts, respectively, of the energy bands as a function the configuration parameter $v/v_T$. The vertical dashed lines correspond to the six values of $v/v_T$ highlighted by (b)-(d) and (f)-(h).}
    \label{fig:floq-pt-ssh-bands}
\end{figure*}

For modulation between a $\mathcal{PT}$-symmetric Hamiltonian and a Hermitian one, as in Fig.~\ref{fig:pt-phases}(b), we find that the energy eigenvalues corresponding to the edge states lie in the complex plane, and the system never has localized edge states in the $\mathcal{PT}$-symmetric phase. This is a result of the inability for the imaginary parts of the associated energies to balance one another during a single oscillation. Just as in the static $\mathcal{PT}$-SSH case, the edge states which correspond to the resulting imaginary eigenvalues occur when $v/v_T < 0.5$ and are highly localized. In fact, the result of the static $\mathcal{PT}$-SSH model with gain/loss parameter $\gamma$ is recovered identically in the high-frequency regime when the Floquet driving is between $\gamma_1=2\gamma$ and $\gamma_2=0$, i.e. the average $(\gamma_1 + \gamma_2)/2$ is equal to $\gamma$.

In Fig.~\ref{fig:pt-phases}(d), we can see that in the case of $\mathcal{PT}$-symmetric to $\mathcal{PT}$-symmetric driving, the same $\mathcal{PT}$ phase diagram produces regions of dramatic stability within the topologically nontrivial phase. In this case, the alternation of the source and drain allows for cancellation in a single driving oscillation, i.e. the average of the gain/loss parameter over one period is zero, thus protecting the system from breaking the $\mathcal{PT}$ symmetry in the presence of the edge modes. In this case, in Fig.~\ref{fig:floq-pt-ssh-bands}, we see that the result of the static SSH model can be recovered, where the edge state corresponds to a real, zero-energy eigenvalue, and the system is simultaneously $\mathcal{PT}$-symmetric and has $v/v_T < 0.5$.

In Fig.~\ref{fig:floq-pt-ssh-bands}(b)-(d), we observe the presence of real, mid-gap states (red/yellow). Between (b) and (c), the gap closes, but remains real, and the edge states remain localized inside the bulk spectrum. Between (d) and (f), the bulk states enter the imaginary portion of the complex plane, but again the edge states remain real. In Fig.~\ref{fig:floq-pt-ssh-bands}(f)-(h) we observe the reversal of the changes from (b)-(d), except that the edge state is removed; from (f) to (g), the bulk spectrum becomes real again, and from (g) to (h), a gap opens between the real bulk bands.

Specifically, for the complex energies depicted in Fig.~\ref{fig:floq-pt-ssh-bands}(b)-(d), the localization (IPR) of the corresponding states is indicated by the color of each marker, with the zero-energy edge states depicted with IPR $\approx 0.98$ (red) in (b), $\approx 0.88$ (red) in (c), and $\approx 0.38$ (yellow) in (d). In Fig.~\ref{fig:floq-pt-ssh-bands}(f)-(h), the highest corresponding IPR is below $0.01$.

It is important to remember that the energy spectrum shown here is restricted to values which are between $\pm\omega/2$, and that this spectrum pattern can be repeated infinitely along the real axis in both directions. We also observe that the localized, zero-energy mode is present even after the bulk gap has closed in the complex plane and exists within the bulk.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{floq-pt-ssh-states}
    \caption{State amplitudes (by magnitude) of the Floquet $\mathcal{PT}$ SSH model. In (a), we show an example with $v/v_T = 0.2$, for which it is clear that there exist localized edge states (shown in blue and red) along with the typical bulk states (dark gray). In (b), we show an example with $v/v_T = 0.8$; we have plotted the state with the highest IPR, to demonstrate that it is not localized to either edge. In both cases, the driving frequency is $\omega/v_T = 0.7$ and the gain/loss rate is $\gamma/v_T = 0.2$}
    \label{fig:floq-pt-ssh-states}
\end{figure}

In Fig.~\ref{fig:floq-pt-ssh-states}, we provide a specific example to demonstrate the localization of the edge states in this model. In (a), we show the site amplitudes (by magnitude) of the edge states as well as a bulk state for a realization with $v/v_T < 0.5$. Here, the localization of the left (blue) and right (red) edge states is clear, while the bulk (dark gray) remains highly delocalized. Similarly in (b), for a realization with $v/v_T > 0.5$, we show the site amplitudes for the state with the highest IPR; we specifically note the clear absence of localized edge states in this case.

\section{Conclusion} \label{sec:conclusion}

In conclusion, we have confirmed that Floquet modulation is able to provide a stabilizing effect to the static $\mathcal{PT}$-SSH model by stabilizing the the edge modes. The choice of modulation has nontrivial effects on the presence of the $\mathcal{PT}$-symmetric phase. When the modulation is kept $\mathcal{PT}$ symmetric for half a period and then Hermitian for the other half, in the Floquet picture, we see that, at high enough frequencies, the effective Hamiltonian has an eigenspectrum reminiscent of the static $\mathcal{PT}$-SSH lattice with complex edge modes. However, when we drive between a pair of time-reversed $\mathcal{PT}$-symmetric lattices, the high-frequency spectrum of the effective Hamiltonian resembles that of the static Hermitian SSH lattice with real and stable edge modes.

Above this high-driving frequency, the Floquet effective Hamiltonian has all real eigenvalues for all choices of lattice parameter up to a critical gain/loss rate (the $\mathcal{PT}$-symmetry breaking threshold). Importantly, below this high frequency regime, we find many other regions of stability. Specifically, below the first resonance frequency ($\omega/v_T = 2$) but above the second ($\omega/v_T = 2/3$), there can exist a broad portion of the $\mathcal{PT}$-phase diagram which is still in the $\mathcal{PT}$-unbroken phase and which has real, zero-energy, localized edge states.

Finally, we stress that this model is amenable to experimental implementation. Because of the simple, temporally-periodic switching used as the driving pattern, the burden of matching the gain/loss rates to a smooth, continuous function is lifted. Possibilities include purely loss-only systems where the location of the loss is alternated; after a simple imaginary gauge transformation, such systems map directly onto the one we have discussed here.

\begin{acknowledgments}
    The authors would like to acknowledge the financial support of Japan Society for the Promotion of Science (JSPS) Grants-in-Aid for Scientific Research (KAKENHI) JP19F19321 (A. H.) and JP19H00658 (N. H.). A. H. also acknowledges support as an International Research Fellow of JSPS (Postdoctoral Fellowships for Research in Japan (Standard)). Both authors would also like to thank Kohei Kowabata and Yogesh N. Joglekar for their helpful and insightful discussions and comments.
\end{acknowledgments}


\section{Introduction}

Face recognition technology has become ubiquitous these days, thanks to the advance in deep learning based methods \cite{jain2011handbook}. However, the reliability of face recognition systems (FRS) in the presence of attacks is poor in practical situations, mainly due to the vulnerability against presentation attacks (a.k.a spoofing attacks) \cite{handbook2}, \cite{ISO}. Presenting artifacts like a photograph or video in front of the camera could be enough to fool unprotected FR systems. The artifact used for this attack is known as a presentation attack instrument (PAI).

Presentation attack detection (PAD) tries to protect the FR systems against such attacks. Majority of the PAD methods proposed in the literature focuses on the detection of 2D attacks. Most of these methods use either feature-based methods or CNN based approaches using visible spectrum images. However, they do not emulate a realistic scenario of encountering realistic 3D attacks. Recently, there has been a lot of works focusing on the detection of a wide variety of attacks, including 2D, 3D, and partial attacks. Multi-channel methods have been proposed as a possible alternative to deal with real-world scenarios with a wide variety of attacks \cite{george_mccnn_tifs2019, nikisins2019domain}. Usage of multiple channels makes it harder to fool the PAD systems.

\begin{figure}[t]
     \centering
         \includegraphics[width=0.99\linewidth]{figures/Framework.png}

\caption{Schematic of the proposed framework. The camera captures images in color, depth and infrared channels (all channel coming from a single commercially available device) and the RetinaNet based multi channel face detector performs simultaneous face detection and presentation attack detection using the composite image from these channels. The framework can be adapted to work with RGB-D combination just by stacking RGB and Depth channels at the input side (to make it work with OpenCV AI Kit (OAK-D) \cite{oak_d} and Kinect \cite{zhang2012microsoft}).}
\label{fig:framework}
\end{figure}

For reliable usage of face recognition systems, they must be equipped with a presentation attack detection module. The PAD module can act before, after, or together with the face recognition module.

Typically face recognition frameworks consists of a preprocessing stage including face detection and alignment, followed by the actual face recognition task. In the proposed PAD framework, we propose to use a multi-channel face detector \textit{as} the PAD system. In this way, the preprocessing stage for the face recognition system itself can perform PAD, thanks to the multi-channel information. Furthermore, the face detectors used typically are also CNN based, so if we can replace this face detector with the proposed module, it should not increase the overall complexity while adding PAD capabilities, simplifying the overall face recognition pipeline by removing the redundancy.

Presentation attack detection is achieved by changing the task of the face detector from \textit{detecting faces} to \textit{detecting bonafide faces}. This is a hard challenge when only RGB channels are used and can result in a lot of detection errors. However, the proposed face detection based PAD framework (Fig. \ref{fig:framework}) leverage the multi-channel information to discriminate between \textit{bonafide} and presentation attacks.

Another advantage of the proposed face detection-based framework is the efficient use of the background. Most of the common PAD frameworks utilize only the facial region and completely ignores the information from the background. Notably, in a multi-channel framework, background could give beneficial negative samples to make the system robust. The proposed framework can be configured based on the availability of data in different ways. It can be trained as a purely one class model when no negative samples are available, i.e, when only \textit{bonafide} samples are available. The \textit{bonafide} face location can be used as positive samples, and the entire background can be used as negative samples. When attacks are available, the face detector can be trained to classify between \textit{bonafide} and the attack classes. It is also possible to train the face detector as a multi-class classifier when we have the labels for different kinds of attack classes. While training the face detector, a lot of negative samples contribute to the loss function. To avoid this, we use the focal loss formulation, which focuses on the hard examples. If the face detection network is trained only on \textit{bonafide} samples, most of the negative samples (patches) would be easy to classify. The attack class images can be used as hard negative samples. This introduces greater flexibility as the face detector can be trained with different configurations, for instance, as a one-class classifier, either with only \textit{bonafide} class, or with the addition of attack classes too in training providing the CNN based face detector with harder examples, or as multi class classifier. This could be useful in the cases where only a limited amount of attacks are available in training.

To the best knowledge of the authors, this is the first work using a multi-channel face detector for the task of face presentation attack detection. The main contributions of this work are listed below.

\begin{itemize}
\item Proposes a novel framework for PAD, using a multi-channel face detector, performing simultaneous face localization and presentation attack detection.
\item The channels used in this work comes from an affordable consumer-grade camera, which makes it easy to deploy this system in real-world conditions.
\item The proposed algorithm can be used as a preprocessing stage in face recognition systems, thus reducing the redundancy of an additional PAD system in the face recognition pipeline.
\end{itemize}

\section{Related work}
Most of the prevailing literature in face presentation attack detection deals with the detection of 2D attacks using visible spectrum images.  Majority of them depend on the quality degradation of the recaptured samples for PAD. Methods such as motion patterns \cite{anjos2011counter}, Local Binary Patterns (LBP) \cite{boulkenafet2015face}, image quality features \cite{galbally2014image}, and image distortion analysis \cite{wen2015face} are examples of feature based PAD methods. Also, there are several CNN based methods achieving state of the art performance \cite{liu2018learning,george-icb-2019, atoum2017face}. While there has been a lot of work in the detection of 2D presentation attacks, the assumptions of quality degradation during recapture does not hold for attacks such as realistic silicone masks and partial attacks. Here we limit the discussion to recent and representative methods that handle a wide variety of 2D and 3D attacks.

With the ever-improving quality of attacks, visible spectrum images alone may not suffice for detecting presentation attacks. This problem becomes more severe when there is a wide variety of possible 2D and 3D presentation attacks possible. Multi-channel and multi-spectral methods have been proposed as a solution for this problem \cite{raghavendra2017extended}, \cite{steiner2016reliable}.

Raghavendra \textit{et al}. \cite{raghavendra2017extended} presented an approach to use complementary information from different channels using a multi-spectral PAD framework. Their method used a fusion of wavelet-based features. Score level fusion achieved better performance as compared to feature fusion in detecting attacks prepared using different kinds of printers. Erdogmus and Marcel \cite{erdogmus2014spoofing} showed that 3D masks could fool the PAD systems easily. By combining the LBP features from color and depth channels, they could achieve good performance in the 3DMAD dataset.

Steiner \textit{et al}. \cite{steiner2016reliable} introduced multi-spectral SWIR image-based PAD method, capturing four different wavelengths - 935\textit{nm}, 1060\textit{nm}, 1300\textit{nm} and 1550\textit{nm}. Their method essentially consisted of a skin level classifier in a predefined Region Of Interest (ROI) where the skin was expected to be present. They trained a pixel-level SVM to classify each pixel as skin or not. The percentage of skin detections in the ROI was used as the PAD score. Their method achieved 99.28\% accuracy in pixel-level skin classification.

In \cite{dhamecha2013disguise}, the authors combined visible and thermal image patches for PAD. The patches used for the subsequent face recognition stage were selected by first classifying the patches as either \textit{bonafide} or attacks.

Bhattacharjee \textit{et al}. \cite{Bhattacharjee:256262} showed the vulnerability of CNN based face recognition systems against 3D masks. They also proposed
simple thermal-based features for PAD. Further, in \cite{bhattacharjee2017you}, they presented preliminary experiments on the use of additional channels for PAD.

In \cite{george_mccnn_tifs2019}, George \textit{et al}. presented a multi-channel face presentation attack detection framework. Their approach consisted of extending a pretrained face recognition network to accept multiple channels and to adapt a minimal number of layers to prevent overfitting. The proposed method achieved an error rate of 0.3 \% in the challenging \textit{WMCA} dataset using color, depth, infrared and thermal channels. Further several other works have reported improved performance with the use of multi-channel methods \cite{heusch2020deep,george-tifs-2020}.

Many multi-channel datasets have been made available for PAD recently. However, the variety of attacks is rather limited in most of the available databases.

Typical PAD frameworks perform a preprocessing stage involving face detection with RGB channel, followed by alignment and binary classification for PAD. The PAD classifier, for most of them, has the form of a binary classifier acting on the detected face regions which might result in poor performance in unseen attack scenarios. Such frameworks ignores the useful information coming from the background regions completely. Moreover, the PAD module adds additional complexity to the face recognition pipeline.

\section{Proposed approach}

In the proposed approach, we use an object detection framework for localizing face and classifying between \textit{bonafide} and \textit{non-face} (attacks) labels. Specifically, we formulate the problem of multi-channel presentation attack detection as a two-class object detection problem. We leverage a single-stage object detection framework for this purpose. The task for the proposed object detector is to detect the presence of `\textit{bonafide}' or `\textit{non-face}' classes. While training, face locations detected by the RGB face detector \cite{zhang2016joint} is used as the ground truth bounding box locations. The class labels of the bounding boxes is considered as `\textit{bonafide}' class for real samples and all other attack classes are grouped as `\textit{non-face}' class. In the evaluation phase, the confidence of the \textit{bonafide} class is used for the scoring. If no object classes (`\textit{bonafide}' or `\textit{non-face}') are detected, then the sample is considered as an attack which was not seen in training. We use the RetinaNet \cite{lin2017focal} as the base architecture for the proposed multi-channel face detector, which obtained state of the art performance in object detection tasks.

\subsection{Preprocessing}
The data used in the network comprises of color, depth, and infrared channels. The color channel is converted to gray-scale. The raw depth and infrared channels available from the hardware are in 16-bit format. We first normalize these channels with Median Absolute Deviation (MAD) based normalization to convert them to 8 bit \cite{nikisins2019domain}.

This stage is detailed as follows, let $\mathbf{I}$ be the image; first,  a vector $\mathbf{v}$ containing non-zero elements of $\mathbf{I}$ is obtained to reduce the effect of background holes and dead pixels. The MAD value is computed as follows,
\begin{equation}
\mathbb{MAD} = median(|\mathbf{v}-median(\mathbf{v})|)
\end{equation}
Once the $\mathbb{MAD}$ is computed, the image can be normalized to 8 bit as:
\begin{equation}
\hat{\mathbf{I}}_{i,j} = \frac{(\mathbf{I}_{i,j} - median(\mathbf{v}) + \sigma \cdot \mathbb{MAD})}{2\cdot\sigma\cdot\mathbb{MAD}} \cdot (2^8-1),
\end{equation}

Values $i$ and $j$ denote the coordinates of the pixels in the image. The $\sigma$ value used in our experiments was four.

Once all the channels are available in 8-bit format, they are concatenated to form the grayscale-depth-infrared composite image, which is used as the input in the subsequent CNN pipeline.
An example of the composite \textit{bonafide} image is shown in Fig. \ref{fig:composite}

\begin{figure}[t]
     \centering
         \includegraphics[width=0.99\linewidth]{figures/Composite.jpg}

\caption{The composite image is formed from gray-scale, depth and infrared channels by stacking the normalized images.}
\label{fig:composite}
\end{figure}

\subsection{Network architecture and loss function}

The objective is to treat PAD problem as an object detection problem. In literature, there are several two-stage detectors that can do object detection by first generating the candidate locations and then classifying the candidate locations to get object labels. However, one stage object detectors are advantageous due to their simple design and faster inference. The RetinaNet \cite{lin2017focal} architecture is a one-stage object detection architecture, which is simpler and quicker at the inference stage. One of the main issues with one stage detectors is the poor performance due to heavy class imbalance. Out of a large number of candidate locations, only a few of them contain objects. This could result in poor models as the loss includes a lot of contributions from the background. In RetinaNet \cite{lin2017focal}, this issue is handled by introducing a new loss function called focal loss. This is done by down-weighting easy examples so that their contribution is less in the loss function. In summary, focal loss naturally favors the training on hard instances.

\subsubsection{Architecture}

We use the standard RetinaNet architecture built from a backbone network as the object detector. In our implementation, we used  Feature Pyramid Network \cite{lin2017feature} with ResNet-18 \cite{he2016deep} backbone. There are two subnetworks; one is tasked with regressing the bounding box and the other one is performing the classification. In our case, the classifier needs to classify the presence of `\textit{bonafide}' or `\textit{non-face}' categories. The task of the regression network is to predict the `\textit{bonafide}' (or `\textit{non-face}' ) bounding box. Instead of RGB images, we use the composite image created in the preprocessing stage as the input to the network.

\subsubsection{Loss function}

As shown earlier, the RetinaNet architecture consists of two subnetworks, one for landmark regression and one for classification. For the object classification, the typically used cross entropy loss (CE) has the following form:

\begin{equation}
{\textrm{CE}(p,y) = \begin{cases} -\log(p) &\text{if $y = 1$}\\
 -\log (1 - p) &\text{otherwise.}\end{cases}}
\end{equation}

In the above equation, $y \in \{\pm1\}$ specifies the label, $p \in [0,1]$ is probability computed by the network for the class with label $y=1$. We define $p_\textrm{t}$ as:

 \eqnnm{pt}{p_\textrm{t}=\begin{cases} p &\text{if $y = 1$}\\ 1 - p &\text{otherwise,}\end{cases}}
so that we can rewrite the loss as  $\textrm{CE}(p,y) = \textrm{CE}(p_\textrm{t}) = - \log (p_\textrm{t})$.

After adding the modulation factor {\small $(1 - p_\textrm{t})^\gamma$}, and  $\alpha$-balancing the expression for focal loss (FL) becomes \cite{lin2017focal}:

 \eqnnm{flalpha}{\textrm{FL}(p_\textrm{t}) = - \alpha_\textrm{t} (1 - p_\textrm{t})^\gamma \log (p_\textrm{t}).}

The classification head of the network is trained with this loss function which reduces the effect of easily classifiable background patches in the loss function.

\subsubsection{Scoring method}

The network described above is trained as an object detector to localize and classify face. Here the notation `\textit{face}' denotes the location of face irrespective of the class label. The possible classes are `\textit{bonafide}' and `\textit{non-face}'.
The `\textit{non-face}' category might contain a wide variety of attacks, and its possible that in some cases, no objects are detected. However, for the PAD system, we need to get a score, which is defined as the probability of a \textit{bonafide} class. The network outputs a classification label and class probability for the detected class, above a detection threshold. Once this information is available, the scoring is performed as follows. If the object class predicted is `\textit{bonafide}' then PAD score is computed as $score=P_{bonafide}$ where $P_{bonafide}$ is the probability of `\textit{bonafide}' class. If the object class returned is `\textit{non-face}' then $score=1-P_{non-face}$ where $P_{non-face}$ is the probability of the `\textit{non-face}' class. In the case where no objects are detected, as the case when an unseen attack is present, the sample is considered as an attack and given a predefined low score (i.e., considered as an attack).

\subsubsection{Implementation details}

We used a standard RetinaNet architecture with ResNet-18 \cite{he2016deep} backbone for our experiments. Pretrained weights of a model trained on ImageNet \cite{deng2009imagenet} was used to initialize the ResNet blocks. To train the face detector, $train$ set of the PAD database was used for training, and the $dev$ set was used for validation.  The model was trained as a standard object detector with two classes (`\textit{bonafide}' and `\textit{non-face}' classes).
The face bounding box used in the training was obtained from MTCNN algorithm \cite{zhang2016joint} on color channel images. Adam Optimizer \cite{kingma2014adam} with a learning rate of $2\times10^{-5}$ was used in the training. The model was trained for 50 epochs with ten frames from each video on a GPU grid. All the layers were adapted during the network training. The model corresponding to minimum validation loss was selected for the evaluation. The implementation was done with PyTorch \cite{paszke2017automatic}.

\section{Experiments and Results}

\subsection{Databases used}

The key requirement of the proposed framework is the use of `multi-channel' information.
We have conducted experiments in the publicly available  \textit{Wide Multi-Channel presentation Attack} (\textit{WMCA}) \cite{george_mccnn_tifs2019} database, which contains a wide variety of 2D and 3D presentation attacks. To the best of our knowledge, there are no other publically available PAD datasets that contain synchronized multi-channel data in a usable format for our task. One other notable multi-channel PAD dataset is CASIA-SURF \cite{zhang2020casia}. However, they only provide cropped and preprocessed images in the public distribution which cannot be used in an object detector framework. The \textit{WMCA} dataset contains a total of \textit{1679} video samples from \textit{72} individuals. Mainly, the database contains four different channels of information recorded simultaneously, namely, color, depth, infrared, and thermal channels, collected using two consumer-grade devices,  Intel\textsuperscript{\textregistered} RealSense\texttrademark SR300 (for color, depth and infrared), and Seek Thermal CompactPRO (for the thermal channel). The database contains different types of 2D and 3D attacks, such as print, replay, funny eyeglasses, fake head, rigid mask, flexible silicone mask, and paper masks. The statistics of attack types is shown in Table \ref{tab:BATL-data}. More details of the database can be found in \cite{george_mccnn_tifs2019}. Even though four different channels are available in the database, we use only three channels coming from the consumer-grade Intel RealSense camera (color, depth, and infrared). This makes the practical deployment more feasible as compared to the method in \cite{george_mccnn_tifs2019}. The samples of RGB images for the attack categories are shown in Fig. \ref{fig:pa_wmca}

\begin{figure}[t]
     \centering
         \includegraphics[width=1\linewidth]{figures/face_pais.png}

\caption{Attack categories in \textit{WMCA}, (image taken from \cite{george_mccnn_tifs2019}) dataset (Color channel), (a): glasses (paper glasses), (b): glasses (funny eyes), (c): print (2D), (d): replay (2D), (e): fake head, (f): rigid mask (Obama mask), (g): rigid mask (transparent plastic mask), (h): rigid mask (custom made), (i): flexible mask (custom made), and (j): paper mask.}
\label{fig:pa_wmca}
\end{figure}

\begin{table}[ht]
\centering
\caption{Number of samples in each presentation attack category in \textit{WMCA} database.}
\label{tab:BATL-data}
\begin{tabular}{lcr}
\toprule
Attack Type                        & Category & \#Videos                \\ \midrule
\textit{Bonafide} & -    & 347                            \\
Print                              & 2D   & 200                            \\
Replay                             & 2D   & 348                            \\
Fake head                          & 3D   & 122                            \\
Rigid mask                         & 3D   & 137                            \\
Flexible mask                      & 3D   & 379                            \\
Paper mask                         & 3D   & 71                             \\
Glasses                            & 3D (Partial)   & 75                    \\ \hline
\textbf{TOTAL}    &      & \textbf{1679} \\ \bottomrule
\end{tabular}
\end{table}

\subsection{Protocols}

We use the same \textit{grandtest} protocol defined in \cite{george_mccnn_tifs2019} for our experiments. However, it is to be noted that we only use three channels out of the four channels available from the database. All the channels we used are coming from the same consumer grade device making it suitable for deployment.

\subsection{Metrics}
For the evaluation of the algorithms, we have used the ISO/IEC 30107-3 metrics \cite{ISO}, Attack Presentation Classification Error Rate (APCER), and Bonafide Presentation Classification Error Rate (BPCER) along with the Average Classification Error Rate (ACER) in the $eval$ set. We compute the threshold in the $dev$ set for a BPCER value of 0.2\%. EPC curves and APCER and BPCER are also calculated. Additionally, the APCER of individual attack types and the ACER-AP which is defined as the average of BPCER and the maximum APCER of the attack types are also reported.

\subsection{Baseline methods}

We have implemented feature-based and CNN based baselines to compare with the proposed method. We used the same channels, i.e., gray-scale, depth, and infrared channels in most of the baselines for a fair comparison (except for the use of color channel in some baselines). We compared the proposed approach with competing multi-channel baselines using the same channels. The description of the baseline methods is given below.

\begin{itemize}
\item \textit{Haralick Fusion}: Here we use an extension of RDWT-Haralick-SVM method in \cite{agarwal2017face}. First, we perform a preprocessing stage, which consists of face detection. After that, the face region in all channels are normalized to the size of $128\times128$ pixels and rotated to make the eye-line horizontal. For each channel, Haralick \cite{haralick1979statistical} features computed from $4 \times 4$ grid are concatenated to form the feature vector. This feature vector was used with logistic regression (LR) for individual channels for the PAD task. A score fusion of all the channels was performed to get the final PAD score.

\item \textit{IQM-LBP fusion}: In this baseline, we followed a similar preprocessing as the previous baseline. After the preprocessing stage, features are extracted from different channels, for example, Image Quality Measures (IQM)~\cite{galbally2014image} were extracted for the RGB channel, and variants of Local Binary Patterns (LBP) features for non-RGB channels. For each channel, a logistic regression model is trained, and score level fusion is performed to obtain the PAD scores.

\item \textit{FASNet (Color)}: This is an RGB only CNN \cite{lucena2017transfer} baseline for the PAD task. Similar to the other baselines, this network also utilizes aligned images from the preprocessing stage.

\item \textit{MCCNN-GDI}: Here we implemented \cite{george_mccnn_tifs2019}, which achieved state of the art performance in the \textit{grandtest} protocol in \textit{WMCA} dataset, when used with gray-scale, depth, infrared and thermal channels. Here, to make it comparable with that of the proposed system, we retrained the system with the channels used in this work, i.e., gray-scale, depth, and infrared channels.

\end{itemize}

Additionally, for the feature-based methods, the best performing method from each channel is also added in the evaluations.

\subsection{Experiments}

\subsubsection{Baseline results}

The results of different baselines systems with color, depth, and infrared channels, as well as the fusion methods, are shown in Table \ref{tab:baseline_results}. For the individual channels, only the result from the best one is reported. Of all the channels present, the infrared channel performs the best with an ACER of 11.0\%. Though fusion improves the overall ROC, the performance in low BPCER regions becomes slightly worse. Out of the baselines, the MCCNN-GDI method is clearly superior, achieving an ACER of 3\%.

\begin{table*}[t]
\centering
\caption{Performance of the baseline systems and the components in \textbf{grandtest} protocol of \textit{WMCA} dataset. The values reported are obtained with a threshold computed for BPCER 0.2\% in $dev$ set.}
\label{tab:baseline_results}
\begin{tabular}{@{}lcc|ccc@{}}
\toprule
\multirow{2}{*}{Method} & \multicolumn{2}{c|}{dev (\%)}     & \multicolumn{3}{c}{test (\%)}                                  \\ \cmidrule(l){2-6}
                        & \multicolumn{1}{c|}{APCER} & ACER & \multicolumn{1}{c|}{APCER} & \multicolumn{1}{c|}{BPCER} & ACER  \\ \midrule

Color (Haralick-LR)         &41.1  &20.6  &42.0     &1.0  &21.5 \\
Depth (Haralick-LR)         &42.2  &21.2  &50.6  &0.2  &25.4 \\
Infrared (Haralick-LR)       &24.9  &12.5  &21.9  &0.0  &11.0 \\ \hline
Haralick Fusion  &29.0 &14.6 &27.0 &0.3 &13.6 \\
IQM-LBP fusion  &41.9 &21.0 & 47.7 & 0.3 &24.0 \\\hline
FASNet (Color) &29.3  &14.7  &24.9  &3.8  &14.3 \\
MCCNN-GDI &2.3 & 1.2 & 6.1 & 0.0 & 3.0 \\ \hline
\textbf{Proposed} &1.2 & 0.7 & 1.9 & 1.0 & \textbf{1.5} \\ \bottomrule
\end{tabular}
\end{table*}

\subsubsection{Results with the proposed framework}

From Table \ref{tab:baseline_results}, it can be seen that the proposed face detection based PAD framework achieves an ACER of 1.5\% outperforming the state of the art baselines. Some of the successful detections are shown in Fig. \ref{fig:successful}.

To identify the limitations of the framework, some of the false detections are also shown in Fig. \ref{fig:failure}. The instances of \textit{bonafide} getting classified as the attack could be due to the missing values in the depth channel. Depth data from the device contains holes where depth is not estimated correctly. This can probably be avoided with a preprocessing stage with hole filling and smoothing of the depth channel.

The same could explain the misclassification of the transparent mask as \textit{bonafide}. Some of the attacks with paper glasses are misclassified since it looks identical to \textit{bonafide} with corrective glasses.

\begin{figure}[h]
     \centering
         \includegraphics[width=1\linewidth]{figures/successful.jpg}

\caption{ Image showing the successful detections, green boxes indicate the \textit{bonafide} detections and red boxes denote attack detections.
First row shows the successful \textit{bonafide} detections and second row shows the successful presentation attack detections.}
\label{fig:successful}
\end{figure}

\begin{figure}[h]
     \centering
         \includegraphics[width=0.8\linewidth]{figures/failures.jpg}

\caption{ Image showing the failure case in detections, green boxes indicate the \textit{bonafide} detections and red boxes denote attack detections. The first row shows the failure cases where \textit{bonafide} are classified as attacks and the second row shows the failure cases where attacks are misclassified as \textit{bonafide}.}
\label{fig:failure}
\end{figure}

\subsubsection{Detailed analysis of the scores}

In this section, we perform a detailed analysis of the scores to get insights into the misdetections and failure cases.
For this analysis, we compute the decision threshold in the development set with BPCER 0.2\% criterion, and we apply this threshold on the evaluation set. After that, we aggregate the APCER for individual attack categories. The results for different attack categories are tabulated in Table \ref{tab:detailed_pa}.

\begin{table}[h]
\centering
\caption{Detailed analysis of the performance for different PA's, at threshold computed for BPCER 0.2\% in $dev$ set.}
\label{tab:detailed_pa}
\begin{tabular}{@{}lccc@{}}
\toprule
Attack Type       & Category & Dev Set & Eval Set \\ \midrule
Prints        &2D& 0.0\%   & 0.0\%    \\
Replay        &2D& 0.0\%   & 0.0\%    \\
Flexible mask &3D & 2.0\%   & 0.8\%    \\
Rigid mask    &3D& 3.7\%   & 3.2\%    \\
Fake head     &3D& 0.0\%   & 0.4\%    \\
Glasses      &3D (Partial)& 2.5\%   & 11.9\%   \\ \hline
APCER-AP         & & 3.7\%   & 11.9\%   \\
BPCER             & & 0.2\%   & 1.0\%    \\\hline
\textbf{ACER-AP}  & &1.9\%    &6.5\% \\\bottomrule
\end{tabular}
\end{table}

From the Table \ref{tab:detailed_pa}, it can be seen that the proposed method achieves perfect performance in 2D attack classification. This is expected since the information from the depth channel alone might suffice for this task. The performance in 3D masks are also satisfactory, considering the very low BPCER values. The glasses attacks are the most challenging ones, and they contribute to most of the misclassifications. The glasses category consists of `Paper glass` (Fig. \ref{fig:pa_wmca} (a)) and `Funny eyes` (Fig. \ref{fig:pa_wmca} (b)) sub categories. Moreover, these are the only partial attacks present in the database. Especially in the channels considered, the funny eyes and the paper glass appear very similar to \textit{bonafide} samples with medical glass. This results in misclassification of these attacks.

The ROC plots provide information about the performance in the evaluation set only. However, in realistic conditions, the decision threshold has to be set a priori. Expected performance curves (EPC) \cite{bengio2005expected} provide an unbiased estimate of the performance of  (two-class) classifiers. As opposed to the ROC curve, where the evaluation is done only on the test set, EPC combines both development and test group to obtain the expected performance. The EPC curve shows the performance changes as a trade-off between False Matching Rate (FMR) and False Non Matching Rate (FNMR). This trade-off is controlled by the parameter \mbox{$0 \leq \alpha \leq1$} as follows:
\begin{equation}
\label{eqn:weighted_error}
WER = \alpha\times(FMR) + (1 - \alpha)\times(FNMR)
\end{equation}
where $WER$ (weighted error rate) is computed as a weighted combination of FMR and FNMR.

For each value of $\alpha$, the score-threshold which minimize the WER is selected based on the development group. The threshold chosen is used to compute the Half Total Error Rate (HTER) in the evaluation set.

Fig. \ref{fig:epc_WMCA} shows the EPC plots of the proposed approach and the baseline methods. From the EPC curves, it can be seen that, in the lower BPCER regions (\mbox{$\alpha \leq0.5$}), the proposed approach outperforms the state of the art methods, indicating the robustness of the approach.

\begin{figure}[ht]
\centering
\includegraphics[width=0.8\linewidth,page=5]{figures/grandtefacedet_evaluate.pdf}
\caption{EPC plots for the reproducible baselines and proposed method in the \textit{eval} set of \textit{grandtest} protocol in \textit{WMCA} dataset }
\label{fig:epc_WMCA}
\end{figure}

\subsection{Discussions}

From the experimental results, it can be seen that the proposed framework for PAD with multi-channel data works better than binary classifier based approaches in the low BPCER range. This framework is ideal for typical usage settings since face detection is a part of the face recognition pipeline. Instead of using conventional face detection methods which detect and pass both \textit{bonafide} and attacks to the subsequent stages, the proposed method uses a multi-channel face detector, specifically trained to localize and distinguish \textit{bonafide} and \textit{attack} ('\textit{non-face}') classes simultaneously. Also, it is worth noting that the hardware used in the proposed pipeline is a single off-the-shelf Intel RealSense SR300 camera. To sum it up, the proposed approach provides face PAD as a face detector itself. Currently, the multi channel presentation attack databases available are collected in controlled conditions. Availability of data collected `in the wild' conditions could make it possible to train the face detector using just the \textit{bonafide} samples. The harder negatives from the background might improve the performance of the detector greatly.

\section{Conclusions}

The proposed method shows a simple yet efficient method for PAD using multi-channel information. The method does not add additional complexity to a face recognition pipeline
since it can be a drop-in replacement for the face detector in the preprocessing stage. Even though the method is multichannel, the channels used are readily available from cheap consumer devices which makes it feasible for deployment. The Intel RealSense family of devices, Kinect and the upcoming OpenCV AI Kit (OAK-D) \cite{oak_d} are suitable candidates for such system as this model can be replaced with RGB-D channels as well.

Addition of a PAD module is essential to secure face recognition systems against spoofing attempts. Typical PAD modules adds another layer of computational complexity to the face recognition pipeline. In this work, we propose a simple yet effective solution to add PAD capabilities to a face recognition system without adding any additional computational complexity. This is done by swapping the face detector in typical FR pipelines with the proposed multichannel face detector. The (multi-channel) face detector is used \textit{as} a presentation attack detector, thereby reducing the redundancy in the face recognition pipeline. To our best knowledge, this is the first method to use multi-channel face detection \textit{as} a PAD system. Specifically, a multi-channel face detector is trained to localize and classify \textit{bonafide} faces and presentation attacks. The proposed method utilizes color, depth, and infrared channels available from a commercially available sensors for PAD. The proposed method was compared with feature-based methods and state of the art CNN based methods for comparison in the \textit{WMCA} dataset and was found to outperform the state of the art method while obviating the requirement of additional face detection stage.

\section*{Acknowledgment}

Part of this research is based upon work supported by the Office of the
Director of National Intelligence (ODNI), Intelligence Advanced Research
Projects Activity (IARPA), via IARPA R\&D Contract No. 2017-17020200005.
The views and conclusions contained herein are those of the authors and
should not be interpreted as necessarily representing the official
policies or endorsements, either expressed or implied, of the ODNI,
IARPA, or the U.S. Government. The U.S. Government is authorized to
reproduce and distribute reprints for Governmental purposes
notwithstanding any copyright annotation thereon.

{\small
\bibliographystyle{IEEEtran}