\section{Introduction}

The 8.2-m Subaru telescope has  been running
a 300-night Strategic Survey Program (SSP)
over 5 years since 2014 March\footnote{http://www.naoj.org/Projects/HSC/surveyplan.html},
in order to elucidate the mystery of dark matter and dark energy
as well as the evolution of galaxies.
The survey utilizes Hyper Suprime-Cam (HSC; \cite{Miyazaki+2012})
with a wide field of view of 1.77 square degrees.
In 2016--18, we observe two ultra deep fields, COSMOS and SXDS,
for six-month each.
The observations will be performed around the new moon
in each month, with a typical cadence of 3--4 days.
Note that the HSC is not installed on Subaru for two weeks
around the full moon.
Among other transient surveys,
the Subaru HSC/SSP survey is the deepest for this survey area
($1.77 \times 2$ square degrees), and thus,
will provide the unique dataset for transients.
For example, the HSC/SSP survey will
triple the number of type-Ia supernovae (SNe) beyond redshift $z>1$
and will also discover a few tens of superluminous SNe at $z>1$.

Difference imaging is the standard method to search for optical transient
objects, and so we use it in this study.
We define a transient object as the one that  appears only
in the image at the later epoch (newer image),
but not in the one at the earlier epoch (reference image or template image), out of the two images taken at different epochs.
After the standard data reduction is  made,
two images are astrometrically aligned, and
the reference image is subtracted from the newer image
by matching point spread functions (PSFs) of the two images.
The source-finding algorithm is applied to the difference image, and  the detected sources are
the candidates of transients
(section \ref{sec: Data Analysis and Feature Extraction}).

\begin{figure}
 \begin{center}
   \includegraphics[width=7cm, angle=0]{bogus_real_sample.eps}
 \end{center}
\caption{Examples of real and bogus objects obtained  with Subaru-HSC.
The left,  middle, and right columns show
 the reference, new, and difference images, respectively.
The first, second, and third rows show
the cosmic ray (a--c), ghost near a bright star (d--f), and
inaccurate image convolution or astrometric alignment (g--i),
respectively.
The bottom row shows a real transient located in a galaxy (j--l).
}\label{fig:bogus_sample}
\end{figure}

In an ideal situation, all the sources detected in difference images
 would be  transient/variable astronomical sources, such as
supernovae, variable stars, moving objects, and so on.
In reality, however, they also include artifacts
(see panels a--i in figure \ref{fig:bogus_sample}),
such as cosmic-ray events, spikes around bright stars,
and residuals related to inaccurate image convolution or astrometric alignment.
These artifacts are  present in every
optical survey project \citep{Bailey+2007, Bloom+2012, Brink+2013}.
Hereafter, we call them ``bogus'' \citep{Bloom+2012}.

In the HSC/SSP survey, not only a few hundred transients, including SNe,
but also $\sim 10^5 - 10^6$ bogus objects are expected to be detected each night.
After the scheduled 300 nights, the number of candidates of transients,
 real and bogus combined, will reach $\sim 10^8$,
which is  well qualified as Big Data.
We need to filter out bogus objects  to select SNe and other real transients.
 Processing of filtering must be  performed swiftly in order to increase
 the chances of new findings in an early phase of
transient phenomena.

 The primary method to  distinguish transient/bogus objects is, traditionally, visual inspection by human checkers,
as many surveys have  been adopting.
However, the expected size of our data  is so big that human checkers would
not be able to go through all the data in a reasonable time.
We have decided to introduce
machine-learning (ML) techniques to select real transients.

In the filtering process,
we should not miss  real objects,
while a vast number of bogus objects are filtered out.
Throughout the development,
we try to minimize the false positive rate (FPR),
while we maintain the true positive rate (TPR) of 90 \% or larger
(namely the false negative rate, FNR $< 10\%$).

We performed two HSC observations in 2015 May and August
\citep{tominaga15atel1, tominaga15atel2}, aiming to detect short transients
with a time scale of a few hours to a few days
(e.g. \citet{Tanaka+2016}; \citet{Morokuma+2016}).
An example of such short transients is an optical flash at the time of
shock breakout of a supernova,
of which the time-variance would be detectable
during an observation for a single night.
We use three kinds of ML methods,
AUC Boosting, Random Forests (RF), and Deep Neural Network (DNN), both individually and in a combined way.
We verify the performance of these machines
by making the receiver operation characteristic (ROC) curves.

Conditions of observations (the noise and seeing)
vary every night.  Hence, the ML classifiers must be robust
against the change of environment.
We use normalized features to reduce the influence of the variation.
To validate the performance of our method, we show the
ROC curves of our machine trained
with the data on one-night observation,
applied to the data on the other night observation.
The results show the proposed classifier is robust.

In optical surveys, ML techniques
have been introduced by \citet{Bailey+2007}
for the data from the Nearby Supernova Factory.
They applied Boosted Decision Trees, RF, and
Support Vector Machines (SVM) and succeeded  in
reducing the number of bogus candidates by a factor of ten.
The Palomar Transient Factory team \citep{Bloom+2012}
used RF, and  achieved the TPR of 92.3\%
at the FPR of 1\%
\citep{Brink+2013}.
The Pan-STARRS1 Medium Deep Survey used
Artificial Neural Network, SVM and RF,  and  achieved
TPR of 90\% at FPR of 1\% \citep{Wright+2015}.
\citet{Goldstein+2015} applied the RF for Dark Energy Survey Supernova program (DES-SN),
and reduced the number of transient candidates  by a factor of 13.4, which were then fed to human scanning.
\citet{du Buisson+2015} applied the RF, $k$-nearest neighbor, and the SkyNet artificial
neural net algorithm, using features trained from eigen-image analysis for
the Sloan Digital Sky Survey supernova survey.

The rest of this paper is  structured as follows.
In section~\ref{sec: Data Analysis and Feature Extraction},
we explain the HSC data reduction and feature
extractions.
In section~\ref{sec: ML}, we introduce three machine learning methods we used.
In section~\ref{sec: experiment}, the applications to the actual Subaru data are presented, and then
we discuss the result of real vs bogus segregation in section~\ref{sec: Discussion and Conclusion}.

Prior to the forthcoming HSC/SSP Transient survey,
this paper  will provide a `path finder'
to identify real  astronomical
objects and to demonstrate the power of  machine learning.

\section{Data Analysis and Feature Extraction}
\label{sec: Data Analysis and Feature Extraction}

\begin{longtable}{ll}
\caption{List of the features}
\label{table:feature}
\hline\hline
\multicolumn{1}{l}{Feature variable} &
\multicolumn{1}{l}{Description} \\
\hline
\endfirsthead
mag     & Magnitude \\
magerr  & Error of magnitude \\
elongation\_norm & Elongation normalized by nearby stars \\
fwhm\_norm     & FWHM normalized by nearby stars \\
significance\_abs & Significance obtained  with the PSF fit \\
residual       & Residual of PSF fit \\
psffit\_sigma\_ratio & See text \\
psffit\_peak\_ratio  & See text \\
frac\_det & See text \\
density & Number of objects around the target within a square with 120 $\times$ 120 pixels
centered  at the target \\
density\_good & ``Density'' after weak screening \\
bapsf & Elongation of nearby stars \\
psf & FWHM of nearby stars \\
\hline
\end{longtable}

We describe the flow of HSC data reduction and feature extraction for machine learning.
The pipeline processing,
using the on-site data analysis system \citep{furusawa11},
and  then transient finding \citep{tominaga15}, are performed immediately after the data acquisition.

The HSC data are reduced  with the HSC pipeline
(version 3.6.1), which  has been developed based on the LSST pipeline
\citep{ivezic08,axelrod10}.
First, the pipeline  performs the standard reduction,
such as bias subtraction and flat fielding.
Then, astrometric and photometric calibrations are made,
using the Sloan Digital Sky Survey DR8 data \citep{aihara11}.
Finally, mosaic solution is derived and images are warped
to align astrometrically
in preparation for co-adding (if needed) and image subtraction
for difference-imaging at the next stage.

For image subtraction we use the HSC pipeline.
The algorithm to match the PSFs of the two images
is the same  as that by \citet{alard98} and \citet{alard00}
\footnote{http://www2.iap.fr/users/alard/package.html},
 and adopts a position-dependent convolution kernel.
The optimal convolution kernel is derived so that
the difference between convoluted PSFs becomes the smallest.
The algorithm has been implemented also in the HOTPANTS package
\footnote{http://www.astro.washington.edu/users/becker/v2.0/hotpants.html}.

Source detection is performed in  difference images,
using the HSC pipeline instead of the standard tool for it,
SExtractor \citep{bertin96}, because
the former outperforms the latter
 in reducing  bogus detections.\newline
For each source, the features listed
in table \ref{table:feature}  are computed as follows. We fit the  image with the PSF of two-dimensional Gaussian
and subtract the best-fit  PSF from  it.
Then, we measure the residuals within 2 $\times$ FWHM radius
($\sigma_{\rm on}$) and in the surrounding region (3--4 $\times$ FWHM, $\sigma_{\rm off}$).
We define ``psffit$\_$sigma$\_$ratio'' $= \sigma_{\rm on}/\sigma_{\rm off}$;
 it is expected to be much larger than unity
when the shape of  a source is  very different from PSF.
We also define ``psffit$\_$peak$\_$ratio'',
the ratio of the actual peak of the source
to the peak of the best-fit Gaussian in order to quantify the degree of deviation of
 the source-image shape  from the PSF.

To  obtain features sensitive to mis-alignment,
we count the number of positive pixels within $3 \times$ FWHM,
 and define ``frac$\_$posi'' as  their fraction.
This feature frac$\_$posi should  be close to unity
for   good detection with a high S/N ratio.
Similarly, we count the number of negative pixels and define ``frac$\_$nega'' for its ratio to detect ``sources'' with negative counts.
To use both the cases in a consistent manner, we define
``frac$\_$det'', which is the same  as ``frac$\_$posi'' for
the positive detection  but ($1 -$ ``frac$\_$nega'') for the negative detection.
Following \citet{Bloom+2012} and \citet{Brink+2013},
we also count the number of the detected sources inside
a 120 $\times$ 120 pixel box centered  at the source (``density'').

In the HSC pipeline,
PSF fitting for the detected sources  is
performed  with the position-dependent PSF.
We also compute the significance of  detection (significance$\_$abs)
by comparing the fitted PSF with the noise level around   a source.

In total, 13 features  are used for machine learning (table \ref{table:feature}).

\section{Methods of Machine Learning for Real-Bogus Separation}
\label{sec: ML}

\subsection{AUC Boosting}
\label{sec: AUC boosting}

Boosting is a method to classify the data
by majority-voting of weak classifiers
\citep{Hastie-Tibshirani-Friedman_2009}.
Among various boosting methods,
we employed the AUC boosting
developed by \citet{Komori_2011},
which is trained to maximize the empirical
area under the ROC curve (AUC).
The AUC boosting classifies objects according to
a score function, where the objects with larger scores are regarded as real.
This method uses only one hyper-parameter $\lambda$
to control the smoothness of the score function,
which  is optimized through cross-validation.

Once the machine  has been trained,
classification for  a new set
of feature variables will be fast.
Hence, when the trained machine is installed in
the pipeline process of HSC,
real and bogus are classified quickly.
Note that the computation speed for the classification
is fast for the following two methods, too.

\subsection{Random Forest}

Random Forest (RF; \cite{Breiman_2001})
is an ensemble-learning method using an ensemble of decision trees.
Each decision tree is trained with a subset of training data.
Classification is based on majority-voting of the decision trees.
For a training data set $X=\{{\bf x}_1, {\bf x}_2,\ldots,{\bf x}_n\}$
with labels $Y=\{y_1, y_2,\ldots,y_n\}$,
($y=0$ or 1, corresponding to bogus or real, respectively),
each of $B$ training sets  $X_b$ ($b=1,...,B$)
is selected by bootstrap sampling (i.e., randomly drawn with replacement
from $X$). Then, the $b$-th decision tree $f_b$ is trained
with $X_b$ with the corresponding label set $Y_b$.
When  each decision tree is trained, a subset of the features
is also randomly selected.

RF is considered to be robust  against outliers and noise
because its majority-voting of low-correlated classifiers
decreases the variance of the model.
The meta parameters of RF are the number of trees, $B$,
and the number of subsampled features, $m_f$.
Following the standard practice, we use cross-validation
 to evaluate the performance and  determine the appropriate $B$ for our observation data,
and set $m_f$  $=$ $\sqrt{p}$, where $p$ is the number of all the features.
We employ the RF implementation in scikit-learn\footnote{An open-source machine learning library for python, http://scikitlearn.org}.
The standard setting for the classification task is used.

\subsection{Deep Neural Network}

Deep Neural Network (DNN) or Deep Learning \citep{LeCun15}
is the current state-of-the-art technique,
which achieves the best performance in speech and image recognitions.
The network consists of multiple layers
of neurons with directional connections.
The neural network is trained to tune the connecting weights
between neurons and parameters of emission functions,
so that it approximates a mapping
from the input observations (feature vector; ${\bf x} \in {\bf R}^p$)
fed into the input layer
to the output observations emitted from the output layer
(bogus or real; $y = 0, 1$).
Efficient stochastic optimization algorithms
are used to tune millions of weights and parameters.
In our case, we use a fully-connected feed-forward network.

We use the Chainer library\footnote{
It is provided by Preferred Network Inc., http://chainer.org/}.
For the parameter estimation, we use
the stochastic gradient descent (SGD) methods.
The SGD computes the noisy gradient of the objective function
based on the mini-batch, a subset with  $M (\ll N)$,
out of   $N$ samples.
DNNs have many meta parameters to be tuned.
We choose them via preliminary cross-validations.

\section{Experiment}
\label{sec: experiment}

We performed HSC observations on 2015 May 24 UT and August 19 UT with a high cadence
\citep{tominaga15atel1, tominaga15atel2},
aiming to detect short-timescale transients, of which the time-variance would be detectable
during an observation for a single night.
From the observational data of May, optical transients were detected and screened
by conventional visual inspection, and 48 definitely real transient objects, most likely supernovae,
were identified (see the next sub-section).
Note that they were not used for training machines but for performance validation of the machines
(section \ref{sec: result sel SN}).
We created data set for training based on the May data (section \ref{sec: Artificial real objects})
and trained machines with them. The machines were applied to the August data
and transients events were identified (sections \ref{sec: result sel SN} and \ref{sec: quick sel SN}).

To validate the robustness of the machines against variation of environment,
we also did the reverse; i.e., we trained the machines with the August data
and  applied them to  the May data and their
 48  definite transients.
 In the training, we added artificial  transient objects
because the number of real transient objects in the data is  much smaller than
 that of the bogus. We explain how we created the training data
in section \ref{sec: Artificial real objects}.

\subsection{Dataset on 2015 May}

In this observing run,
eight HSC field-of-views (about 14.2 deg$^2$)
were repeatedly observed with  roughly 1-hr interval
in $g$-band (3--4 visits) and $r$-band (1 visit).
One visit consisted of 3 frames of 2-min exposure
(3 images are co-added for each epoch).

To detect short-timescale transients,
image subtraction was performed by using the  first-visit data as the references.
Therefore, source detection in  difference images
was sensitive to intranight variability.
 From the difference images, 54,672 sources were detected in total.
Features listed in table \ref{table:feature}
were computed for these sources.
The average 5-sigma detection limit was $24.7$ mag.

Furthermore, we performed independent image subtraction,
using the reference data taken on 2014 July 2 and 3 UT.
The difference images between the reference
and this observing run revealed many transients.
Among them, 48 definite transients were
discovered with visual inspection
and were reported as supernova candidates
by \citet{tominaga15atel1}.
We used the 48 supernova candidates as real transient objects
in section \ref{sec: result sel SN}.

\subsection{Dataset on 2015 Aug}

In this observing run, the weather condition was  poor and
only 3-hr data were taken under poor seeing of 1.1--1.5 arcsec.
Then, it was not possible to take sufficient baseline
for the image subtraction within the same night.
Hence, we  made difference images
for three HSC field-of-views (5.3 deg$^2$)
with the reference images taken in 2014 Jul.
 From the difference images, 45,019 sources were detected in total.
Features listed in table \ref{table:feature}
were computed for those sources.
The average 5-sigma detection limit was $25.3$ mag.

\subsection{Artificial real objects}
\label{sec: Artificial real objects}

\begin{figure}
 \begin{center}
   \includegraphics[width=8cm, angle=0]{mag_dist.eps}
 \end{center}
\caption{Distribution of magnitudes of bogus objects (left) and
artificial ``real'' objects (right)  in the training data
made from the May observation.
}\label{fig:mag_dist_artificial_bogus}
\end{figure}

The small number of real sources, as well as
imbalance between the real and bogus, is a major obstacle
in effective machine-learning.
To improve the performance, we generated
artificial transient objects to  use
as the real sample  in  training  machines.

 The procedure is as follows.
We constructed spatially-varying PSFs based on the detected stars,
using the HSC pipeline, and
generated a number of artificial sources from them
with different brightness  at random positions.
Magnitude ($m$) distribution of artificial sources were  set
to follow  $N(m) = 10^{0.6 m}$, which is expected
when the density of sources with the same
luminosity is constant in the universe.
The faintest magnitude was set  at 27.0 mag, and
1,000 artificial sources were generated
per  CCD chip. These artificial sources were
added to the observed images
to mimic the actual observing conditions.

We made the training sample containing
33,742 real and 25,468 bogus sources.
Figure \ref{fig:mag_dist_artificial_bogus}
shows the distribution of magnitudes
 of bogus and artificial ``real'' objects.
This has made up the number of ``real'' objects
comparable to that of the bogus, and accordingly
 enabled training with a sufficient number of ``real'' objects.
The data with artificial sources were reduced and
features were also extracted in the same manner
as for the real sources.  Note that this strategy was  taken also by \citet{Wright+2015}.

\subsection{Training of Machines}
\subsubsection{AUC Boosting}
\label{AUC Boosting}

\begin{figure*}
 \begin{center}
   \includegraphics[width=6cm, angle=-90]{task1_fake31822_roc3.ps}
 \end{center}
\caption{ROC curves of the three machines:
  AUC Boosting (left), Random Forest (middle), and
  Deep Neural Network (right), trained  with the artificial data,
   based on the May observation.
  Horizontal and vertical axes are FPR and TPR, respectively.
  The red solid and black dotted lines
   show  the curves  for the optical magnitude ranges of  22.0--25.0,
  and  25.0--26.0 mag, respectively.
  Insets show the magnification of the lines  for  22.0--25.0 mag.
}\label{fig:roc3}
\end{figure*}

In the training step, we first determined
the optimal hyper-parameter $\lambda$
by performing cross-validation as follows.
We made 30 partitions of randomly sampled data,
each of which contained the training, validation and test data
with the same size.
For each partition of data, we searched for the optimal $\lambda$
parameter,  evaluating the performance of the machine
with the AUC value for the corresponding validation data.
Next,  fixing the $\lambda$ parameter to the derived optimal,
we trained the machine using the training and validation data,
and evaluated the performance  with the test data.

Figure \ref{fig:roc3} (left\ panel) shows the result: the average of
30 ROC curves  of two magnitude slices for the data of May.
We obtained the average FPR of 3.0\% at the point of
the TPR of 90\% in the magnitude range of 22.0--25.0 mag.

\subsubsection{Random Forest}

We determined the meta parameter $B=100$
by the cross-validation in the same way as
in the case of AUC boosting (previous subsection). We chose $m_f = 4$, based on the standard
setting for classification task.
The ROC curves in the magnitude slices are shown
in figure~\ref{fig:roc3}~(middle panel).
We obtained the average FPR of 0.95\% at the point of
the TPR of 90\% for the same range
as in AUC Boosting.

\subsubsection{Deep Neural Network}
\label{sec:DNN May}

We chose the meta-parameters by preliminary cross validations, as follows.
The number of hidden layers was $3$, and thus our DNN was $5$-layered.
All emission functions of the hidden units were
ReLU (the rectified linear unit).
The emission function of the output layer was the soft-max function.
We used a classical momentum SGD method to optimize the network
 with the mini-batch size $M=4000$.
The number of maximum iterations was $20000$,
but the early stopping rule was employed,
following the standard practice of DNNs.

For the cross-validation,
we split the sample equally in data size into
training, development, and test dataset,
keeping the ratio between the real and bogus objects,
and prepared $30$ partitions.
For each partition, we performed $30$ repetitions of
training and evaluation to search for the best combination of meta-parameters
 of the number of neural units in hidden layers and the step size of SGD, among exhaustive combinations of them.
Then, we obtained the average test score
of these repetitions.

We obtained the average FPR of 0.85\% at the point of
the TPR of 90\% for the same range as in AUC Boosting.
The ROC curves  are shown in figure \ref{fig:roc3} (right\ panel).

\subsection{Combined machine}

\begin{figure*}
 \begin{center}
   \includegraphics[width=6cm, angle=-90]{task1_fake31822_merge.ps}
 \end{center}
\caption{ROC curves of the three ML machines and the combinations of them for the data of 22.0 -- 25.0 mag on 2015 May.
Those of AUC boosting, RF, and DNN are shown in blue, cyan, and magenta dotted curves,
respectively.
Those of ``unanimous voting'', ``majority voting'', and ``safeguard minority opinions''
are shown in black, red, and green solid curves, respectively.
}\label{fig:roc_all}
\end{figure*}

Finally, we combined the results
of the three machines to have a single decision.
There are logically three ways to combine
three machines  indiscriminately (A, B, and C):  ``unanimous voting'', ``majority voting'',
and ``safeguard minority opinions'',  which are basically ``A \& B \& C'',
``(A \& B) $|$ (B \& C) $|$ (C \& A)'',
and ``A $|$ B $|$ C'', respectively, where ``\&'' and ``$|$''  denote logical ``AND'' and ``OR'', respectively.
Figure \ref{fig:roc_all} summarizes their ROC curves, together with the ones with individual machines.
We found the FPR of the combined machines at the TPR of 90\%  to be
1.0\%, 1.0\%, and 2.1\% for  the above-mentioned three combinations, respectively.
Among them, the performance of
``safeguard minority opinions'' is not good,
 whereas the performances of ``unanimous voting'' and ``majority voting'' are  good.
We have decided to use ``majority voting'', because
the ``unanimous voting'' would give too tight constraints.

\subsection{Robustness for variation of environment}

\begin{figure*}
 \begin{center}
   \includegraphics[width=6cm, angle=-90]{cross_fake31822_22.0-25.0_fake38744_all_21.7-24.7_2and.ps}
 \end{center}
\caption{Combined ROC curves as ``majority voting'' of the three ML machines.
Left: ROC curve made by evaluating the data of 22.0 -- 25.0 mag on 2015 May,
using the machine trained  with the same data (black dot) and the data on 2015 Aug (red solid).
Right: ROC curve made by evaluating the data of 21.7 -- 24.7 mag on 2015 August,
using the machine trained  with the same data (black dot) and the data on 2015 May (red solid).
}\label{fig:roc_all_cross}
\end{figure*}

To validate that our machines  are robust against  variations of
environment, we made the following ROC curves.
Figure \ref{fig:roc_all_cross} (left) shows two ROC curves
of the combined ``majority voting'' machine
applied to the artificial data on May;
the black dotted curve shows the result of the machine trained  with the same data,
while the red solid curve is the result of the machine trained  with artificial data on August.
Figure \ref{fig:roc_all_cross} (right) shows the ROC curves of the data  for which May and August are swapped.
In both the panels, the red solid curves are slightly worse than the black dotted curves,
but the differences are small.
We conclude that the machines trained  with the  adopted method  are robust against  changes of
the conditions.

\subsection{Results of selection of supernovae}
\label{sec: result sel SN}

\begin{longtable}{lr}
\caption{Number of selected objects  in the 2015 August data}
\label{table:number of selected objects}
\hline\hline
\multicolumn{1}{c}{Selection} &
\multicolumn{1}{c}{\# of objects} \\
\hline
\endfirsthead
Total         &  45,019 \\ \hline
AUC Boosting  &  21,487 \\
Random Forest &  16,307 \\
Deep Neural Network  &   11,645 \\ \hline
Two or three machines
              &  16,888 \\
All machines  &   8,514 \\ \hline
\end{longtable}

Table \ref{table:number of selected objects}  summarizes
the numbers of selected objects obtained  with the three individual
machines and with the combined  machines, where
 the thresholds are set at the points corresponding to the TPR of 90\%.
 In order not to miss the real objects, we adopted
the ``majority voting'' for the combined machines,  and it reduced the number of objects to $16,888$.

Among 48 supernovae obtained  with the May observation,
26 are in the magnitude range of 22.0 -- 25.0 mag.
We applied the combined machine of ``majority voting'' with
the threshold of TPR of 90\% at this magnitude range.
Then, 22 sources, or 85\% of them, were  accepted.

\subsection{Quick selection of supernovae}
\label{sec: quick sel SN}

Before the August observation, we had installed the trained machines
in order to run transient search  immediately after the observation
(see section \ref{sec: result sel SN}).
We then visually inspected every object detected, selected ten clean sources as  definite candidates of supernovae,
and reported the list to Astronomers' Telegram \citep{tominaga15atel2}.
This  whole process was  carried out within the same day as the observation.

\section{Discussion and Conclusion}
\label{sec: Discussion and Conclusion}

\begin{figure*}
 \begin{center}
   \includegraphics[width=16cm, angle=0]{score_scat3_sel10.eps}
 \end{center}
\caption{Cross-relations of scores obtained  with the three machines
  for  the entire sample (black) and the selected sample of
  ten candidates of supernovae
  reported to Astronomers' Telegram \citep{tominaga15atel2} (red)
  for the observation data on 2015 August\ 19.
  The left,  middle, and right panels  plot
  those of Random Forest v.s. AUC Boosting,
  AUC Boosting v.s. Deep Neural Network,
  and Deep Neural Network v.s. Random Forest, respectively.
  Horizontal and vertical dashed lines show the thresholds
  corresponding to the TPR of 90\%.}
\label{fig:score_scat3_sel10}
\end{figure*}

We performed Subaru/HSC observations in 2015 May and August, made transient searches,
and  applied machine-learning techniques to the result to reduce the
bogus transient objects.
We  have developed real-bogus classifiers as the core function for it, using
the three machine-learning methods of AUC Boosting,
Random Forests, and Deep Neural Network, and
 then made the combined classifier as  their majority voting.
We have installed  our machines
in the analysis pipeline of the HSC,
and  successfully found real supernovae within the same day as
the observation, demonstrating the power of our method.
Now, we have  completed the preparation for
the forthcoming HSC/SSP transient survey observation.

In training our machines, we used artificial objects,
because  the data are  highly imbalanced between real/bogus objects,
and it was found to be crucial to make
good machines efficiently.
Although the HSC survey data are technically more difficult to deal with than other survey data,  given that the HSC
survey is deeper than other surveys, we have achieved  the results  comparable to  the similar studies in other surveys.

The cross-relations between three machines
in figure \ref{fig:score_scat3_sel10} show that
 none of the three is significantly better than the other two.
Therefore, combining the machines is beneficial.
We used a moderate selection with all the combinations
of multiple machines, ``majority voting'',
to avoid missing some real objects.
For the machine, robustness against variation of environment was confirmed.

In this paper we have focused on real-bogus separation with
machine-learning methods,
and  have demonstrated that they were indeed  useful.  Their use in extracting scientific results
out of Big-data in astronomy is promising.
As the next step, we will use machine-learning  for classification of types of transients,
by combining timing and color information,
as well as the shape of objects.

\begin{ack}
We thank  the Subaru Hyper Suprime-Cam team.
This work is supported by Core Research for Evolutionary Science
and Technology (CREST), Japan Science and Technology Agency (JST).
It is also supported by the research grant program
of Toyota foundation (D11-R-0830)
and was in part supported by
Grants-in-Aid for Scientific Research of JSPS (15H02075),
MEXT (15H00788),
and the World Premier International Research Center Initiative, MEXT, Japan.
This paper makes use of software developed for the LSST.
We thank the LSST Project for  allowing their code available
as free software at http://dm.lsstcorp.org.
\end{ack}


\section{Introduction}
\label{introduction}

The discovery of the time-reversal invariant topological insulators~\cite{hasan2010,qi2011} greatly inspired the study of symmetry-protected Dirac semimetals (DSMs)~\cite{young2012,wangzj2012,wangzj2013,young2014,yang2014,liu2014a,liu2014b,borisenko2014,xiong2015,wieder2016,bradlyn2016,watanabe2016}. In a DSM, the low-energy physics is well described by pseudorelativistic Dirac fermions with linear energy dispersions along all momentum directions around Dirac points (DPs) in the Brillouin zone (BZ). To guarantee the robustness of three-dimensional (3D) DPs, the crystalline symmetry protection is a necessitate~\cite{yang2014}, similar to the 2D DPs in graphene~\cite{geim2009}. The 3D DSM materials~\cite{young2012,young2014,wangzj2012,wangzj2013,yang2014,liu2014a,liu2014b,borisenko2014,xiong2015} have both the time-reversal symmetry $\mathcal{T}$ and inversion symmetry $\mathcal{P}$. The $\mathcal{T}$ or $\mathcal{P}$ breaking in general leads to Weyl semimetals~\cite{murakami2007,wan2011,xu2011,burkov2011,liu2014,huang2015,weng2015,xu2015,lv2015,bernevig2015,sun2015,ruan2016}, where the DPs split into Weyl points. However, under certain conditions, such 3D DPs still remain stable with $\mathcal{T}$-broken~\cite{tang2016}. This motivates us to ask whether stable DPs could exist in 2D in the absence of $\mathcal{T}$ and $\mathcal{P}$. The positive answer to this question will lead to many interesting physical phenomena and novel topological states, as expected from the great success in the field of graphene.

Recently, the concept of 2D DSMs in the presence of spin-orbit coupling (SOC) has been introduced by Young and Kane~\cite{young2015}, where the nonsymmorphic space group symmetries play a key role~\cite{young2015,wieder2016b,chen2014,venderbos2016,zhao2016}. In this paper we study the 2D Dirac fermions in antiferromagnetic (AFM) semimetals. A general framework is proposed to classify stable 2D DSMs in SOC systems with $\mathcal{PT}$ symmetry. In sharp contrast to $\mathcal{T}$-invariant 3D DSMs which fall into two distinct classes~\cite{yang2014}, DSMs in 2D with $\mathcal{PT}$ symmetry only have one class. The DPs of 2D DSMs reside at the BZ boundary which are protected by the nonsymmorphic space group symmetries. We further provide a tight-binding model in AFM semimetals which supports symmetry-protected 2D DPs and symmetry breaking in such systems leads to exotic chiral topological states in 2D. We conclude with a brief discussion on the possible material venues for such phases.

The 2D DP with fourfold degeneracy is generated when two doubly-degenerate energy bands avoidably cross. The most common symmetry for crystalline solids to have double degeneracy of bands is $\mathcal{PT}$ symmetry satisfying $(\mathcal{PT})^2=-1$, which includes two cases with both $\mathcal{T}$, $\mathcal{P}$ conserved and $\mathcal{T}$, $\mathcal{P}$ broken. Under this condition, the low-energy physics can be described by the minimal four-band effective Hamiltonian,
\begin{equation}\label{minimal}
\mathcal{H}(\mathbf{k}) = \epsilon_0(\mathbf{k})+\sum\limits_{a=1}^5d_a(\mathbf{k})\Gamma_a,
\end{equation}
where $\Gamma_a$ are $4\times4$ Dirac $\Gamma$ matrices satisfying $\{\Gamma_a,\Gamma_b\}=2\delta_{ab}$, and the specific representation of $\Gamma_a$ depends on the crystalline symmetry. The particle-hole asymmetry term $\epsilon_0(\mathbf{k})$ is neglected for simplicity, $d_a(\mathbf{k})$ ($a=1,\ldots,5$) are real functions of $\mathbf{k}$, and $\mathbf{k}=(k_x,k_y)$. The energy spectrum is
\begin{equation}
E_\pm(\mathbf{k})=\pm\sqrt{\sum_{a=1}^5d_a^2(\mathbf{k})}.
\end{equation}
The DPs can be generated only when all $d_a(\mathbf{k})=0$ for certain $\mathbf{k}$ in the BZ. A simple case for this condition exists at the quantum critical point of the phase transition between a quantum spin Hall (QSH) and a normal insulator (NI) in the presence of $\mathcal{T}$ and $\mathcal{P}$~\cite{murakami2007,buttner2011}. Such DPs are accidental degeneracies and not robust against perturbations. Therefore, crystalline symmetry is needed to protect the 2D DPs.

The organization of this paper is as follows. After this introductory section, Sec.~\ref{theory} describes the general framework to classify stable 2D DSM with $\mathcal{PT}$ symmetry. Section~\ref{model}
presents a magnetic tight-binding model which supports symmetry-protected 2D DPs, and the symmetry-breaking phases. Section IV presents discussion on the possible material venues for 2D AFM DSMs, and concludes this paper. Some auxiliary materials are relegated
to Appendixes.

\section{General Theory}
\label{theory}

The basic mechanism for the avoidable crossing of energy bands is to let the bands have different symmetry representations~\cite{yang2014,young2015,fang2015}. We start with classifying stable DSMs of 2D SOC systems in the presence of $\mathcal{PT}$ and crystalline symmetry such as rotational and mirror symmetries. The generic Hamiltonian has the form $\mathcal{H}(\mathbf{k})=\sum_{i,j=0}^3d_{ij}(\mathbf{k})\sigma_i\otimes\tau_j$, where $\sigma_i$ and $\tau_i$ ($i=1,2,3$) are Pauli matrices acting on the spin and orbital, respectively. $\sigma_0$ and $\tau_0$ are $2\times2$ identity matrices. $d_{ij}(\mathbf{k})$ are real functions of $\mathbf{k}$. The invariance of the system will set many terms to be zero or nonindependent, which reduces to the effective Hamiltonian in Eq.~(\ref{minimal}).

\subsection{Case A. $\mathcal{T}$, $\mathcal{P}$ conserved}
The $\mathcal{T}$ operator is $\mathcal{T}=i\sigma_2\mathcal{K}$, where $\mathcal{K}$ is complex conjugation. $\mathcal{H}(-\mathbf{k})=\mathcal{T}\mathcal{H}(\mathbf{k})\mathcal{T}^{-1}$. $\mathcal{P}$ is independent of the spin-rotation, where its specific form is determined by $\mathcal{P}^\dag\mathcal{P}=1$, $[\mathcal{T},\mathcal{P}]=0$, and $(\mathcal{PT})^2=-1$. Thus $\mathcal{P}=\pm\tau_0$, $\pm\tau_3$, or $\pm\tau_1$. $\mathcal{H}(-\mathbf{k})=\mathcal{P}\mathcal{H}(\mathbf{k})\mathcal{P}^{-1}$.
If $\mathcal{P}=\pm\tau_1$, $\Gamma^{(1,2,3,4,5)}=(\sigma_3\otimes\tau_3,\tau_2,\sigma_1\otimes\tau_3,\sigma_2\otimes\tau_3,\tau_1)$,
and $d_a(\mathbf{k})=-d_a(-\mathbf{k})$ for $a=1,\ldots,4$, $d_5(\mathbf{k})=d_5(-\mathbf{k})$. All $d_a(\mathbf{k})=0$ satisfied simultaneously leads to DPs, however, the number of equations is larger than the number of variables $k_x$, $k_y$ and the external parameter $m$. Therefore, the band crossing will not happen at generic $\mathbf{k}$. The crystalline symmetry will further set constraints on $d_a(\mathbf{k})$. In 2D, the relevant symmetries are $C_{2\hat{x}}$, $C_{2\hat{y}}$, $M_{\hat{z}}$, and $C_{n\hat{z}}$.

\begin{table}[b]
\caption{The classification table of 2D DSMs for \emph{case A} with both $\mathcal{T}$ and $\mathcal{P}$ symmetries. $\ell=2,4,6$. The DPs are induced by the topological band crossing at the BZ boundary.}
\begin{center}\label{table1}
\renewcommand{\arraystretch}{1.5}
\begin{tabular*}{3.4in}
{@{\extracolsep{\fill}}cccc}
\hline
\hline
Symmetry & $\mathcal{P}$ & $\mathcal{T}$ &  Dispersion
\\
\hline
$C_{2\hat{x}}=i\sigma_3\otimes\tau_3$ & $\pm\tau_1$ & $i\sigma_2\mathcal{K}$ &  Linear
\\
$M_{\hat{z}}=\sigma_3\otimes\tau_2$ & $\pm\tau_1$ & $i\sigma_2\mathcal{K}$ &  Linear
\\
$C_{\ell\hat{z}}=e^{i\frac{\pi}{\ell}\sigma_3}\otimes\tau_3$ & $\pm\tau_1$ & $i\sigma_2\mathcal{K}$ &  Linear
\\
$C_{6\hat{z}}=e^{i\frac{\pi}{2}\sigma_3}\otimes\tau_3$ & $\pm\tau_1$ & $i\sigma_2\mathcal{K}$ &  Cubic
\\
\hline
\hline
\end{tabular*}
\end{center}
\end{table}

\emph{$C_{2\hat{x}}$ symmetry.} The invariance of the system under $\pi$ rotation along the $x$ axis requires $C_{2\hat{x}}\mathcal{H}(k_x,k_y)C_{2\hat{x}}^{-1}=\mathcal{H}(k_x,-k_y)$. $k_y=0$ along the $k_x$ axis, $[C_{2\hat{x}},\mathcal{H}(k_x,0)]=0$. Therefore we can choose a basis to make both $C_{2\hat{x}}$ and $\mathcal{H}(k_x)$ diagonal. In such a basis, the explicit forms of $C_{2\hat{x}}$ is obtained by the constraint $[\mathcal{T},C_{2\hat{x}}]=0$. Therefore $C_{2\hat{x}}=\text{diag}[\alpha_p,\alpha_q,\alpha_p^*,\alpha_q^*]$, where $\alpha_p=\exp[i\pi(p+\frac{1}{2})]$ with $p=0,1$. Thus we have $C_{2\hat{x}}=i\sigma_3$, or $i\sigma_3\otimes\tau_3$.
With the explicit representations of $C_{2\hat{x}}$ and $\mathcal{P}$ operators, we can get the symmetry allowed forms of $d_a(\mathbf{k})$. Take $\mathcal{P}=\pm\tau_1$ for example, when $C_{2\hat{x}}=i\sigma_3\otimes\tau_3$, to the leading order in $\mathbf{k}$, the symmetry constraints leads to
$d_{1,2,3,4,5}(\mathbf{k})=(v_1k_x,v_2k_y,v_3k_y,v_4k_y,v_5k_xk_y)$,
where $v_i$ are coefficients. This is nothing but the 2D DSM with anisotropic linear dispersions. Considering the full periodic structure of the BZ, the time-reversal-invariant momentum (TRIM) $(\pi,0)$, $(0,\pi)$ and $(\pi,\pi)$ are possible locations of the 2D DPs. The $(0,0)$ point is simply excluded because the symmetry-protected DPs are not compatible with threefold rotation, while the space groups admit the 4D representations at $(0,0)$ contain $C_{3\hat{z}}$~\cite{young2012,bradley1972}. Similarly, when $C_{2\hat{x}}=i\sigma_3$,
$d_{1,2,3,4,5}(\mathbf{k})=(v_1k_x,v_2k_x,v_3k_y,v_4k_y,m+u_1k_x^2+u_2k_y^2)$, where $v_i$, $u_i$ and $m$ are constants. This is not DSM Hamiltonian. However, it is worth to mention that critical DPs indeed exist at TRIM at the transition between QSH and NI~\cite{murakami2007} as shown in Fig.~\ref{fig4}(b).
At TRIM point, $\mathbf{k}$ and $-\mathbf{k}$ are equivalent, and all odd functions in $\mathcal{H}(\mathbf{k})$ vanishes. By tuning $m$, $d_5(\mathbf{k},m)=0$ is satisfied.

The classification of 2D DSMs for \emph{case A} is shown in Table~\ref{table1}. The details on 2D DPs protected by $M_{\hat{z}}$ and $C_{n\hat{z}}$ are discussed in Appendix~\ref{appendix_A}. The stable DPs exist at TRIM of the BZ boundary and are protected by the lattice symmetry, including symmorphic and nonsymmorphic symmetries~\cite{young2015}. In nonsymmorphic symmetric systems with two sublattices,
$\mathcal{P}=\pm\tau_1$ is the inversion operator with sublattices interchanged~\cite{yang2014}, which is
realized by the inversion with respect to a lattice site followed by a fraction translation.

\begin{figure}[t]
\begin{center}
\includegraphics[width=3.3in]{fig1}
\end{center}
\caption{(color online). A nonsymmorphic symmetry $\{g|\mathbf{t}\}$ ensures band crossing on a $g$ invariant line at the BZ boundary. The doubly-degenerate bands are artificially split for clarity. (a) With $\mathcal{P}$ and $\mathcal{T}$, and $\mathcal{P}=\tau_1$, Bloch states  $(|u^+_{\mathbf{k}}\rangle,\mathcal{T}|u^-_{\mathbf{k}}\rangle)$ and $(|u^-_{\mathbf{k}}\rangle,\mathcal{T}|u^+_{\mathbf{k}}\rangle)$ carry different representations of $\{g|\mathbf{t}\}$ ($C_{2\hat{x}}$ for example) and cross at TRIM. (b) Without $\mathcal{P}$ and $\mathcal{T}$, but with combined $\mathcal{PT}$, $(|u^+_{\mathbf{k}}\rangle,\mathcal{PT}|u^+_{\mathbf{k}}\rangle)$ and $(|u^-_{\mathbf{k}}\rangle,\mathcal{PT}|u^-_{\mathbf{k}}\rangle)$ have different eigenvalues of $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}$ and cross avoidably.}
\label{fig1}
\end{figure}

\subsection{Case B. $\mathcal{T}$, $\mathcal{P}$ broken}

The anti-unitary $\mathcal{PT}$ symmetry reverses spins and keeps $\mathbf{k}$ invariant. Without loss of generality, in this section we set $\mathcal{PT}=i\sigma_2\mathcal{K}$. With $[\mathcal{PT},\mathcal{H}(\mathbf{k})]=0$, one has $\mathcal{H}(\mathbf{k})=d_1\tau_1+d_2\sigma_3\otimes\tau_2+d_3\sigma_1\otimes\tau_2+d_4\sigma_2\otimes\tau_2+d_5\tau_3$.
The crystalline symmetries fall into two classes: symmorphic and nonsymmorphic space group symmetries, the operations of which is denoted as $\{g|\mathbf{t}\}$ and can be constructed by the point group operations $g$ with translation $\mathbf{t}$ that are a full and a fraction of a Bravais lattice vector, respectively. In 2D, the representative symmorphic group symmetries are $\{C_{2\hat{x}}|00\}$, $\{C_{2\hat{y}}|00\}$,  $\{M_{\hat{z}}|00\}$ and $C_{n\hat{z}}$. The typical 2D nonsymmorphic group operations include screw rotation $\{C_{2\hat{x}}|\mathbf{t}\}$, $\{C_{2\hat{y}}|\mathbf{t}\}$, glide mirror lines $\{M_{\hat{x}}|\mathbf{t}\}$, $\{M_{\hat{y}}|\mathbf{t}\}$, and glide mirror plane $\{M_{\hat{z}}|\mathbf{t}\}$, where $\mathbf{t}$ is a half-translation which satisfying $g\mathbf{t}=\mathbf{t}$ and $e^{i\mathbf{G}\cdot\mathbf{t}}=-1$ for odd reciprocal lattice vector $\mathbf{G}$. Explicitly, in units of Bravais lattice constant, $\mathbf{t}=(\frac{1}{2},0),(0,\frac{1}{2}),(\frac{1}{2},\frac{1}{2})$. In the following, we will show that stable 2D DPs can only exist in the presence of nonsymmorphic symmetries.

(i). For the \emph{symmorphic group}, we take $\{M_{\hat{z}}|00\}$ for example, $\mathbf{k}$ is invariant and $(\sigma_1,\sigma_2,\sigma_3)\rightarrow(-\sigma_1,-\sigma_2,\sigma_3)$. The symmetry constraints on $d_a(\mathbf{k})$ is $d_{3,4}(\mathbf{k})=0$. To the lowest order in $\mathbf{k}$, the general form is $d_{a}(\mathbf{k})=m_a+v_a^xk_x+v_a^yk_y$ for $a=1,2,5$, and $m_a,v_a^x,v_a^y$ are coefficients. The DPs are not guaranteed and can only be generated by fune tuning. Even with the extra $\{C_{2\hat{x}}|00\}$ or $\{C_{2\hat{y}}|00\}$ symmetry, it only constrains $d_2(\mathbf{k})=v_2^yk_y$, still the DPs are not guaranteed.

\begin{table}[b]
\caption{The classification table of 2D DSMs for \emph{case B} with $\mathcal{T}$ and $\mathcal{P}$ broken. The nonsymmorphic symmetries with corresponding half-translation $\mathbf{t}$ are indicated.}
\begin{center}\label{table2}
\renewcommand{\arraystretch}{1.5}
\begin{tabular*}{3.4in}
{@{\extracolsep{\fill}}ccc}
\hline
\hline
Symmetry & Half-translation $\mathbf{t}$ & Dispersion
\\
\hline
$\{C_{2\hat{x}}|\mathbf{t}\}$ & $\mathbf{t}=(0,\frac{1}{2})$, $(\frac{1}{2},\frac{1}{2})$ & Linear
\\
$\{C_{2\hat{y}}|\mathbf{t}\}$ & $\mathbf{t}=(\frac{1}{2},0)$, $(\frac{1}{2},\frac{1}{2})$ & Linear
\\
$\{M_{\hat{x}}|\mathbf{t}\}$ & $\mathbf{t}=(\frac{1}{2},0)$, $(\frac{1}{2},\frac{1}{2})$ & Linear
\\
$\{M_{\hat{y}}|\mathbf{t}\}$ & $\mathbf{t}=(0,\frac{1}{2})$, $(\frac{1}{2},\frac{1}{2})$ & Linear
\\
\hline
\hline
\end{tabular*}
\end{center}
\end{table}

(ii). For the \emph{nonsymmorphic group}, we take $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}$ for example.
Under $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}$, $(x,y)\rightarrow(x+\frac{1}{2},-y+\frac{1}{2})$, $(k_x,k_y)\rightarrow(k_x,-k_y)$ and $(\sigma_1,\sigma_2,\sigma_3)\rightarrow(\sigma_1,-\sigma_2,-\sigma_3)$.
$\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}\mathcal{H}(k_x,k_y)\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}^{-1}=\mathcal{H}(k_x,-k_y)$.
Thus $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}^2=-e^{-ik_x}$ is just a full-translation along $x$ axis with a $2\pi$ rotation of spins.
$\mathcal{H}(\mathbf{k})$ is invariant under $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}$ along $k_y=0,\pi$ lines. Therefore both $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}$ and $\mathcal{H}(k_x)$ can be chosen as diagonal.
The exact representation of $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}$ operator is obtained by noticing
\begin{equation}\label{anticommutation} \{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}\mathcal{PT}=e^{-ik_x}e^{-ik_y}\mathcal{PT}\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}.
\end{equation}
Thus $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}=ie^{-ik_x/2}\sigma_3$, or $ie^{-ik_x/2}\sigma_3\otimes\tau_3$ when $k_y=0$; and $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}=ie^{-ik_x/2}\tau_3$ when $k_y=\pi$. However, there is only one case that stable band crossing and Dirac Hamiltonian can appear. Namely, $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}=ie^{-ik_x/2}\tau_3$ along $k_y=\pi$ line, only $d_5(k_x,\pi)$ term is present. If $d_5(k_x,\pi)=0$ at $\mathbf{k}_0=(k_{x0},\pi)$, around $\mathbf{k}_0$ point, to the leading order, the effective Hamiltonian is Dirac like and expanded as $d_{1,2,3,4,5}(\mathbf{k})=(v_1k_y,v_2k_y,v_3k_y,v_4k_y,v_5k_x)$, where $v_i$ are real. There is an alternative way to see why topological band crossing can only exist along $k_y=\pi$ line. The two Bloch states $|u^+_{\mathbf{k}}\rangle$, $|u^-_{\mathbf{k}}\rangle$ and their $\mathcal{PT}$ partners $\mathcal{PT}|u^+_{\mathbf{k}}\rangle$, $\mathcal{PT}|u^-_{\mathbf{k}}\rangle$ are eigenstates of $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}$, which is $\pm ie^{-ik_x/2}$. $|u^\pm_{\mathbf{k}}\rangle$ has the same energy as its $\mathcal{PT}$ partner. To have a stable band crossing, ($|u^+_{\mathbf{k}}\rangle,\mathcal{PT}|u^+_{\mathbf{k}}\rangle$) and ($|u^-_{\mathbf{k}}\rangle,\mathcal{PT}|u^-_{\mathbf{k}}\rangle$) should carry different representations of $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}$ as shown in Fig.~\ref{fig1}(b). This is only possible as $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}=ie^{-ik_x/2}\tau_3$ along $k_y=\pi$.

As shown in Table~\ref{table2}, the stable DPs can indeed exist in certain AFM semimetals respecting $\mathcal{PT}$ symmetry. The DPs in both $\mathcal{T}$-invariant and $\mathcal{T}$-broken semimetals are generated by the avoidable crossing of energy bands with distinct representations, which are protected by nonsymmorphic symmetries. However, different from the $\mathcal{T}$-invariant semimetals where the DPs reside at TRIMs of the BZ boundary, here the DPs in $\mathcal{T}$-broken systems can only appear at the BZ boundary but not at TRIMs in general. Moreover, the glide mirror plane symmetry $M_{\hat{z}}$ can not protect the DPs with $\mathcal{T}$-broken as listed in Appendix~\ref{appendix_B}, contrary to its protection of DPs with $\mathcal{T}$. More importantly, the DPs in the $\mathcal{T}$-broken AFM DSMs are locally permitted by crystalline symmetries, while the DPs in $\mathcal{T}$-invariant DSMs are essential and cannot be gapped without lowering the specific space group symmetries.

\section{Model}
\label{model}

\subsection{Tight-binding model}

Now we study a tight-binding model for a 2D AFM DSM to illustrate the nonsymmorphic symmetries protected DPs listed in Table~\ref{table2}. For direct comparison with $\mathcal{T}$ invariant case, we adopted the lattice similar to Ref.~\onlinecite{young2015} as shown in Fig.~\ref{fig2}(a), the system has a layered AFM structure with two atoms in one unit cell, denoted as $A$ and $B$. The lattice vectors are denoted as $\vec{\mathbf{a}}_1=(1,0,0)$, $\vec{\mathbf{a}}_2=(0,1,0)$. In each unit cell, the $A$ and $B$ atoms are shifted along the $z$ axis with the position $r_A=(-\frac{1}{4},-\frac{1}{4},-\frac{c}{2})$ and $r_B=(\frac{1}{4},\frac{1}{4},\frac{c}{2})$. Without AFM ordering, this structure is in space group No.~129 ($P4/nmm$). Each lattice site contains a $d_{z^2}$ orbital, where the square pyramidal crystal field splits $d_{z^2}$ from other $d$ orbitals. The Hamiltonian is
\begin{eqnarray}\label{H}
\mathcal{H} &=& t\sum\limits_{\langle
ij\rangle}c_{i}^{\dag}c_{j}+\sum\limits_{\langle\langle
ij\rangle\rangle}c_{i}^{\dag}\left[t_2+
i\lambda_{\text{so}}(\hat{\mathbf{d}}_1\times\hat{\mathbf{d}}_2)\cdot\mathbf{s}\right]c_{j}
\nonumber
\\
&&+\Delta\sum\limits_i\xi_ic_i^{\dag}\mathbf{s}\cdot\hat{\mathbf{n}}c_i,
\end{eqnarray}
where $\langle ij\rangle$ and $\langle\langle ij\rangle\rangle$
denote the nearest and next-nearest neighbor sites, respectively. $\lambda_{\text{so}}$ is SOC which involves spin-dependent next-nearest neighbor hopping, where $\hat{\mathbf{d}}_1$ and $\hat{\mathbf{d}}_2$ are unit vectors along two nearest neighbor bonds hopping from site $j$ to $i$~\cite{kane2005b}. $\mathbf{s}$ describes the electron spins. The third term is the staggered Zeeman term ($\xi_i=\pm1$) which describes the AFM ordering (along $\hat{\mathbf{n}}$ direction).

\begin{figure}[t]
\begin{center}
\includegraphics[width=3.3in]{fig2}
\end{center}
\caption{(color online). Energy band with DPs protected by nonsymmorphic symmetries in an AFM lattice respecting $\mathcal{PT}$. The DPs are marked as green dots in the BZ. (a) DPs are gapped and cannot be protected by $\{C_{2\hat{x}}|\frac{1}{2}0\}$ and $\{C_{2\hat{y}}|0\frac{1}{2}\}$ when $\hat{\mathbf{n}}=\hat{z}$. (b) Two DPs located along $X_1$-$M$ line are protected by $\{M_{\hat{x}}|\frac{1}{2}0\}$ solely when $\hat{\mathbf{n}}=\hat{x}$. (c) Distortion in the $\langle11\rangle$ direction eliminates $\{M_{\hat{x}}|\frac{1}{2}0\}$ but keeps $\{M_{\hat{z}}|\frac{1}{2}\frac{1}{2}\}$, gapping the DPs. (d) Alternatively, $B$ sites shifted in $y$ direction breaks $\{M_{\hat{z}}|\frac{1}{2}\frac{1}{2}\}$, but DPs still remain protected.}
\label{fig2}
\end{figure}

The symmetry of the system depends on the AFM order direction. In Fig.~\ref{fig2}(a), if $\hat{\mathbf{n}}=\hat{z}$, the system breaks $\mathcal{T}$ and $\mathcal{P}$ but respects $\mathcal{PT}$. However, from the energy bands, the screw axes $\{C_{2\hat{x}}|\frac{1}{2}0\}$ and $\{C_{2\hat{y}}|0\frac{1}{2}\}$ can not protect the DPs, consistent with the analysis in Table~\ref{table2}. This is quite different from the $\mathcal{T}$-, $\mathcal{P}$-invariant case (dashed line), where the DPs at $X_1$, $M$ and $X_2$ are protected by $\{C_{2\hat{x}}|\frac{1}{2}0\}$ and $\{C_{2\hat{y}}|0\frac{1}{2}\}$~\cite{young2015}. In Fig.~\ref{fig2}(b), if $\hat{\mathbf{n}}=\hat{x}$, two DPs are located along the $X_1$-$M$ line, and they are at different energies in the presence of $t_2$ term. The Hamiltonian
\begin{eqnarray}
\mathcal{H}(\mathbf{k}) &=& 4t\tau_1\cos\frac{k_x}{2}\cos\frac{k_y}{2}+2t_2(\cos k_x+\cos k_y)+
\nonumber
\\
&&2\lambda_{\text{so}}\sigma_2\otimes\tau_3\sin k_x+(\Delta-2\lambda_{\text{so}}\sin k_y)\sigma_1\otimes\tau_3.
\nonumber
\end{eqnarray}
The DPs are at $\mathbf{k}_1$=$(\pi,k_{y0})$ and $\mathbf{k}_2$=$(\pi,\pi-k_{y0})$, where $\sin k_{y0}=\Delta/2\lambda_{\text{so}}$. These two inequivalent DPs are protected by $\{M_{\hat{x}}|\frac{1}{2}0\}$. This can be seen by examining the effective model near these points. Near $\mathbf{k}=\mathbf{k}_1$,
\begin{eqnarray}
\mathcal{H}(\mathbf{k}_1+\mathbf{q})
&=& -2t\cos\frac{k_{y0}}{2}\tau_1q_x-2\lambda_{\text{so}}\sigma_2\otimes\tau_3q_x
\nonumber
\\
&&-2\lambda_{\text{so}}\cos k_{y0}\sigma_1\otimes\tau_3q_y.
\end{eqnarray}
At $\mathbf{k}_1$, the symmetry $\mathcal{PT}=i\sigma_2\mathcal{K}\otimes\tau_1$ allows the mass terms $\tau_2$ and $\sigma_3\otimes\tau_3$. This is forbidden by $\{M_{\hat{x}}|\frac{1}{2}0\}=i\sigma_1\otimes\tau_3$, but is allowed by $\{C_{2\hat{y}}|0\frac{1}{2}\}=\sigma_2\otimes\tau_2$. There is a subtle point, it seems $\{M_{\hat{z}}|\frac{1}{2}\frac{1}{2}\}=\sigma_3\otimes\tau_1$ could also forbid the mass terms.
To further clarify the role of $\{M_{\hat{x}}|\frac{1}{2}0\}$ and $\{M_{\hat{z}}|\frac{1}{2}\frac{1}{2}\}$ in protecting the DPs, we distorted the lattice which eliminates $\{M_{\hat{x}}|\frac{1}{2}0\}$ and keeps $\{M_{\hat{z}}|\frac{1}{2}\frac{1}{2}\}$ in Fig.~\ref{fig2}(c). Such distortion adds a term to nearest-neighbor hopping $\mathcal{H}_1=t_3\sin(k_x/2)\sin(k_y/2)\tau_1$, which will gap the DPs at $\mathbf{k}_{1,2}$, consistent with that DPs are not protected by $\{M_{\hat{z}}|\frac{1}{2}\frac{1}{2}\}$. In Fig.~\ref{fig2}(d), $B$ site is shifted along $y$ axis, which breaks $\{M_{\hat{z}}|\frac{1}{2}\frac{1}{2}\}$ but keeps $\{M_{\hat{x}}|\frac{1}{2}0\}$. Such reduced symmetry allows a term $\mathcal{H}_2=\cos(k_x/2)\sin(k_y/2)\left(t_4\tau_1+t_5\tau_2\right)$, where the DPs located along the $X_1$-$M$ line remain protected. However, it is noted that these two inequivalent DPs always appear or disappear together.

\subsection{Symmetry breaking}

\begin{figure}[b]
\begin{center}
\includegraphics[width=3.3in]{fig3}
\end{center}
\caption{(color online). Symmetry breaking phases. (a) Breaking $\mathcal{PT}$ by SIA, while preserving $\{M_{\hat{x}}|\frac{1}{2}0\}$, leads to Weyl points (marked as yellow and green dots) on the $k_x=\pi$ line. (b) Breaking both $\mathcal{PT}$ and $\{M_{\hat{x}}|\frac{1}{2}0\}$ results in the QAH insulator. In the edge spectrum, the green and red lines are edge states located at upper and lower edges, respectively.}
\label{fig3}
\end{figure}

Lowering the symmetry by external perturbations provides a toolbox to explore a wealth of topological phases. We start with the system as shown in Fig.~\ref{fig2}(d). The $\{M_{\hat{x}}|\frac{1}{2}0\}$ breaking by distortion in Fig.~\ref{fig2}(c) leads to a NI with gapped DPs. Displacing the $B$ sites along $[11]$ or $[1\bar{1}]$ always lead to NI phases. Quite different from $\mathcal{T}$-, $\mathcal{P}$-invariant 2D DSM where it lies at the boundary between QSH and NI~\cite{young2015}, the 2D AFM DSM lies in the NI phases deeply. This is simply because the edge of a physical system always breaks $\mathcal{PT}$ and leads to gapped edge states, even though a nontrivial bulk $Z_2$ index may be defined. Similarly, the structure inversion asymmetry (SIA) by applying an electric field along $z$ axis naturally breaks $\mathcal{PT}$ and adds a term $V\tau_3$, results in nondegenerate bands where DPs split into Weyl points as shown in Fig.~\ref{fig3}(a). Such Weyl points are located along the $k_x=\pi$ line and are not protected by symmetry. More interestingly, if both $\mathcal{PT}$ and $\{M_{\hat{x}}|\frac{1}{2}0\}$ are broken, a quantum anomalous Hall (QAH) state can be realized~\cite{dong2016} and electrically controllable as shown in Fig.~\ref{fig3}(b). Furthermore, with the proximity effect to an $s$-wave superconductor, the realization of a chiral topological superconductor with external tunability is expected~\cite{qi2010b,wang2015c,he2016,wang2016b}.

\section{Discussion and Conclusion}
\label{discussion}

The AFM long range order could exist in 2D. It should be interesting to study how the magnetic fluctuation and Coulomb interaction affect the stability of DPs, which is left to future work. In terms of realistic materials,
the actual existence of such phases in known materials remains an open question. However, similar to TI materials, the 2D magnetic DSMs may exist in materials with both strong SOC and specific space group symmetries. Furthermore, we comment on the possible candidate in the \emph{A}MnBi$_2$ family of compounds (\emph{A} = Ca, Sr, Eu, Yb)~\cite{kim2011,lee2013,guo2014,borisenko2015,petrovic2016,masuda2016}. The transport experiments have already shown the 2D Dirac fermion behaviors in the bulk materials~\cite{masuda2016,petrovic2016}. Take SrMnBi$_2$ for example, its bulk structure is in space group No.~139 ($I4/mmm$). For a monolayer of this material, the AFM order on Mn atoms is along the [001] direction which breaks both $\mathcal{T}$ and $\mathcal{P}$ whereas $\mathcal{PT}$ still holds. The Bi atoms form two layers of square lattice which is similar to the structure in Fig.~\ref{fig2}(a), which determine most of the electronic structure~\cite{kim2011,lee2013,guo2014}. The nonsymmorphic symmetries are $\{C_{2\hat{x}}|\frac{1}{2}\frac{1}{2}\}$ and $\{C_{2\hat{y}}|\frac{1}{2}\frac{1}{2}\}$. Therefore, the 2D DPs if exist in this material, can be protected by these symmetries in principle. Unfortunately, this system shows \emph{massive} Dirac fermions around the Fermi level, with Dirac mass gap along the line $\Gamma$-$M$ and $M$-$X$~\cite{kim2011,lee2013}. In fact, it is similar to the model in Eq.~(\ref{H}) with AFM ordering $\hat{\mathbf{n}}$ along $[001]$ direction.

In summary, we show that the 2D Dirac fermions could exist in AFM semimetals with $\mathcal{PT}$ symmetry, where the nonsymmorphic space group symmetries play an essential role. The realistic AFM DSM materials remain unknown. However, considering the ongoing rapid progress in the field of 2D materials~\cite{novoselov2016}, together with the nonsymmorphic symmtries listed in Table~\ref{table1} as a guidance, we are thus optimistic about the material search of the 2D AFM DSMs.

{\it Note added}: Recently, we became aware of an independent work on a similar problem~\cite{young2016}. However, their approach to 2D magnetic DSMs is different from our results.

\begin{acknowledgments}
The author is grateful to Haijun Zhang, Biao Lian and Shoucheng Zhang for insightful discussions. This work is supported by the National Thousand-Young-Talents Program; the National Key Research Program of China under Grant No.~2016YFA0300703; the Natural Science Foundation of Shanghai under Grant No.~17ZR1442500; the Open Research Fund Program of the State Key Laboratory of Low-Dimensional Quantum Physics, through Contract No.~KF201606; and by Fudan University Initiative Scientific Research Program.
\end{acknowledgments}

\begin{appendix}

\section{$C_{n\hat{z}}$ and $M_{\hat{z}}$ symmetries in \emph{case A}}
\label{appendix_A}

The representations of the $\Gamma$ matrices depends on the crystalline symmetry. If $\mathcal{P}=\pm\tau_0$,
we get $\Gamma^{(1,2,3,4,5)}=(\tau_1,\sigma_3\otimes\tau_2,\sigma_1\otimes\tau_2,\sigma_2\otimes\tau_2,\tau_3)$
and $d_a(\mathbf{k})=d_a(-\mathbf{k})$ for $a=1,\ldots,5$. Similarly, if $\mathcal{P}=\pm\tau_3$,
$\Gamma^{(1,2,3,4,5)}=(\sigma_3\otimes\tau_1,\tau_2,\sigma_1\otimes\tau_1,\sigma_2\otimes\tau_1,\tau_3)$,
and $d_a(\mathbf{k})=-d_a(-\mathbf{k})$ for $a=1,\ldots,4$, $d_5(\mathbf{k})=d_5(-\mathbf{k})$.

\begin{figure}[b]
\begin{center}
\includegraphics[width=3.3in]{fig4}
\end{center}
\caption{(color online). (a) Nonsymmorphic symmetries protected 2D DSMs with DPs located at BZ boundary. (b) 2D DPs with accidental degeneracy appear at the phase transition between QSH and NI controlled by parameter $m$.}
\label{fig4}
\end{figure}

Here we give a short summary for the representations of $M_{\hat{z}}$ and $C_{n\hat{z}}$, and then give a detailed analysis for $C_{n\hat{z}}$ symmetry.

(i). \emph{$M_{\hat{z}}$ symmetry.} $M_{\hat{z}}$ is the mirror symmetry with respect to the $xy$ plane and requires $[M_{\hat{z}},\mathcal{H}(\mathbf{k})]=0$. It can be viewed as $M_{\hat{z}}=C_{2\hat{z}}\mathcal{P}$, where $C_{2\hat{z}}$ is the $\pi$ rotation along the $z$ axis. $M_{\hat{z}}$ follows the constraints $[\mathcal{T},M_{\hat{z}}]=0$, $M_{\hat{z}}^\dag M_{\hat{z}}=1$ and $M^2_{\hat{z}}=e^{i\phi}$.
The explicit form of $M_{\hat{z}}$ is (i) $M_{\hat{z}}=i\sigma_3\otimes\tau_1$, or $\sigma_3\otimes\tau_2$ when $|\mathcal{P}|=\tau_1$; (ii) $M_{\hat{z}}=i\sigma_3$, or $i\sigma_3\otimes\tau_3$ when $|\mathcal{P}|=\tau_0,\tau_3$. The stable DPs can appear with $M_{\hat{z}}=\sigma_3\otimes\tau_2$ and $\mathcal{P}=\pm\tau_1$.

(ii). \emph{$C_{n\hat{z}}$ symmetry.} $C_{n\hat{z}}$ is $n$-fold rotation along the $z$ axis, with $n=2,3,4,6$. The stable DPs can appear with $C_{\ell\hat{z}}$ ($\ell=2,4,6$) and $\mathcal{P}=\pm\tau_1$.

$C_{n\hat{z}}$ is the $n$-fold rotation along the $z$ axis, with $n=2,3,4,6$. In order to obtain the explicit form of $C_{n\hat{z}}$, we can choose a basis where $C_{n\hat{z}}$ is diagonal as $C_{n\hat{z}}=\text{diag}[u_{A}^{\uparrow},u_{B}^{\uparrow},u_{A}^{\downarrow},u_{B}^{\downarrow}]=\text{diag}[\alpha_p,\alpha_q,\alpha_p^*,\alpha_q^*]$, where $\alpha_p=\exp[i\frac{2\pi}{n}(p+\frac{1}{2})]$ with $p=0,1,\ldots,n-1$. The invariance of the system under $C_{n\hat{z}}$ requires
\begin{equation}\label{constraint}
C_{n\hat{z}}\mathcal{H}(k_+,k_-)C_{n\hat{z}}^{-1}=\mathcal{H}(k_+e^{i\frac{2\pi}{n}},k_-e^{-i\frac{2\pi}{n}}),
\end{equation}
where $k_\pm=k_x\pm ik_y$. Therefore, at the TRIM, the Hamiltonian commutes with $C_{n\hat{z}}$, i.e.,
\begin{equation}
[C_{n\hat{z}},\mathcal{H}(\mathbf{k}_{\text{TRIM}})]=0.
\end{equation}
This is different from $C_{2\hat{x}}$ and $M_{\hat{z}}$ cases, where invariant lines or plane exist in the BZ. We find that stable DPs can appear and are protected by $C_{2\hat{z}}$, $C_{4\hat{z}}$, $C_{6\hat{z}}$ together with $\mathcal{P}=\pm\tau_1$. The classification of 2D DSM protected by $C_{n\hat{z}}$ symmetry in \emph{case A} is listed in Table.~\ref{table1}. From Table~\ref{table1}, we can see that the $C_{4\hat{z}}$ case is consistent with the results obtained in Ref.~\onlinecite{young2015}. The $C_{2\hat{z}}$ case is not independent from $M_{\hat{z}}$ case in the main text, for in the presence of $\mathcal{P}$, $C_{2\hat{z}}=M_{\hat{z}}\mathcal{P}$.

Below we show the basic steps to see when the stable DPs can appear. As we discussed in the main text, the basic mechanism for the avoidable crossing of energy bands is to let the bands have different symmetry representations. The four Bloch states can be chosen as $|u^\uparrow_A\rangle$, $|u^\uparrow_B\rangle$, $|u^\downarrow_A\rangle$, and $|u^\downarrow_B\rangle$. Here $|u^\uparrow_A\rangle$ and $|u^\uparrow_B\rangle$ have the $C_{n\hat{z}}$ eigenvalues $u_{A}^{\uparrow}$ and $u_{B}^{\uparrow}$, respectively. There are two distinct cases,

(i). When $\mathcal{P}=\pm\tau_0,\pm\tau_3$, the inversion operation will not flip the orbitals, thus the energy $E_{A,\uparrow}=E_{A,\downarrow}$, and $E_{B,\uparrow}=E_{B,\downarrow}$ at $\mathbf{k}_{\text{TRIM}}$. Therefore, the four Bloch states are grouped into $(u_{A}^{\uparrow},u_{A}^{\downarrow})$ and $(u_{B}^{\uparrow},u_{B}^{\downarrow})$. In order to have a stable band crossing, one \emph{necessary} condition is that these two groups must have different eigenvalues of $C_{n\hat{z}}$. Namely, $(u_{A}^{\uparrow},u_{A}^{\downarrow})$ are different from $(u_{B}^{\uparrow},u_{B}^{\downarrow})$. Due to $[C_{n\hat{z}},\mathcal{H}({\mathbf{k}_{\text{TRIM}}})]=0$, both $C_{n\hat{z}}$ and $\mathcal{H}({\mathbf{k}_{\text{TRIM}}})$ can be diagonal. The only possible form is $\mathcal{H}({\mathbf{k}_{\text{TRIM}}})=d_5\tau_3$. However, $d_5$ is even under $\mathcal{P}$, thus $d_5=m$, where $m$ is the external parameters. In general, $m\neq0$, and Eq.~(\ref{constraint}) can not constrain $m$ to be zero, therefore, the stable Dirac points are not possible in this case.

(ii). Similarly, when $\mathcal{P}=\pm\tau_1$, inversion operation will switch the orbitals, therefore $E_{A,\uparrow}=E_{B,\downarrow}$, and $E_{B,\uparrow}=E_{A,\downarrow}$. Thus, $(u_{A}^{\uparrow},u_{B}^{\downarrow})$ must be different from $(u_{B}^{\uparrow},u_{A}^{\downarrow})$. This necessary condition rules out $C_{3\hat{z}}$. Moreover, we consider the only possible form of $\mathcal{H}$ at $\mathbf{k}_{\text{TRIM}}$ is $\mathcal{H}({\mathbf{k}_{\text{TRIM}}})=d_1\sigma_3\otimes\tau_3$. However, $d_1(\mathbf{k})$ is an odd function, thus $d_1$ vanishes. Furthermore, using Eq.~(\ref{constraint}) one can get the possible forms of $d_a(\mathbf{k})$. Here we take $C_{4\hat{z}}$ for example, from the above symmetry analysis, the only possible form is $C_{4\hat{z}}=e^{\pm i\frac{\pi}{4}\sigma_3}\otimes\tau_3$. At the TRIM, all odd functions $d_{1,2,3,4}$ vanish, while $d_5(\mathbf{k}_{\text{TRIM}})$ is an even function. But $C_{4\hat{z}}$ also force $d_5(\mathbf{k}_{\text{TRIM}})$ to be odd. Therefore $d_5$ vanishes at the TRIM, and DPs can appear the such TRIMs. At the TRIM, the Hamiltonian is expanded as $d_{1,2,3,4,5}(\mathbf{k})=(0,0,v k_x,v k_y,0)$. More interestingly, $C_{6\hat{z}}$ symmetry could support the cubic Dirac fermions.

Here we point out, in \emph{case A} with $\mathcal{P}$ symmetry, one can always construct the mirror symmetry $M_{\hat{n}_{\perp}}$ ($\hat{n}_\perp\perp\hat{z}$) from $C_{2\hat{n}_{\perp}}$ and $\mathcal{P}$ as $M_{\hat{n}_{\perp}}=C_{2\hat{n}_{\perp}}\mathcal{P}$. Here $\hat{n}_{\perp}$ can chosen to be $\hat{x}$ or $\hat{y}$. Therefore, $M_{\hat{n}_{\perp}}$ is not independent from $C_{2\hat{n}_{\perp}}$, and is not listed in Table~\ref{table1}. However, in \emph{case B} without $\mathcal{P}$, $M_{\hat{n}_{\perp}}$ is different from $C_{2\hat{n}_{\perp}}$, as we can see from Table~\ref{table1}.

\section{Nonsymmorphic symmetry with no DPs protection in \emph{case B}}
\label{appendix_B}

Here we list in Table~\ref{table3} the nonsymmorphic symmetries which cannot protect the DPs in \emph{case B}.

\begin{table}[htbp]
\caption{The nonsymmorphic symmetries which cannot protect the DPs in \emph{case B} with $\mathcal{T}$ and $\mathcal{P}$ broken.}
\begin{center}\label{table3}
\renewcommand{\arraystretch}{1.5}
\begin{tabular*}{3.4in}
{@{\extracolsep{\fill}}cccc}
\hline
\hline
Symmetry & half-translation $\mathbf{t}$ & DP protection
\\
\hline
$\{C_{2\hat{x}}|\mathbf{t}\}$ & $\mathbf{t}=(\frac{1}{2},0)$ & No
\\
$\{C_{2\hat{y}}|\mathbf{t}\}$ & $\mathbf{t}=(0,\frac{1}{2})$ & No
\\
$\{M_{\hat{x}}|\mathbf{t}\}$ & $\mathbf{t}=(0,\frac{1}{2})$ & No
\\
$\{M_{\hat{y}}|\mathbf{t}\}$ & $\mathbf{t}=(\frac{1}{2},0)$ & No
\\
$\{M_{\hat{z}}|\mathbf{t}\}$ & $\mathbf{t}=(\frac{1}{2},0), (0,\frac{1}{2}), (\frac{1}{2},\frac{1}{2})$ & No
\\
\hline
\hline
\end{tabular*}
\end{center}
\end{table}

\section{Electrically controllable QAH state}
\label{appendix_C}

The parameters in Fig.~\ref{fig3}(b) of main text is $t_2=-0.1$, $\lambda_{\text{so}}=0.5$, $\Delta=0.3$, $V=0.9$, and the four nearest-neighbor hopping are $t_{AB}=(1.5, 1.3, 1.3, 0.7)$. Meanwhile, the QAH state is electrically tunable, i.e., small $V$ leads to NI phase as shown in Fig.~\ref{fig5}.

\begin{figure}[b]
\begin{center}
\includegraphics[width=2.8in]{fig5}
\end{center}
\caption{The electrically tunable QAH state in both $\mathcal{PT}$ and $\{M_{\hat{x}}|\frac{1}{2}0\}$ broken AFM DSMs. The parameters $V=0.2$ and others are the same as in Fig.~\ref{fig3}(b) of main text, which leads to a NI phase.}
\label{fig5}
\end{figure}

\end{appendix}