# Enhanced Multi-Parameter Cognitive Architecture for Future Wireless
Communications

Kaiqing Zhang, Feifei Gao, and Qihui Wu K. Zhang and F. Gao are with Tsinghua
National Laboratory for Information Science and Technology, Beijing, 10084,
China<EMAIL_ADDRESS>feifeigao@ieee.org)Q. Wu is with
the College of Communications Engineering, PLA University of Science and
Technology, Nanjing, 211101, China<EMAIL_ADDRESS>
###### Abstract

The very original concept of cognitive radio (CR) raised by Mitola targets at
all the environment parameters, including those in physical layer, MAC layer,
application layer as well as the information extracted from reasoning. Hence
the first CR is also referred to as “full cognitive radio”. However, due to
its difficult implementation, FCC and Simon Haykin separately proposed a much
more simplified definition, in which CR mainly detects one single parameter,
i.e., spectrum occupancy, and is also called as “spectrum sensing cognitive
radio”. With the rapid development of wireless communication, the
infrastructure of a wireless system becomes much more complicated while the
functionality at every node is desired to be as intelligent as possible, say
the self-organized capability in the approaching 5G cellular networks. It is
then interesting to re-look into Mitola’s definition and think whether one
could, besides obtaining the “on/off” status of the licensed user only,
achieve more parameters in a cognitive way. In this article, we propose a new
cognitive architecture targeting at multiple parameters in future cellular
networks, which is a one step further towards the “full cognition” compared to
the most existing CR research. The new architecture is elaborated in detailed
stages, and three representative examples are provided based on the recent
research progress to illustrate the feasibility as well as the validity of the
proposed architecture.

###### Index Terms:

Multiple parameters cognition, detection and estimation, machine learning,
future communication systems.

##  I Introduction

The rapid development of wireless communications has engendered quick
proliferation of media-rich mobile devices as well as significant enhancement
of communication ability. The actual mobile traffic in 2010 is five times
greater than an official forecast made by the International Telecommunication
Union (ITU) in 2005 and would continuously increase about 1000 times till the
end of 2020 [1].

The system capacity could be enhanced by enlarging the network coverage (via
relays, femtocell), increasing the space dimension (via massive MIMO),
improving the bandwidth efficiency and energy efficiency (via cognitive radio
technology, green communications) [2]. On the other side, dynamic networks are
adopted to speed up the service innovation in a more intelligent way. For
example, the self-organization networks (SON) with artificial intelligence
serves as a promising solution to address the challenges imposed by large-
scale networks, e.g., the high cost of configuring and managing networks, the
fluctuation nature of the available spectrum, and the diverse QoS requirements
of various applications[3, 4]. In order to realize these key technologies,
cross-layer parameters, e.g., spectrum occupancy, transmit power, modulation,
constellation, channel coding, location, and cell edge, should be achieved and
shared among nodes in the network. However, the current information exchanging
scheme still relies on the cooperative feedback among nodes, which causes
severe transmission overhead especially when the scale of the network and the
amount of the data traffic are large. Hence, an intelligent way to cognitively
obtain as many parameters as possible at every node could greatly enhance the
network efficiency.

The terminology “full cognitive radio” has already been proposed by Mitola [5]
in which the cognition targets at every possible parameter observable by a
wireless node or network. However, this ultimate goal of cognitive radio (CR)
only stays in the conceptual level, while the recent 10 years’ studies in
parameter achieving mainly focused on one single parameter, i.e., the
occupancy of a particular spectrum. Various parameters depicting the network
status, e.g., channel occupancy, transmit power level, signal constellation,
modulation scheme, channel coding, as well as cell coverage, network topology,
user preferences, communication protocols, and sensing policies, have not been
cognized under a unified framework. It is then necessary to expand the
parameter space from a single “spectrum hole” to “multiple parameters” space
in order to make the network more “intelligent”.

Meanwhile, the multi-parameter space also shares the highly structured and
mutually coupled characteristics, which could be further utilized to enhance
the cognition performance and reduce cognition complexity. If sufficient prior
knowledge are available at each node, the cognition of multiple parameters
will be similar to the conventional estimation and detection problem for
single parameter cognition but with an enlarged dimension of parameter space.
In this case, modified signal processing techniques, e.g., multiple hypotheses
testing, could be applied [6]. Nonetheless, in complicated networks, the prior
knowledge is normally unavailable and conventional signal processing
techniques can be hardly executed. A solution for multiple parameters
cognition is then to introduce the machine learning theory that is able to
establish the prior knowledge base automatically from the proliferation of
traffic data. The amount of data could be accumulated shortly considering the
signal transmission speed over Mega or Giga bits. In fact, machine learning
has already emerged as a booming research area these years, whose core idea is
to mine the “patterns” behind massive data and further identify or predict
unknown “patterns”. The inherent coupling and inter-relationship of network
parameters can be recognized as the “patterns” through learning techniques in
this context.

Practically, the parameters that determine the control and optimization
policies in networks change frequently, while a small change in environmental
parameters may result in a big change in the systematic behaviors [4]. The
parameter cognition should then be as adaptive and predictive as possible to
cope with the dynamic characteristics of parameters. If the cognition is able
to provide accurate prediction of the network status from the patterns behind
multi-parameter space, then better quality of service can be achieved via
resource pre-distribution.

In this article, we propose a new cognitive architecture based on either the
conventional signal processing or the machine learning techniques, targeting
at multiple parameters in future wireless communication networks. The
remainder of the article is organized as follows. We first describe the core
idea and the framework of the proposed cognitive architecture. We also provide
three concrete examples that demonstrate the validity and efficacy of the
architecture. Summary and some prospects are highlighted at the end of the
article.

##  II Architecture of Enhanced Multi-Parameter Cognitive Techniques

The architecture of the proposed cognitive technique is elaborated as three
main stages: parameter structuralization, multi-parameter cognition, and
parameter prediction.

###  II-A Parameter Structuralization

Figure 1: Enormous number of parameters involved in cognitive communication
networks.

As illustrated in Fig. 1, a communication network is described by gigantic
amount of parameters in different layers, and cognition for as many as
possible parameters at each node could greatly release the expense of the
feedback and enhance network efficiency.

Quite many efforts have been made in parameter cognition such as channel
estimation and spectrum detection, where some fine results have been
established in the past decades. However, these efforts each concerns the
extraction of only an individual parameter, ignoring the potential
relationship among different parameters, i.e., parameter structure. The
structure of parameters in the same layer can be achieved through direct
signal processing techniques while structure in different layers can be
achieved through pattern recognition and learning techniques. For example, the
cognized symbol constellation will directly exhibit the spectrum occupancy
status, and could also indicate the transmission behaviours or preferences of
users through certain data mining approach.

There are three main advantages of parameter structuralization: (i)
Structuralization could improve the cognition accuracy since various
parameters complement one another in the light of the “structure”; (ii)
Complexity of cognition could be reduced since smaller number of parameters
are necessary to be cognized due to the correlation among them; (iii) Range of
parameter cognition could be expanded because some non-cognitable parameters
could be achieved by reasoning from the cognitable parameters via machine
learning techniques.

###  II-B Multi-Parameter Cognition

Figure 2: Framework of the multi-parameter space cognitive architecture.

The proposed multi-parameter cognitive architecture is shown in Fig. 2, where
the most popular parameters include spectrum occupancy, transmit power,
modulation, and user behaviour.

In general, parameter cognition can be categorized into two cases: (i) prior-
sufficient case; (ii) prior-deficient case, depending on how much the prior
knowledge is available before cognition. For the prior-sufficient case, where
the key information, e.g., interference level, noise characteristics, channel
statistics are known to users, the multiple parameters cognition problem could
be solved by conventional signal processing techniques. In this case, closed-
form solutions or numerical solutions normally exist, which sometimes serve as
the performance bound for other low-complex sub-optimal algorithms.

When the prior knowledge is insufficient to execute conventional signal
processing algorithms, one may resort to learning techniques, e.g., data-
driven algorithms from machine learning and pattern recognition, to
intelligently establish the cognition knowledge base. Indeed, some machine
learning algorithms have been proposed for various tasks in CR parameter
cognition. For example, in [7], a spectrum sensing engine based on support
vector machine (SVM) was designed, advancing the sensing performance with
smaller samples compared to the energy detector. In cooperative spectrum
sensing paradigm, [8] presents reinforce learning methods to reduce the
sensing overhead by releasing network flow congestion. Most of these learning
techniques, however, are mainly designed to obtain a single parameter, i.e.,
spectrum occupancy, and do not consider the learning of other system
parameters.

Learning-based cognition techniques seem to be more applicable in real CR
setups and are subject to little performance downgrade in the above cases.
Nonetheless, there are still some deficiencies that are worth further
attention for there is no “free lunch”. First, it is challenging to collect
sufficient “clean” data for learning and recognition. Extra overhead and high
complexity for establishing and maintaining the database during the whole
cognition are unavoidable. Second, in some centralized control and
optimization scenarios, learning techniques may not be adequate to cope with
large-scale traffic data sets. Hence, it may be necessary to incorporate non-
parametric learning algorithms with low computational complexity in the
cognition framework. Moreover, data-driven algorithms can be prone to data
falsification and attack from malicious users, and more efforts are necessary
for safety analysis and protection operations.

###  II-C Parameter Prediction

Following the multiple parameters cognition, users are able to utilize
internal structure as well as underlying pattern of parameters to further
understand the parameter evolution. For example, channel occupancy status can
be predicted by learning the traffic characteristics of licensed systems using
the neural network model [9]. The reliable prediction of channel status
considerably reduces the overhead consumed in spectrum sensing for the reason
that only those channels predicted to be idle need to be sensed in the next
time-slot. Moreover, if the unlicensed user can predict other parameters of
licensed users, e.g., power level, modulation, and coding-scheme, it could
then adjust the corresponding parameters such that the interference caused to
the licensed user can be optimally reduced. In intelligent communication
networks, e.g., SON, predicability normally serves as the counterpart of
scalability but is difficult to model due to the complexity of systems [4]. In
this case, learning-based cognition techniques is much more applicable for
parameter prediction because of its extrapolative nature and simple
implementation.

##  III Examples

In the following, we provide three representative examples of the latest
research progress to show how the signal processing and machine learning
techniques could be used to cognize multiple parameters.

###  III-A Example 1: Cognition for Transmit Power Levels

Let us first consider a practical CR scenario, where the licensed user could
work under more than one discrete power levels, as opposed to single power
level in conventional CR. The parameter that could be cognized is then not
only the on/off status of the licensed user, but also its transmit power
levels, where non-zero power levels indicates the “on” of licensed user and
vice versa.

If the candidate power levels and noise statistics are known in prior, the
spectrum occupancy as well as the current transmit power levels can be
obtained via multiple hypotheses detection [6]. Besides the false alarm
probability and detection probability that are used to evaluate the
conventional spectrum sensing, one may further define a new performance metric
called discrimination probability, which describes the capability of correctly
recognizing each power pattern. A numerical result is shown in Fig. 3, where
the theoretical curves can serve as an upper bound of performance for any sub-
optimal detectors.

Figure 3: The theoretical performance metric of cognition versus number of
samples with system parameters known in prior.

However, the assumption that the unlicensed user fully knows all prior
information cannot be true in realistic setup. In this case, the machine
learning based cognition techniques, such as clustering analysis and
classifier construction, could be exploited to replace the multiple hypotheses
testing approach. The unlicensed users can collect energy statistics and
formulate energy feature vector through multi-slot sensing scheme [10], and
licensed transmission patterns can be discovered by clustering analysis from a
sufficient number of energy feature vectors. Specifically, feature vectors
that share the same transmit power of PU can be grouped together so that the
number of power states and the corresponding vector clusters are determined.
Hence, the transmit power levels can be evaluated by the average of the
vectors in each cluster, from which the knowledge about transmit patterns and
preferences, e.g., the power states of the licensed user, the average power
level of each state, etc., are learned. Moreover, each energy feature vector
would be labeled by the cluster number it is partitioned to. Classifiers and
decision boundaries between different power states can then be constructed
through supervised learning on the basis of these labeled data. A numerical
example where SVM is used to train the classifiers and perform sensing tasks
is shown in Fig. 4.

Figure 4: Clustering results and SVM trained decision boundary on normalized
energy feature vectors when average signal noise ratio (SNR) is -12 dB.

Once the power level is detected, the unlicensed users are able to adjust
their own transmit power such that the interference temperature for that
particular power level is fulfilled,111FCC has regulated different protection
for different powered services in their recent reports [11]. by which means
the unlicensed users could fully squeeze the tolerance of the licensed user
and maximize their own throughput.

###  III-B Example 2: Cognition for Modulation Pattern

In this example, we address the problem of modulation recognition (MR) in CR
scenario. For conventional MR, there are two main genres of methods: the
theoretical maximum-likelihood-based methods and the statistical pattern-
recognition-based methods [12, 13].

However, MR in CR differs from the convention in several ways. Firstly, the
existing MR methods assume a prior knowledge of the candidate signal
constellations, which can be denoted as the modulation dictionary. It is
obligatory for the dictionary to contain as many potential modulation types as
possible while a redundant dictionary would definitely degrade the recognition
performance, especially in low signal-to-noise ratio region. Hence the
modulation dictionary should be pruned tactfully before recognition. Secondly,
the existing methods always assume the transmitter to be “on”, which is not
the case in CR scenario where the channel occupancy status is also unknown and
its detection is coupled with MR. Thirdly, the existing methods do not
consider different transmit power levels, which can, in fact, be coupled with
modulation types to characterize the transmission behavior of the licensed
users.

We provide one numerical example in Fig. 5 to shed lights on how to obtain the
coupled parameters from unsupervised learning approach. Higher order
statistics (HOS) are used as the feature for recognition because they
characterize the distribution shape of noisy baseband samples with low
complexity [12]. We introduce a new concept “Modulation Pattern” to denote the
combination of the modulation type and the transmit power level. We then
organize the estimates of cumulants of different orders and lags as a multiple
cumulant vector, serving as the input feature of the machine learning
algorithms. It can be proved that the feature vector is a multivariate
asymptotic Gaussian approximation of cumulants [14]. Specially, the first
cumulant in the vector should be the second-order cumulant in order to
indicate the energy information.

Figure 5: Pattern discovery and clustering results of multiple cumulants
vectors before normlization using Dirichlet Process Gaussian Mixture
Model.$\mu$ denotes the mean vector of each Gaussian component and $P_{i}$
denotes the transmit power level where
$P_{1}^{2}:P_{2}^{2}:P_{3}^{2}:P_{4}^{2}=2.5:5:5.3:4$ when overall average SNR
is 10 $dB$ and the number of samples is 100.

In order to construct the modulation dictionary and identify different
“modulation patterns”, we exploit Dirichlet Process Gaussian Mixture Model
(DPGMM), a type of unsupervised clustering analysis, to construct dynamic and
flexible statistical representations for the training data. We formulate a
mixture model in which the number of mixture components is infinite and is not
required at the beginning of clustering. After the convergence of DPGMM [15],
the cumulant vectors that aggregate around the origin of coordinate system
represent the “noise constellation”. Detection of cumulant vectors belonging
to “noise constellation” reveals that only noise signals are received and no
symbols are transmitted. The variance of noise can be evaluated from the
second-order cumulant of “noise constellation”. Furthermore, the average
vector of any other cluster is evaluated. The average vectors are used to
identify the alive modulation types and establish the minimal dictionary
following the correspondence of HOS and modulations. Consequently, new
observations of cumulant vector can be efficiently classified and predicted
using the predictive probability distribution constructed by DPGMM. The
dictionary can also be updated along with the update of posterior predictive
distributions adaptively.

The simulation results in Fig. 5 show the efficacy of the modulation
recognition and indicate superior discrimination capability compared with the
pure pattern recognition approach [12]. Moreover, Fig. 5 also demonstrates
that four modulation types as well as the channel idle status can be specified
automatically, along which the transmit power level can be cognized as well.
Furthermore, the modulation preference of certain user is understood from each
component of the predictive distribution.

###  III-C Example 3: Prediction for Spectrum Environment

In the last example, we evaluate the prediction of multi-channels spectrum
occupancy using machine learning techniques. Due to energy and hardware
constraints, unlicensed users may not be able to perform spectrum sensing at
all channels. This can be relieved by predicting the channel occupancies
before each sensing time slot starts. As modeled in [16], the state of
different channels and the time the licensed users spend dwelling on each
state are assumed to be independent. Moreover, the change of spectrum
environment can be considered as the process of license channel vacancy and
occupancy states appear alternately. We introduce throughput as the evaluation
criterion of the proposed learning-based cognitive strategy.

The multi-channels spectrum sensing based on learning and prediction can be
performed via three steps. Primarily, the parameters of the learning model are
estimated by the historical information embedded in the previous observations
database. Subsequently, given a new observation of the channel occupancy, the
next time slot can be predicted using the learning model. In the meanwhile the
probability of vacancy status of each channel can be ordered from high to low.
Lastly the unlicensed user is able to efficiently perform spectrum sensing in
accordance with the order and update the channel status database with the true
channel status detected. The average throughput can be derived in terms of the
mean error prediction probability and the probability of vacancy state for
objective channels in closed-form [16]. Simulation results in Fig. 6
demonstrate that the multi-channels spectrum sensing and prediction based on
machine learning techniques achieve better performance in terms of throughput,
compared with random selection of license channels.

Figure 6: The performance of mean throughput versus the probability of vacancy
state when the number of channels is $25$ and the normalized channel capacity
is $1$bit/s.

Besides the enhancement of spectrum sensing efficiency, there are a few other
benefits gained from learning and prediction of multiple channels occupancy.
Firstly, prediction results can also be utilized to foresee the traffic flow
of cellular networks consist of certain users at certain spectrum sub-bands.
Hence the network congestion may be detected in advance. Moreover, prediction
of channel occupancies makes it possible for predicting transmission behaviour
of users.

##  IV Summary and Prospects

In this article, we introduce a concrete multi-parameter cognitive
architecture for future wireless communication systems that contains three key
stages. Our efforts prompt a compromising but necessary way towards the
ultimate goal of CR, i.e., “full cognition”, which still stays in the
conceptual level. Examples demonstrate that the cognition in multiple
parameters space can be rather different from conventional single parameter
cognition and reveal fundamental insights in the proposed cognitive
architecture. The future research could include how to utilize the cognized
parameters to preserve quality of service of the network in terms of data
traffic, latency, overhead, etc.

## References

  * [1] S. Singh and P. Singh, “Key concepts and network architecture for 5G mobile technology,” Int. J. of Scientific Research Eng. & Tech. (IJSRET), vol. 1, no. 5, pp. 165–170, Aug. 2012.
  * [2] C.-X. Wang et al., “Cellular architecture and key technologies for 5G wireless communication networks,” IEEE Commun. Mag., vol. 52, no. 2, pp. 12–30, Feb. 2014.
  * [3] Z. Zhang, K. Long, and J. Wang, “Self-organization paradigms and optimization approaches for cognitive radio technologies: a survey,” IEEE Wireless Commun., vol. 20, no. 2, pp. 36–42, Apr. 2013.
  * [4] Z. Zhang, K. Long, and J. Wang, “On swarm intelligence inspired self-organized networking: its bionic mechanisms, designing principles and optimization approaches,” IEEE Commun. Surveys & Tutorials, vol. 16, no. 1, First Quarter, 2014.
  * [5] J. Mitola and J. Maguire, “Cognitive radio: making software radios more personal,” IEEE Personal Commun., vol. 6, no. 4, pp. 13–18, Aug. 1999.
  * [6] J. Li, F. Gao, T. Jiang, and W. Chen, “A new spectrum sensing strategy when primary user has multiple power levels,” to appear in IEEE GLOBECOM, 2014. Available at http://arxiv.org/abs/1311.5681.
  * [7] Y. Huang, H. Jiang, and H. Hu, “Design of learning engine based on support vector machine in cognitive radio,” in IEEE Int. Conf. on Computational Intelligence and Software Engineering, Wuhan, China, Dec. 2009, pp. 1–4.
  * [8] B. Lo and I. Akyildiz, “Reinforcement learning-based cooperative sensing in cognitive radio ad hoc networks,” in IEEE Int. Symp. on Personal Indoor and Mobile Radio Commun., Instanbul, Turkey, Sept. 2010, pp. 2244–2249.
  * [9] V. Tumuluru, P. Wang, and D. Niyato, “A neural network based spectrum prediction scheme for cognitive radio,” in IEEE Int. Conf. on. Commun., Capetown, Southafrica, May 2010, pp. 1–5.
  * [10] J. Guo, G. Zhong, and D. Qu, “Multi-slot spectrum sensing with backward SPRT in cognitive radio networks,” in IEEE Int. Conf. on Wireless Commun. & Signal Proc., Nanjing, China, Nov. 2009, pp. 1–5.
  * [11] _In the Matter of Unlicensed Operation in the TV Broadcast Bands_ , ET Docket No. 04–186, Notice of Proposed Rulemaking, FCC OET, May 2004.
  * [12] A. Swami and B. Sadler, “Hierarchical digital modulation classification using cumulants,” IEEE Trans. Commun., vol. 48, no. 3, pp. 416–429, Mar. 2000.
  * [13] O. A. Dobre, A. Abdi, Y. Bar-Ness, and W. Su, “Survey of automatic modulation classification techniques: classical approaches and new trends,” IET Commun., vol. 1, no. 2, pp. 137–156, Apr. 2007.
  * [14] A. Dandawate and G. Giannakis, “Asymptotic theory of mixed time averages and kth-order cyclic-moment and cumulant statistics,” IEEE Trans. Inf. Theory, vol. 41, no. 1, pp. 216–232, Jan. 1995.
  * [15] G. Rasmussen, “Dirichlet process Gaussian mixture models: choice of the base distribution,” J. of Comp. Science and Tech., vol. 25, no. 4, pp. 653–664, July 2010.
  * [16] Z. Ye, Q. Feng, and K. Shen, “Spectrum environment learning and prediction in cognitive radio,” in IEEE Int. Conf. on Signal Proc., Commun. and Comput., Xi’an, China, Sept. 2011, pp. 1–6.


# Topologically protected boundary discrete time crystal for a solvable model

Peng Xu  School of Physics, Zhengzhou University, Zhengzhou 450001, China
Tian-Shu Deng<EMAIL_ADDRESS>Institute for Advanced Study, Tsinghua
University, Beijing,100084, China

###### Abstract

Floquet time crystal, which breaks discrete time-translation symmetry, is an
intriguing phenomenon in non-equilibrium systems. It is crucial to understand
the rigidity and robustness of discrete time crystal (DTC) phases in a many-
body system, and finding a precisely solvable model can pave a way for
understanding of the DTC phase. Here, we propose and study a solvable spin
chain model by mapping it to a Floquet superconductor through the Jordan-
Wigner transformation. The phase diagrams of Floquet topological systems are
characterized by topological invariants and tell the existence of anomalous
edge states. The sub-harmonic oscillation, which is the typical signal of the
DTC, can be generated from such edge states and protected by topology. We also
examine the robustness of the DTC by adding symmetry-preserving and symmetry-
breaking perturbations. Our results on topologically protected DTC can provide
a deep understanding of the DTC when generalized to other interacting or
dissipative systems.

##  I Introduction

Periodic driving is a powerful tool to engineer variety of unique and
fascinating phenomena in many-body systems, such as Mott-insulator-superfluid
transitionZenesini  _et al._ (2009), dynamical gauge fieldXu  _et al._ (2021);
Clark  _et al._ (2018), many-body echoChen  _et al._ (2020a, b); Lv  _et al._
(2020), the realization of topological band structuresJotzu  _et al._ (2014);
Aidelsburger  _et al._ (2015). Recently, Floquet systems also bring new
possibilities to simulate time-translation symmetry broken phases, which are
called discrete time crystals (DTC), and this has attracted much attention
from both experimentalists Zhang  _et al._ (2017); Choi  _et al._ (2017);
Rovny  _et al._ (2018a, b); Pal  _et al._ (2018); Keßler  _et al._ (2021);
Randall  _et al._ (2021); Kyprianidis  _et al._ (2021); Bluvstein  _et al._
(2021); Mi  _et al._ (2022); Zhang  _et al._ (2022); Frey and Rachel (2022)
and theoristsKhemani  _et al._ (2016a); Else  _et al._ (2016a); Yao  _et al._
(2017); Huang  _et al._ (2018); Gong  _et al._ (2018); Zhu  _et al._ (2019);
Chinzei and Ikeda (2020); Huang  _et al._ (2022). The concept of the time
crystal was originally proposed by WilczekWilczek (2012). However, it was
ruled out by the no-go theorem in thermal equilibrium systemsBruno (2013);
Watanabe and Oshikawa (2015). Then Shivaji Sondhi’s Khemani  _et al._ (2016b)
and Chetan Nayak’s Else  _et al._ (2016b) groups generalized the concept of
time crystal and proposed the DTC, which exhibits a unique property that the
expectation values of generic observables manifest a sub-harmonic oscillation.
For example, the kicked Ising chain model with disordered interaction, where
spins collectively flip after one period and back to their initial state after
two periods, is a canonical realization of the DTC.

In a non-interacting spin chain system, taking
$\hat{U}=\exp{(-i\theta\sum_{j}\hat{X}_{j})}$ with $\theta=\pi/2$ as a Floquet
evolution operator is a straightforward method to flip all spins in one
period. Here, $\hat{X}_{j}$ is a spin operator acting on site $j$. However,
when $\theta$ deviates slightly from $\pi/2$, the period of observables also
deviates from twice of the Floquet period. It means the sub-harmonic response
induced by $\hat{U}=\exp{(-i\theta\sum_{j}\hat{X}_{j})}$ is a fine-tuned
result and easily destroyed by perturbations. This simple example implies that
the robustness of sub-harmonic response is a crucial property of the DTC.
According to the previous results, many-body localization and pre-
thermalization may provide two mechanisms to stabilize the sub-harmonic
response. Besides, topologically protected anomalous edge states Rudner  _et
al._ (2013); Titum  _et al._ (2016); Zhou and Gong (2018) also suggest another
mechanism of generating the robust DTC and the reasons are listed in the
following. Firstly, edge states in topological insulators (superconductors)
are protected by symmetries. As long as the symmetries are not broken and the
gap does not close, topological edge states are stable and robust. Secondly,
Floquet topological insulators (superconductors) host anomalous edge states
with quasi-energy $\pi/T$ which can generate a sub-harmonic response to
driving frequency $2\pi/T$. Although the relation between $\pi$-mode edge
states in Floquet topological system and the DTC has been discussed in
previous researchBomantara and Gong (2018); Pan and Wang (2020); Chew  _et
al._ (2020); Shtanko and Movassagh (2020); Bomantara  _et al._ (2021); Wang
_et al._ (2022), detailed and systematic of topologically protected DTC in a
more general spin chain model are still lacking. Bridging Floquet topological
superconductors to the topologically protected DTC in a general periodic
driven spin chain model and explicitly analyzing the dynamics of observables
are significantly helpful for deeply understanding the existence and
robustness of the DTC.

In this work, we investigate the existence of the DTC phase in a general
Floquet spin chain model. Such a periodically driven spin chain can be mapped
to a Floquet superconductor through the Jordan-Wigner transformation, after
which the model becomes the form of Majorana fermion. This model is intrinsic
with particle-hole symmetry. Furthermore, this system can be classified into D
class or BDI class, which is dependent on whether chiral symmetry is
preserved. The topological Floquet superconductor exhibit a special kind of
topologically non-trivial phase, where anomalous edge states with quasi-energy
$\pi/T$ exist. In order to observe a robust DTC, the observable should be
selected as the anomalous edge mode or the end spin. We numerically
demonstrate that both observables manifest a sub-harmonic response with a
generic initial product state. Finally, we also confirm the robustness of the
DTC by adding symmetry-preserving and symmetry-breaking perturbations.

The paper is organized as follows. In section II, we describe a general
periodically driven spin chain model and map it to a Floquet superconductor.
In section III, we discuss the topological classification of this model,
calculate the topological invariants and obtain the phase diagrams. In section
IV, we demonstrate the existence of the DTC resulting from the Floquet
superconductors, by selecting the observable as the anomalous edge mode or the
end spin. Finally, in section V, we examine the robustness of the DTC is
protected by topologically non-trivial phase by adding symmetry-preserving and
symmetry-breaking perturbations.

##  II Model: Periodically Driven Spin Chain

We consider a periodically driven spin-$1/2$ chain as illustrated in Fig. 1
(a),

$\displaystyle\hat{H}(t)=\begin{cases}\hat{H}_{1},&\text{for}\quad nT\leqslant
t<nT+t_{1}\\\ \hat{H}_{2},&\text{for}\quad nT+t_{1}\leqslant
t<(n+1)T\end{cases},$ (1)

where

$\displaystyle\hat{H}_{m}=J_{m}^{xx}\sum_{j}\hat{X}_{j}\hat{X}_{j+1}+J_{m}^{yy}\sum_{j}\hat{Y}_{j}\hat{Y}_{j+1}$
(2)
$\displaystyle+J_{m}^{xy}\sum_{j}\hat{X}{}_{j}\hat{Y}_{j+1}+J_{m}^{yx}\sum_{j}\hat{Y}_{j}\hat{X}_{j+1}+h_{m}^{z}\sum_{j}\hat{Z}{}_{j}.$

$\hat{X}_{j}$, $\hat{Y}_{j}$ and $\hat{Z}_{j}$ are spin operators (with the
form of Pauli matrices) acting on the $j$-th site and the chain contains $N$
spins. $J_{m}^{xx}$, $J_{m}^{yy}$, $J_{m}^{xy}$, $J_{m}^{yx}$ represent the
strengths of nearest spin interactions and $h_{m}^{z}$ the transverse field
during $m$-th time interval. By employing the Jordan-Wigner transformation

$\displaystyle\hat{X}_{j}$
$\displaystyle=(\hat{c}_{j}^{\dagger}+\hat{c}_{j})e^{i\pi\sum_{l<j}\hat{c}_{l}^{\dagger}\hat{c}_{l}},$
(3) $\displaystyle\hat{Y}_{j}$
$\displaystyle=-i(\hat{c}_{j}^{\dagger}-\hat{c}_{j})e^{i\pi\sum_{l<j}\hat{c}_{l}^{\dagger}\hat{c}_{l}},$
$\displaystyle\hat{Z}_{j}$
$\displaystyle=2\hat{c}_{j}^{\dagger}\hat{c}_{j}-1.$

where $\hat{c}_{j}^{\dagger}$ and $\hat{c}_{j}$ are the creation and
annihilation fermion operators on $j$-th site. Then the Hamiltonian in Eq.
(LABEL:HiSpin) can be mapped to a fermionic system,

$\displaystyle\hat{H}_{m}$
$\displaystyle=\sum_{j}\left(J_{m}^{xx}-J_{m}^{yy}-iJ_{m}^{xy}-iJ_{m}^{yx}\right)\hat{c}_{j}^{\dagger}\hat{c}_{j+1}^{\dagger}$
(4)
$\displaystyle+\sum_{j}\left(J_{m}^{xx}+J_{m}^{yy}+iJ_{m}^{xy}-iJ_{m}^{yx}\right)\hat{c}_{j}^{\dagger}\hat{c}_{j+1}$
$\displaystyle+\frac{h_{m}^{z}}{2}\sum_{j}(2\hat{c}_{j}^{\dagger}\hat{c}_{j}-1)+h.c..$

Figure 1: (a) Schematic of spin chain. (b) Schematic of Majorana chain. (c)
Quasi-energy spectrum $\exp(-i\epsilon_{\mu}T)$ with $\epsilon_{\mu}$ the
quasi-energy excitation. $0$ and $\pi/T$ represent two edge states quasi-
energy excitation. (d) Eigenvalues of the Floquet evolution operator
$\hat{U}_{T}$. Quasi-particle operator $\hat{\beta}^{\dagger}_{\pi}$ excites a
state $|E\rangle$ with eigenenergy $E$ to another state $|E+\pi/T\rangle$.

Furthermore, by defining
$\hat{\gamma}_{j,1}=\frac{1}{\sqrt{2}}\left(\hat{c}_{j}^{\dagger}+\hat{c}_{j}\right)$
and
$\hat{\gamma}_{j,2}=i\frac{1}{\sqrt{2}}\left(\hat{c}_{j}-\hat{c}_{j}^{\dagger}\right)$,
we obtain a Hamiltonian in Majorana representaion, as illustrated in Fig. 1
(b),

$\displaystyle\hat{H}_{m}$
$\displaystyle=2J_{m}^{xx}\sum_{j}i\hat{\gamma}_{j,2}\hat{\gamma}_{j+1,1}-2J_{m}^{yy}\sum_{j}i\hat{\gamma}_{j,1}\hat{\gamma}_{j+1,2}$
(5)
$\displaystyle+2J_{m}^{xy}\sum_{j}i\hat{\gamma}_{j,2}\hat{\gamma}_{j+1,2}-2J_{m}^{yx}\sum_{j}i\hat{\gamma}_{j,1}\hat{\gamma}_{j+1,1}$
$\displaystyle+2J_{m}^{z}\sum_{j}i\hat{\gamma}_{j,2}\hat{\gamma}_{j,1}.$

Both Hamiltonians $\hat{H}_{1}$ and $\hat{H}_{2}$ in Eq. (5) are quadratic, so
they can be rewritten as

$\hat{H}_{m}=\frac{1}{2}\hat{\Psi}^{\dagger}H_{m}\hat{\Psi},$ (6)

with

$\hat{\Psi}^{\dagger}=\left(\begin{array}[]{ccccc}\hat{\gamma}_{1,1}&\hat{\gamma}_{1,2}&\hat{\gamma}_{2,1}&\hat{\gamma}_{2,2}&...\end{array}\right),$
(7)

where $H_{m}$ is a $2N*2N$ anti-symmetric matrix. Here and after, $\hat{H}$
($\hat{U}$) represents an operator and $H$ ($U$) represents a matrix. Applying
Fourier transformation
$\hat{\gamma}_{j,m}=\frac{1}{\sqrt{N}}\sum_{k}\hat{\gamma}_{k,m}e^{ikj}$,
$\hat{H}_{m}(k)$ has such a form

$\displaystyle\hat{H}_{m}=$
$\displaystyle\sum_{k}\left(\begin{array}[]{cc}\hat{\gamma}_{k,1}^{\dagger}&\hat{\gamma}_{k,2}^{\dagger}\end{array}\right)\left[d^{0}_{m}(k)I+\bm{d}_{m}(k)\cdot\bm{\sigma}\right]\left(\begin{array}[]{c}\hat{\gamma}_{k,1}\\\
\hat{\gamma}_{k,2}\end{array}\right),$ (11)

where $\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ denote Pauli matrices.
$d^{0}_{m}(k)$ and $\bm{d}_{m}(k)$ are given by

$\displaystyle d_{m}^{0}(k)$
$\displaystyle=\left(J_{m}^{yx}-J_{m}^{xy}\right)\sin(k),$ (12) $\displaystyle
d_{m}^{x}(k)$ $\displaystyle=\left(J_{m}^{yy}-J_{m}^{xx}\right)\sin(k),$
$\displaystyle d_{m}^{y}(k)$
$\displaystyle=J_{m}^{z}+\left(J_{m}^{xx}+J_{m}^{yy}\right)\cos(k),$
$\displaystyle d_{m}^{z}(k)$
$\displaystyle=\left(J_{m}^{xy}+J_{m}^{yx}\right)\sin(k).$

For the periodically driven Hamiltonian in Eq. (1), the Floquet evolution
operator $\hat{U}_{T}$ and the effective Hamiltonian $\hat{H}_{\text{eff}}$
are defined by

$\displaystyle\hat{U}_{T}$
$\displaystyle=\hat{\mathcal{T}}\exp\left(-i\int_{0}^{T}\hat{H}(t)dt\right)=\exp(-i\hat{H}_{2}t_{2})\exp(-i\hat{H}_{1}t_{1})$
$\displaystyle=\exp(-i\hat{H}_{\rm{eff}}T),$ (13)

where $t_{2}+t_{1}=T$. Based on the quadratic form of $\hat{H}_{1}$ and
$\hat{H}_{2}$ in Eq. (5), we have proven in Appendix A that the Floquet
effective Hamiltonian is also quadratic, which means $\hat{H}_{\rm{eff}}$ also
takes the form

$\hat{H}_{\rm{eff}}=\frac{1}{2}\hat{\Psi}^{\dagger}H_{\rm{eff}}\hat{\Psi}.$
(14)

In order to obtain the eigenenergies of $\hat{H}_{\rm{eff}}$, we diagonalize
the above matrix $H_{\rm{eff}}$ in real space as

$H_{{\rm eff}}=V\Lambda^{\epsilon}V^{\dagger},$ (15)

where

$V=\left(\begin{array}[]{cccc}|V_{1}\rangle&|V_{2}\rangle&...&|V_{2N}\rangle\end{array}\right),$
(16) $\Lambda=\left(\begin{array}[]{cccc}\epsilon_{1}\\\ &\epsilon_{2}\\\
&&\ddots\\\ &&&\epsilon_{2N}\end{array}\right).$ (17)

$|V_{\mu}\rangle$ is the eigenvector of $H_{\rm{eff}}$ with eigenvalue
$\epsilon_{\mu}$ such that $H_{{\rm
eff}}|V_{\mu}\rangle=\epsilon_{\mu}|V_{\mu}\rangle$. $\epsilon_{\mu}$ is the
Bogoliubov quasi-particle excitation spectrum of $\hat{H}_{\rm{eff}}$, which
is restricted to the Floquet Brillouin Zone $[-\pi/T,\pi/T)$. Then, the
Floquet effective Hamiltonian $\hat{H}_{\rm{eff}}$ in this quasi-particle
representation can be rewritten as

$\hat{H}_{{\rm
eff}}=\frac{1}{2}\hat{\Psi}^{\dagger}V\Lambda^{\epsilon}V^{\dagger}\hat{\Psi}=\frac{1}{2}\hat{\Phi}^{\dagger}\Lambda^{\epsilon}\hat{\Phi},$
(18)

where $\hat{\Phi}$ is composed of Bogoliubov quasi-particles
$\hat{\Phi}^{\dagger}=\left(\begin{array}[]{cccc}\hat{\alpha}_{1}^{\dagger}&\hat{\alpha}_{2}^{\dagger}&\cdots&\hat{\alpha}_{2N}^{\dagger}\end{array}\right)$
and $\hat{\alpha}^{\dagger}_{\mu}=\hat{\Psi}^{\dagger}|V_{\mu}\rangle$.

Due to the particle-hole symmetry constraint, the eigenvalues and eigenvectors
of $H_{\rm{eff}}$ must come in pair

$\displaystyle H_{\rm{eff}}|V_{\mu}\rangle$
$\displaystyle=\epsilon_{\mu}|V_{\mu}\rangle,$ (19) $\displaystyle
H_{\rm{eff}}|V_{\mu}^{\ast}\rangle$
$\displaystyle=-\epsilon_{\mu}|V_{\mu}^{\ast}\rangle.$

Substituting these relations into Eq. (18), the effective Hamiltonian
$\hat{H}_{\rm{eff}}$ becomes

$\hat{H}_{{\rm
eff}}=\sum_{\epsilon_{\mu}>0}\epsilon_{\mu}(\hat{\alpha}_{\mu}^{\dagger}\hat{\alpha}_{\mu}-\frac{1}{2}).$
(20)

The eigenstate of $\hat{H}_{\rm{eff}}$ is given by
$\prod_{\mu}|n_{\mu}\rangle=(\alpha_{\mu}^{\dagger})^{n_{\mu}}|0\rangle$.
$n_{\mu}=0,1$ represents the number of occupations on quasi-particle mode.
Then the energy spectrum of this system can be expressed as

$\displaystyle E$ $\displaystyle=\langle\hat{H}_{{\rm eff}}\rangle$
$\displaystyle=\sum_{\epsilon_{\mu}>0}\epsilon_{\mu}(\langle\alpha_{\mu}^{\dagger}\alpha_{\mu}\rangle-\frac{1}{2})$
$\displaystyle=E_{0}+\sum_{\mu}\epsilon_{\mu}{n}_{\mu},$ (21)

where $E_{0}=-\frac{1}{2}\sum_{\epsilon_{\mu}>0}\epsilon_{\mu}$ is the energy
of ground state and different configurations of ${n}_{\mu}$ generate $2^{N}$
eigenenergies of the spin chain system.

To arrive at our final destination of observing the DTC in this periodically
driven spin chain, we firstly concentrate on the evolution of operator
$\hat{\alpha}^{\dagger}_{\mu}$

$\displaystyle\hat{\alpha}^{\dagger}_{\mu}(nT)$
$\displaystyle=\hat{U}_{T}^{-n}\hat{\alpha}^{\dagger}_{\mu}(0)\hat{U}_{T}^{n}=\hat{\alpha}^{\dagger}_{\mu}(0)e^{inT\epsilon_{\mu}}.$
(22)

For a given initial state, the expectation value of the operator
$\hat{\alpha}^{\dagger}_{\mu}$ oscillates with period
${2\pi}/{\epsilon_{\mu}}$. Generally, this oscillating behavior is not rigid
and easily perturbed, so the performance for operators
$\hat{\alpha}^{\dagger}_{\mu}(nT)$ cannot be regarded as a signal of the DTC
phase. However, there is an exception due to anomalous edge states in Floquet
topological systems. Such systems exhibit two kinds of edge states with
eigenenergies $\epsilon_{\mu}=0$ and $\epsilon_{\mu}=\pi/T$, protected by
topological invariants $\nu_{0}$ and $\nu_{\pi}$, respectively. Here, we
mainly focus on the edge states with eigenenergy $\epsilon_{\mu}=\pi/T$,
denoted as $H_{\rm{eff}}|W\rangle=\pi/T|W\rangle$ and we define
$\hat{\beta}^{\dagger}_{\pi}=\hat{\Psi}^{\dagger}|W\rangle$. The expectation
value of the operator $\hat{\beta}^{\dagger}_{\pi}$ exactly oscillates with
period $2T$ as shown in Eq. (22). More importantly, this oscillation period is
protected by the Floquet topologically non-trivial phases such that this
phenomenon of the DTC is robust against symmetry-preserving perturbations.

To further demonstrate the DTC induced by Floquet topological superconductors,
we systematically discuss the phase diagrams of Floquet topological
superconductors in section III and discuss the DTC for different observables,
including anomalous edge operators and end spin operators in section IV.

##  III Floquet topological phases

Following the standard approach of studying topological superconductors, we
calculate the topological invariants of the effective Hamiltonian in momentum
space. $\hat{H}_{\rm{eff}}$ is expressed as

$\hat{H}_{{\rm
eff}}=\sum_{k}\left(\begin{array}[]{cc}\hat{\gamma}_{k,1}&\hat{\gamma}_{k,2}\end{array}\right)H_{\rm{eff}}(k)\left(\begin{array}[]{c}\hat{\gamma}_{k,1}\\\
\hat{\gamma}_{k,2}\end{array}\right),$ (23)

where $H_{\rm{eff}}(k)$ are obtained from

$\displaystyle U(k)$
$\displaystyle=\exp[-iH_{2}(k)t_{2}]\exp[-iH_{1}(k)t_{1}]$ (24)
$\displaystyle=\exp[-iH_{{\rm eff}}(k)T]$

according to Appendix A. Here, $d_{{\rm eff}}^{0}(k)$ and
$\boldsymbol{d}_{{\rm eff}}(k)$ are defined using $H_{{\rm eff}}(k)=d_{{\rm
eff}}^{0}(k)I+\boldsymbol{d}_{{\rm eff}}(k)\cdot\boldsymbol{\sigma}$. In the
following, we request $\epsilon(k)=d_{{\rm eff}}^{0}(k)\pm\sqrt{d_{{\rm
eff}}^{2}(k)}\neq 0$ or $\pi/T$ for any $k\in[-\pi,\pi)$, due to the
topological superconductors should be gapped.

For a one-dimensional two-band system with intrinsic particle-hole constraint,
topologically non-trivial phases can only be classified into classes BDI (with
chiral symmetry) and D (without chiral symmetry) foo ; Chiu  _et al._ (2016).
So, in the following, we discuss the calculations of topological invariants
for classes D and BDI, respectively.

###  III.1 D Class

For simplicity, we take $J_{1}^{yy}=J_{2}^{yy}=0,J_{2}^{xx}=h_{1}^{z}=0$, non-
zero coupling strengths $J_{1}^{xy}=J_{2}^{xy}=J_{1}^{yx}=J_{2}^{yx}$ to break
chiral symmetry. Besides, here and after, we take $t_{1}=t_{2}=t$. Then
$H_{1}(k)$ and $H_{2}(k)$ is given by

$\displaystyle
H_{1}(k)=-J_{1}^{xx}(\sin(k)\sigma_{x}-\cos(k)\sigma_{y})+(J_{1}^{xy}+J_{1}^{yx})\sin(k)\sigma_{z},$
$\displaystyle
H_{2}(k)=h_{2}^{z}\sigma_{y}+(J_{1}^{xy}+J_{1}^{yx})\sin(k)\sigma_{z}.$ (25)

Following the approach in Ref.Jiang  _et al._ (2011), topological invariants
$\nu_{0}$ and $\nu_{\pi}$ in Floquet superconductors have such forms

$\displaystyle\nu_{0}\nu_{\pi}$
$\displaystyle=\operatorname{sgn}\left[{\rm{Pf}}(M_{0})\right]\operatorname{sgn}\left[{\rm{Pf}}(M_{\pi})\right],$
(26) $\displaystyle\nu_{0}$
$\displaystyle=\operatorname{sgn}\left[{\rm{Pf}}(N_{0})\right]\operatorname{sgn}\left[{\rm{Pf}}(N_{\pi})\right],$

where $\operatorname{Pf}\left[X\right]$ is the Pfaffian number of a skew marix
$X$ and $M_{k}=\log\left[U(k)\right],N_{k}=\log\left[\sqrt{U(k)}\right]$.
Since both $\log({X})$ and $\log(\sqrt{X})$ are multivalued functions, for
$k=0$ and $k=\pi$ we have

$M_{k}=-i\xi_{M}(k)\sigma_{y},\quad N_{k}=-i\xi_{N}(k)\sigma_{y},$ (27)

with

$\displaystyle\xi_{M}(k)+2z\pi$ $\displaystyle=J_{1}^{xx}t\cos(k)+h_{2}^{z}t,$
(28) $\displaystyle\xi_{N}(k)+2z\pi$
$\displaystyle=\frac{J_{1}^{xx}t\cos(k)+h_{2}^{z}t}{2}.$

Here, $z$ is taken as an appropriate integer so that $\xi_{M/N}(k)$ can be
constrained to the interval $[-\pi,\pi)$. Finally, we obtain

$\displaystyle\nu_{0}\nu_{\pi}$ $\displaystyle={\rm
sgn}[\sin(P_{+})\sin(P_{-})],$ (29) $\displaystyle\nu_{0}$ $\displaystyle={\rm
sgn}[\sin(\frac{P_{+}}{2})\sin(\frac{P_{-}}{2})].$

with

$\displaystyle P_{+}$ $\displaystyle=h_{2}^{z}t+J_{1}^{xx}t,$ (30)
$\displaystyle P_{-}$ $\displaystyle=h_{2}^{z}t-J_{1}^{xx}t.$

In Fig. 2 (a, b), we show the phase diagrams of topological invariants
$\nu_{0}$ and $\nu_{\pi}$ on the $\theta_{1}$-$\theta_{2}$ plane by taking
$\theta_{1}=J_{1}^{xx}t$ and $\theta_{2}=h_{2}^{z}t$. For the D class, the
system is classified by a $\mathbb{Z}_{2}$ topological invariant. $\nu_{0}$
and $\nu_{\pi}$ take the value $\pm 1$, and $\nu_{\pi}=-1$ indicates anomalous
topologically non-trivial phases. The quasi-energy spectrum as a function of
$\theta_{1}\in[-\pi,\pi)$ is plotted in Fig. 2 (c). The gapless $0$-mode and
$\pi$-mode must exist where the topological invariants $\nu_{0}=-1$ and
$\nu_{\pi}=-1$, respectively. Consequently, the topologically protected
discrete time crystal could be observed in the regime $\nu_{\pi}=-1$.

Furthermore, because of particle-hole symmetry, the eigenvalues and
eigenvectors of $H_{\rm{eff}}$ always come in pair as shown in Eq. (19). Thus
for the edge states with quasi-energies $\epsilon_{\mu}=0$ or $\pi/T$, the
eigenvectors $|V_{\mu}\rangle$ and $|V_{\mu}^{*}\rangle$ span a degenerate
subspace, from which we can construct a pure real vector and a pure imaginary
vector as $|V_{\mu}^{R}\rangle=|V_{\mu}\rangle+|V_{\mu}^{*}\rangle$ and
$|V_{\mu}^{I}\rangle=|V_{\mu}\rangle-|V_{\mu}^{*}\rangle$. In Fig. 2 (d), we
have shown the probability distribution of an anomalous edge state
$|W^{R}\rangle=(|W\rangle+|W^{*}\rangle)/\sqrt{2}$. For convenience, here and
after, $\\{j,s\\}$ ($j=1,2,...L$ and $s=1,2$) in Eq. (7) is relabelled as
$\\{2j+s-2\\}$ for the index of $\hat{\gamma}$, as well as $|W^{R}\rangle$.

Figure 2: (a, b) Phase diagrams of the D class topological superconductor
(without chiral symmetry). The areas of topological invariants with value $-1$
indicate a topologically non-trivial superconductor, otherwise topologically
trivial. Topological invariant $\nu_{0}$ is shown in (a) and $\nu_{\pi}$ in
(b). (c) Quasi-energy spectrum of the topological superconductor under open
boundary condition with parameters $\theta_{1}\in[-\pi,\pi),\theta_{2}=3\pi/4$
along the black solid line shown in (a, b). (d) $\pi$-mode edge state
$|W^{R}\rangle$ with parameters $(\theta_{1},\theta_{2})=(\pi/2,3\pi/4)$.
Here, we set $J_{1}^{xx}t=\theta_{1},h_{2}^{z}t=\theta_{2}$.

###  III.2 BDI Class

For simplicity, we take non-zero coupling strengths $J_{1}^{xx}=J_{2}^{xx}$,
$h_{2}^{z}$ and other parameters zero so that the system obeys chiral
symmetry. With given parameters, the chiral operator is taken as
$U_{S}=\sigma_{z}$ and we have

$\displaystyle U_{S}^{\dagger}U^{a,b}U_{S}=(U^{a,b})^{\dagger}.$ (31)

where

$\displaystyle{U}^{a}$
$\displaystyle=\exp(-iH_{1}t_{1}/2)\exp(-iH_{2}t_{2})\exp(-iH_{1}t_{1}/2)$
(32) $\displaystyle=\exp({-iH_{\rm{eff}}^{a}}T),$ $\displaystyle U^{b}$
$\displaystyle=\exp(-iH_{2}t_{2}/2)\exp(-iH_{1}t_{1})\exp(-iH_{2}t_{2}/2)$
$\displaystyle=\exp({-iH_{\rm{eff}}^{b}}T),$

and

$\displaystyle H_{1}(k)=-J_{1}^{xx}(\sin(k)\sigma_{x}-\cos(k)\sigma_{y}),$
(33) $\displaystyle
H_{2}(k)=h_{2}^{z}\sigma_{y}-J_{1}^{xx}(\sin(k)\sigma_{x}-\cos(k)\sigma_{y}).$

Here, we introduce $U^{a}$ and $U^{b}$, which are both unitarily transformed
from $U(k)$, in order that the chiral operator is explicitly $\sigma_{z}$.
Besides, based on $U^{a}$ and $U^{b}$ the topological invariants $\nu_{0}$ and
$\nu_{\pi}$ can be calculated as Asbóth and Obuse (2013) follows,

$\nu_{0,\pi}=\frac{\nu_{a}\pm\nu_{b}}{2},$ (34)

where $\nu_{a}$ and $\nu_{b}$ are the winding numbers defined as

$\nu_{a,b}=\frac{1}{2\pi}\int
dk\frac{-d_{a,b}^{x}(k)\partial_{k}d_{a,b}^{y}(k)+d^{y}_{a,b}(k)\partial_{k}d^{x}_{a,b}(k)}{\left(d_{a,b}^{x}(k)\right)^{2}+\left(d_{a,b}^{y}(k)\right)^{2}},$
(35)

and $d_{a,b}^{x,y}(k)$ is given by

$\displaystyle H_{\rm{eff}}^{a,b}(k)$
$\displaystyle=d_{a,b}^{x}(k)\sigma_{x}+d_{a,b}^{y}(k)\sigma_{y}.$ (36)

In Fig. 3 (a) and (b), we show the phase diagrams of topological invariants
$\nu_{0}$ and $\nu_{\pi}$ on the $\theta_{1}$-$\theta_{2}$ plane. For the BDI
class, the system is classified by a $\mathbb{Z}$ topological invariant. Both
$\nu_{0}$ and $\nu_{\pi}$ take integer values, and $\nu_{\pi}\neq 0$ indicates
anomalous topologically non-trivial phases. The quasi-energy spectrum as a
function of $\theta_{1}\in[-\pi,\pi)$ is plotted in Fig. 2 (c). In contrast to
D class, the topological invariant in class BDI $\nu_{0}$ and $\nu_{\pi}$
could be $\pm 2,\pm 3,\cdots$ leading to larger degeneracy of edge states. For
example, for $(\theta_{1},\theta_{2})=(3\pi/8,7\pi/8)$, we have $\nu_{\pi}=2$,
which results in four $\pi$-mode edge states. These four edge states are
labeled as $|W_{\eta}\rangle$ and $|W_{\eta}^{*}\rangle$ with $\eta=1,2$. In
Fig. 3 (d), we show the the probability distribution of the anomalous edge
states as $|W^{R}_{\eta=1}\rangle$ and $|W^{R}_{\eta=2}\rangle$, respectively.

Figure 3: (a, b) Phase diagrams of the BDI class topological superconductor
(with chiral symmetry). The areas of topological invariants with non-zero
values indicate a topologically non-trivial superconductor, otherwise
topologically trivial. Topological invariant $\nu_{0}$ is shown in (a) and
$\nu_{\pi}$ in (b). (c) Quasi-energy spectrum of the topological
superconductor under open boundary condition with parameters
$\theta_{1}\in[-\pi,\pi),\theta_{2}=7\pi/8$ along the black solid line shown
in (a, b). (d) $\pi$-mode edge state $|W^{R}\rangle$ with parameters
$(\theta_{1},\theta_{2})=(3\pi/8,7\pi/8)$. Here, we set
$\theta_{1}=J_{1}^{xx}t=J_{2}^{xx}t,\theta_{2}=h_{2}^{z}t$.

##  IV Sub-harmonic oscillations of the DTC

Figure 4:  (a, b) Dynamics of anomalous edge mode $\hat{\beta}_{\pi}^{R}$. (c,
d) Distribution of response function $A_{\beta}(\omega)$. The frequency
distributions in (c) and (d) are the Fourier transformation of dynamics in (a)
and (b), respectively. (a, c) with the same parameters as Fig. 2 (d)
correspond to D class, while (b, d) with the same parameters as Fig. 3 (d)
correspond to BDI class. The initial state is taken as a product state
polarized along $x$ direction. Figure 5: (a, b) Upper panels: the values of
topological invariants in different regimes. Lower panels: frequency
distributions of the response function for observable $\hat{X}_{1}$ and the
color represents the amplitude of $A_{X}(\omega)$. Parameters of (a) are the
same as the black solid line in Fig. 2, while parameters of (b) are the same
as the black solid line in Fig. 3.

In this section, we’ll examine the sub-harmonic oscillations of the DTC in
topologically non-trivial phases with $\nu_{\pi}=-1$ in the D class and
$\nu_{\pi}\neq 0$ in the BDI class, by selecting the anomalous edge mode and
the end spin as observables.

We first take the anomalous edge mode as the observable for the DTC,

$\hat{\beta}_{\pi}^{R}=\frac{\hat{\beta}_{\pi}+\hat{\beta}_{\pi}^{\dagger}}{\sqrt{2}}.$
(37)

As we have derived in Eq. (22), the time evolution of observable
$\hat{\beta}_{\pi}^{R}$ is expressed as

$\displaystyle\langle\hat{\beta}_{\pi}^{R}(nT)\rangle$
$\displaystyle=\langle\psi(0)|\hat{U}_{T}^{-n}\hat{\beta}_{\pi}^{R}\hat{U}_{T}^{n}|\psi(0)\rangle$
(38)
$\displaystyle=(-1)^{n}\langle\psi(0)|\hat{\beta}_{\pi}^{R}|\psi(0)\rangle.$

It is obvious that the dynamics of $\langle\hat{\beta}_{\pi}^{R}(nT)\rangle$
manifests an exact sub-harmonic response for a generic initial state with
$\langle\psi(0)|\hat{\beta}_{\pi}^{R}|\psi(0)\rangle\neq 0$. After the Fourier
transformation of $\langle\hat{\beta}_{\pi}^{R}(nT)\rangle$, we obtain the
response function

$A_{\beta}(\omega)=\frac{1}{n_{\rm{max}}}\sum_{n=1}^{n_{\rm{max}}}\langle\hat{\beta}_{\pi}^{R}(nT)\rangle
e^{i\omega nT}.$ (39)

where $n_{\text{max}}$ represents the maximal number of periods. In Fig. 4 (a)
and (b), we show the expectation values of $\hat{\beta}^{R}_{\pi}$ as a
function of evolution time for D class and BDI class, respectively.
Correspondingly the response functions $A_{\beta}(\omega)$ are plotted in Fig.
4 (c) and (d). $\langle\hat{\beta}^{R}_{\pi}\rangle$ oscillates with double
periods and $A_{\beta}(\omega)$ peaks at $\omega=\pi/T$, which confirm sub-
harmonic response for the observable $\hat{\beta}^{R}_{\pi}$ in the regime of
Floquet topologically non-trivial phases with anomalous edge states. However,
we have to admit that $\hat{\beta}^{R}_{\pi}$ is difficult to engineer and
observe in experiments, and the concrete form of operator
$\hat{\beta}^{R}_{\pi}$ in spin basis depends on coupling parameters.
Therefore, it is important to choose a more appropriate operator that is
easily observed in experiments.

Another observable for the DTC is taken as the end spin operator
$\hat{X}_{1}$, since the anomalous edge states finitely occupy the end sites.
Actually, the end spin $\hat{X}_{1}$ is nothing but $\hat{\gamma}_{1}$ after
the Jordan-Wigner transformation, which could be decomposed as

$\hat{\gamma}_{1}=\sum_{\eta=1}^{n_{\rm{edge}}}|W_{\eta}\rangle_{1}\hat{\beta}_{\pi,\eta}+\sum_{\epsilon_{\mu}\neq\pi/T}|V_{\mu}\rangle_{1}\hat{\alpha}_{\mu}+h.c.,$
(40)

where $|W_{\eta}\rangle_{1}$ and $|V_{\mu}\rangle_{1}$ represent the first
element of vectors $|W_{\eta}\rangle$ and $|V_{\mu}\rangle$, respectively.
$n_{\rm edge}$ represents the number of pairs of edge states. Dynamics of the
observable $\hat{X}_{1}$ is then given by

$\displaystyle\langle\hat{X}_{1}(nT)\rangle$
$\displaystyle=\langle\psi(0)|\hat{U}_{T}^{-n}\hat{\gamma}_{1}\hat{U}_{T}^{n}|\psi(0)\rangle$
(41)
$\displaystyle=(-1)^{n}\sum_{\eta=1}^{n_{\rm{edge}}}2{\rm{Re}}(|W_{\eta}\rangle_{1}\langle\psi(0)|\hat{\beta}_{\pi,\eta}|\psi(0)\rangle)$
$\displaystyle+\sum_{\epsilon_{\mu}\neq\pi/T}2{\rm
Re}(|V_{\mu}\rangle_{1}e^{2i\epsilon_{\mu}}\langle\psi(0)|\hat{\alpha}_{\mu}|\psi(0)\rangle).$

In D class, there is only one pair of edge states with $n_{\rm{edge}}=1$,
which is localized around $\hat{\gamma}_{1}|0\rangle$ and
$\hat{\gamma}_{2N}|0\rangle$. Therefore, $|W_{1}\rangle_{1}$ is finitely
distributed so that $\langle\hat{X}_{1}(nT)\rangle$ manifests a sub-harmonic
response at $\omega=\pi/T$ as is shown in the second line of Eq. (41). In Fig.
5 (a), we show the response function $A_{X}(\omega)$ of $\hat{X}_{1}$, which
peaks at $\omega=\pi/T$ when $\nu_{\pi}=-1$. These numerical results indicate
the presence of the DTC for the end spin observable. But in BDI class,
according to the bulk-edge correspondenceAsbóth  _et al._ (2016),
$n_{\rm{edge}}=|\nu_{\pi}|$ and $\nu_{\pi}>0$ ($\nu_{\pi}<0$) represents the
number of edge states occupying $\hat{\gamma}_{1}|0\rangle$
($\hat{\gamma}_{2}|0\rangle$). Therefore, $|W_{\eta}\rangle_{1}$ is finite
when $\nu_{\pi}>0$, while $|W_{\eta}\rangle_{1}$ is exactly $0$ when
$\nu_{\pi}<0$. As shown in Fig. 5 (b), $A(\omega)$ peaks at $\omega=\pi/T$
when $\nu_{\pi}>0$, which also indicats the presence of the DTC for BDI class.
Note that the DTC regime, where $A(\omega)$ peaks at $\omega=\pi/T$, doesn’t
exactly match with the topological region with $\nu_{\pi}=-1$ ($\nu_{\pi}>0$)
for D (BDI) class. We think this phenomenon is due to the finite size effect,
which causes the quasi-energy of edge mode to be not perfectly $\pi/T$. In
order to check the finite size effect, the numerical results for different
lengths of the chain are shown in Appendix B.

##  V Robustness of the DTC

In order to demonstrate the robustness of the DTC induced by topological
superconductors, we continue to investigate the frequency distributions of
response function for $\hat{X}_{1}$ after adding symmetry-preserving and
symmetry-breaking perturbations.

For D class or BDI class, the Hamiltonian obeys particle-hole symmetry, which
means the Hamiltonian (4) is invariant after particle-hole transformation
$\hat{c}^{\dagger}_{j}\rightarrow(-1)^{j}\hat{c}_{j}$,
$\hat{c}_{j}\rightarrow(-1)^{j}\hat{c}^{\dagger}_{j}$, and $i\rightarrow{-i}$.
In the following, we add an operator $\sum_{j}\delta_{j}\hat{Z}_{j}$ as the
particle-hole symmetry-preserving perturbation and $\sum\delta_{j}\hat{X}_{j}$
as the symmetry-breaking perturbation, with $\delta_{j}$ disordered
parameters. As shown in Fig. 6 (a) and (b), the frequency distributions of Eq.
(39), for D class and BDI class respectively, are almost unaffected by the
symmetry-preserving perturbation. As a comparison, the sub-harmonic response
is collapsed after adding symmetry-breaking perturbations for both classes.
These results suggest that sub-harmonic response is robust to symmetry-
preserving perturbations, which confirms the DTC is originated from
topologically non-trivial phases.

Figure 6: (a-d) Frequency distributions of the response function for
observable $\hat{X}_{1}$ and the color represents the amplitude of
$A_{X}(\omega)$. In (a) and (b), $\sum_{j}\delta_{j}\hat{Z}_{j}$ is added
which preserves the particle-hole symmetry, while in (c) and (d)
$\sum_{j}\delta_{j}\hat{X}_{j}$ is added which breaks the particle-hole
symmetry. Parameters of (a, c) are the same as the black solid line in Fig. 2,
while parameters of (b,d) are the same as the black solid line in Fig. 3.
$\delta_{j}$ is randomly distributed in $[-0.1h_{2}^{z},0.1h_{2}^{z}]$, which
is only added in the second interval during a period.

##  VI Summary

We have systematically studied the Floquet time crystal in a periodically
driven spin chain model by elaborating the topological classifications and
phase diagrams after mapping it to a Majorana chain. In this paper, we have
shown that the dynamics of the observable edge mode operator or end spin
operator indeed exhibit robust sub-harmonic oscillation, which is a typical
signature of the DTC. Besides, the rigidity and robustness of the DTC build on
topologically non-trivial phases. Furthermore, since the system is a general
spin chain model and can be strictly solved, it is potentially generalized to
other interacting or dissipative systems. Our results might be helpful in
deeply understanding the mechanism for other kinds of DTC. Besides, the model
is easily implemented and the topological DTC for the end spin are readily
realized in experiments.

Acknowledgements We would like to thank Hui Zhai, Wei Zheng, Chengshu Li,
Yanting Cheng, Haifeng Tang and Fan Yang for helpful discussions. This work
has been supported by the National Natural Science Foundation of China (Grant
No. 12204428).

##  Appendix A Proof of Eq. (8)

In this appendix, we prove that the effective Hamiltonian
$\hat{H}_{\rm{eff}}$, defined by

$\exp(-i\hat{H}_{{\rm
eff}}T)=\exp(-i\hat{H_{2}}t_{2})\exp(-i\hat{H_{1}}t_{1}),$ (42)

is quadratic with
$\hat{H}_{\rm{eff}}=\frac{1}{2}\hat{\Psi}^{\dagger}{H}_{\rm{eff}}\hat{\Psi}$,
where $H_{\rm{eff}}$ has such form

$\exp\left(-iH_{\mathrm{eff}}T\right)=\exp\left(-iH_{2}t_{2}\right)\exp\left(-iH_{1}t_{1}\right).$
(43)

First, we rewrite $\hat{H}_{\rm{eff}}$ by means of the
Baker–Campbell–Hausdorff formula as

$\displaystyle\hat{H}_{{\rm
eff}}=\frac{1}{T}\\{\hat{H}_{1}t_{1}+\hat{H}_{2}t_{2}-\frac{it_{1}t_{2}}{2}[\hat{H}_{1},\hat{H}_{2}]$
(44)
$\displaystyle-\frac{t_{1}^{2}t_{2}}{12}[\hat{H}_{1},[\hat{H}_{1},\hat{H}_{2}]]-\frac{t_{1}^{2}t_{2}}{12}[\hat{H}_{2},[\hat{H}_{2},\hat{H}_{1}]]\ldots\\}.$

Second, since the right-hand of Eq. (44) is made up of a series of
commutators, we choose $[\hat{H}_{1},\hat{H}_{2}]$ as a typical term to
analyse. Because both $\hat{H}_{1}$ and $\hat{H}_{2}$ have the quadratic form
as $\hat{H}_{m}=\frac{1}{2}\hat{\Psi}^{\dagger}H_{m}\hat{\Psi}$, we define a
mapping from a matrix to an operator as
$\hat{\Gamma}(H_{m})=\frac{1}{2}\sum_{i,j}(H_{m})_{i,j}\gamma_{i}\gamma_{j}$.

Third, we prove the relation
$\hat{\Gamma}([H_{1},H_{2}])=[\hat{\Gamma}(H_{1}),\hat{\Gamma}(H_{2})]$ in the
following process,

$\displaystyle[\hat{H_{1}},\hat{H}_{2}]$
$\displaystyle=\frac{1}{4}\sum_{ij}\sum_{u,v}(H_{1})_{i,j}(H_{2})_{u,v}(\gamma_{i}\gamma_{j}\gamma_{u}\gamma_{v}-\gamma_{u}\gamma_{v}\gamma_{i}\gamma_{j})$
(45)
$\displaystyle=\frac{1}{4}\sum_{i,j}\sum_{u,v}(H_{1})_{i,j}(H_{2})_{u,v}\gamma_{i}(\delta_{u,j}-\gamma_{u}\gamma_{j})\gamma_{v}-\frac{1}{4}\sum_{i,j}\sum_{u,v}(H_{1})_{i,j}(H_{2})_{u,v}\gamma_{u}(\delta_{v,i}-\gamma_{i}\gamma_{v})\gamma_{j}$
$\displaystyle=\frac{1}{4}\sum_{i,j}\sum_{u,v}(H_{1})_{i,j}(H_{2})_{u,v}\delta_{u,j}\gamma_{i}\gamma_{v}-\frac{1}{4}\sum_{i,j}\sum_{u,v}(H_{2})_{u,v}(H_{1})_{i,j}\delta_{v,i}\gamma_{u}\gamma_{j}$
$\displaystyle-\frac{1}{4}\sum_{i,j}\sum_{u,v}(H_{1})_{i,j}(H_{2})_{u,v}(\delta_{u,i}-\gamma_{u}\gamma_{i})\gamma_{j}\gamma_{v}+\frac{1}{4}\sum_{i,j}\sum_{u,v}(H_{1})_{i,j}(H_{2})_{u,v}\gamma_{u}\gamma_{i}(\delta_{v,j}-\gamma_{j}\gamma_{v})$
$\displaystyle=\frac{1}{2}(H_{1}H_{2}-H_{2}H_{1})_{i,j}\gamma_{i}\gamma_{j}$
$\displaystyle=\hat{\Gamma}([H_{1},H_{2}]),$

where we have used the property that matrices $H_{1}$ and $H_{2}$ are both
anti-symmetric. Using the relation of Eq. (45), we can prove
$[\hat{H}_{1},[\hat{H}_{1},\hat{H}_{2}]]=\hat{\Gamma}([{H}_{1},[{H}_{1},{H}_{2}]])$,
$[\hat{H}_{2},[\hat{H}_{2},\hat{H}_{1}]]=\hat{\Gamma}([{H}_{2},[{H}_{1},{H}_{2}]])$,
and so on. In this way, every term in Eq. (44) can be rewritten as the
quadratic form $\hat{\Gamma}$, so we have
$\hat{H}_{\rm{eff}}=\hat{\Gamma}(H_{\rm{eff}})$, where

$\displaystyle{H}_{{\rm
eff}}=\frac{1}{T}\\{{H}_{1}t_{1}+\hat{H}_{2}t_{2}-\frac{it_{1}t_{2}}{2}[{H}_{1},{H}_{2}]$
(46)
$\displaystyle-\frac{t_{1}^{2}t_{2}}{12}[{H}_{1},[{H}_{1},{H}_{2}]]-\frac{t_{1}^{2}t_{2}}{12}[{H}_{2},[{H}_{2},{H}_{1}]]\ldots\\}.$

Again, we use the Baker–Campbell–Hausdorff formula and $H_{\rm{eff}}$ is
obtained according to
$\exp(-iH_{\rm{eff}}T)=\exp(-iH_{1}t_{1})\exp(-iH_{2}t_{2})$.

##  Appendix B Finite size effect

In Fig. 5, the DTC regime, where $A_{X}(\omega)$ peaks at $\omega=\pi/T$, is
smaller than the region of topologically non-trivial phases. We think this
phenomenon is due to the finite size effect, which causes the quasi-energy of
anomalous edge state is not perfectly $\pi/T$. As an example, we numerically
calculate the DTC regimes for BDI class as the length of chain $N$ increases
to show the finite size effect. The numerical results are shown in Fig. 7,
from which we find the regime of parameter $\theta_{1}$ satisfying
$A_{X}(\omega=\pi/T)>0.05$ increases as the length of chain increases. It is
inferred that the regime of parameter $\theta_{1}$ could match with the region
of topologically non-trivial phases as $N\rightarrow\infty$.

Figure 7: (a-f) Frequency distributions of the response function for
observable $\hat{X}_{1}$ and the color represents the amplitude of
$A_{X}(\omega)$. The length of chian $N$ increases from $10$ to
$12,14,16,18,20$, corresponding to (a-f), respectively. Parameters are the
same as Fig. 5 except for $\theta_{1}\in[3\pi/16,\pi/2]$, which is the
topologically non-trivial regime with $\nu_{\pi}\neq 0$. (g) The regime of
parameter $\theta_{1}$ satisfying $A_{X}(\omega=\pi/T)>0.05$ varies with the
length of chain $N$.

## References

  * Zenesini  _et al._ (2009) A. Zenesini, H. Lignier, D. Ciampini, O. Morsch, and E. Arimondo, Phys. Rev. Lett. 102, 100403 (2009).
  * Xu  _et al._ (2021) P. Xu, T.-S. Deng, W. Zheng, and H. Zhai, Phys. Rev. A 103, L061302 (2021).
  * Clark  _et al._ (2018) L. W. Clark, B. M. Anderson, L. Feng, A. Gaj, K. Levin, and C. Chin, Phys. Rev. Lett. 121, 030402 (2018).
  * Chen  _et al._ (2020a) Y.-Y. Chen, P. Zhang, W. Zheng, Z. Wu, and H. Zhai, Phys. Rev. A 102, 011301 (2020a).
  * Chen  _et al._ (2020b) Y.-Y. Chen, P. Zhang, W. Zheng, Z. Wu, and H. Zhai, Phys. Rev. A 102, 011301 (2020b).
  * Lv  _et al._ (2020) C. Lv, R. Zhang, and Q. Zhou, Phys. Rev. Lett. 125, 253002 (2020).
  * Jotzu  _et al._ (2014) G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, Nature 515, 237 (2014).
  * Aidelsburger  _et al._ (2015) M. Aidelsburger, M. Lohse, C. Schweizer, M. Atala, J. T. Barreiro, S. Nascimbène, N. Cooper, I. Bloch, and N. Goldman, Nat. Phys. 11, 162 (2015).
  * Zhang  _et al._ (2017) J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, _et al._ , Nature 543, 217 (2017).
  * Choi  _et al._ (2017) S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, _et al._ , Nature 543, 221 (2017).
  * Rovny  _et al._ (2018a) J. Rovny, R. L. Blum, and S. E. Barrett, Phys. Rev. Lett. 120, 180603 (2018a).
  * Rovny  _et al._ (2018b) J. Rovny, R. L. Blum, and S. E. Barrett, Phys. Rev. Lett. 120, 180603 (2018b).
  * Pal  _et al._ (2018) S. Pal, N. Nishad, T. S. Mahesh, and G. J. Sreejith, Phys. Rev. Lett. 120, 180602 (2018).
  * Keßler  _et al._ (2021) H. Keßler, P. Kongkhambut, C. Georges, L. Mathey, J. G. Cosme, and A. Hemmerich, Phys. Rev. Lett. 127, 043602 (2021).
  * Randall  _et al._ (2021) J. Randall, C. Bradley, F. van der Gronden, A. Galicia, M. Abobeih, M. Markham, D. Twitchen, F. Machado, N. Yao, and T. H. Taminiau, Science 374, 1474 (2021).
  * Kyprianidis  _et al._ (2021) A. Kyprianidis, F. Machado, W. Morong, P. Becker, K. Collins, D. Else, L. Feng, P. Hess, C. Nayak, G. Pagano, N. Yao, and C. Monroe, Science 372, 1192 (2021).
  * Bluvstein  _et al._ (2021) D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. Wang, A. Michailidis, N. Maskara, W. W. Ho, _et al._ , Science 371, 1355 (2021).
  * Mi  _et al._ (2022) X. Mi, M. Ippoliti, C. Quintana, A. Greene, Z. Chen, J. Gross, F. Arute, K. Arya, J. Atalaya, R. Babbush, J. Bardin, J. Basso, A. Bengtsson, A. Bilmes, A. Bourassa, L. Brill, M. Broughton, B. Buckley, D. Buell, and P. Roushan, Nature 601, 1 (2022).
  * Zhang  _et al._ (2022) X. Zhang, W. Jiang, J. Deng, K. Wang, J. Chen, P. Zhang, W. Ren, H. Dong, S. Xu, Y. Gao, _et al._ , Nature 607, 468 (2022).
  * Frey and Rachel (2022) P. Frey and S. Rachel, Sci. Adv. 8, eabm7652 (2022).
  * Khemani  _et al._ (2016a) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016a).
  * Else  _et al._ (2016a) D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. 117, 090402 (2016a).
  * Yao  _et al._ (2017) N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Phys. Rev. Lett. 118, 030401 (2017).
  * Huang  _et al._ (2018) B. Huang, Y.-H. Wu, and W. V. Liu, Phys. Rev. Lett. 120, 110603 (2018).
  * Gong  _et al._ (2018) Z. Gong, R. Hamazaki, and M. Ueda, Phys. Rev. Lett. 120, 040404 (2018).
  * Zhu  _et al._ (2019) B. Zhu, J. Marino, N. Y. Yao, M. D. Lukin, and E. A. Demler, New J. Phys. 21, 073028 (2019).
  * Chinzei and Ikeda (2020) K. Chinzei and T. N. Ikeda, Phys. Rev. Lett. 125, 060601 (2020).
  * Huang  _et al._ (2022) B. Huang, T.-H. Leung, D. M. Stamper-Kurn, and W. V. Liu, Phys. Rev. Lett. 129, 133001 (2022).
  * Wilczek (2012) F. Wilczek, Phys. Rev. Lett. 109, 160401 (2012).
  * Bruno (2013) P. Bruno, Phys. Rev. Lett. 111, 070402 (2013).
  * Watanabe and Oshikawa (2015) H. Watanabe and M. Oshikawa, Phys. Rev. Lett. 114, 251603 (2015).
  * Khemani  _et al._ (2016b) V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. 116, 250401 (2016b).
  * Else  _et al._ (2016b) D. V. Else, B. Bauer, and C. Nayak, Phys. Rev. Lett. 117, 090402 (2016b).
  * Rudner  _et al._ (2013) M. S. Rudner, N. H. Lindner, E. Berg, and M. Levin, Phys. Rev. X 3, 031005 (2013).
  * Titum  _et al._ (2016) P. Titum, E. Berg, M. S. Rudner, G. Refael, and N. H. Lindner, Phys. Rev. X 6, 021013 (2016).
  * Zhou and Gong (2018) L. Zhou and J. Gong, Phys. Rev. B 97, 245430 (2018).
  * Bomantara and Gong (2018) R. W. Bomantara and J. Gong, Phys. Rev. Lett. 120, 230405 (2018).
  * Pan and Wang (2020) Y. Pan and B. Wang, Phys. Rev. Research 2, 043239 (2020).
  * Chew  _et al._ (2020) A. Chew, D. F. Mross, and J. Alicea, Phys. Rev. Lett. 124, 096802 (2020).
  * Shtanko and Movassagh (2020) O. Shtanko and R. Movassagh, Phys. Rev. Lett. 125, 086804 (2020).
  * Bomantara  _et al._ (2021) R. W. Bomantara, S. Mu, and J. Gong, Phys. Rev. B 103, 235404 (2021).
  * Wang  _et al._ (2022) B. Wang, J. Quan, J. Han, X. Shen, H. Wu, and Y. Pan, Laser Photonics Rev. , 2100469 (2022).
  * (43) The BdG constraint requires $C^{-1}HC=-H$ with $C^{2}=1$, which means only Class D, BDI and DIII are allowed in 1-D system. However, for Class DIII it also requires $T^{-1}HT=-H$ with $T^{2}=-1$, which leads to a trivial result for a two-band model.
  * Chiu  _et al._ (2016) C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * Jiang  _et al._ (2011) L. Jiang, T. Kitagawa, J. Alicea, A. R. Akhmerov, D. Pekker, G. Refael, J. I. Cirac, E. Demler, M. D. Lukin, and P. Zoller, Phys. Rev. Lett. 106, 220402 (2011).
  * Asbóth and Obuse (2013) J. K. Asbóth and H. Obuse, Phys. Rev. B 88, 121406 (2013).
  * Asbóth  _et al._ (2016) J. K. Asbóth, L. Oroszlány, and A. Pályi,  _A short course on topological insulators_ , Vol. 919 (Springer, 2016).