# Emergence of Competing Orders and Possible Quantum Spin Liquid in SU(N)
Fermions

Xue-Jia Yu  International Center for Quantum Materials, School of Physics,
Peking University, Beijing 100871, China  Shao-Hang Shi  Beijing National
Laboratory for Condensed Matter Physics and Institute of Physics, Chinese
Academy of Sciences, Beijing 100190, China  University of Chinese Academy of
Sciences, Beijing 100049, China  Limei Xu  International Center for Quantum
Materials, School of Physics, Peking University, Beijing 100871, China
Collaborative Innovation Center of Quantum Matter, Beijing, China
Interdisciplinary Institute of Light-Element Quantum Materials and Research
Center for Light-Element Advanced Materials, Peking University, Beijing, China
Zi-Xiang Li<EMAIL_ADDRESS>Beijing National Laboratory for Condensed
Matter Physics and Institute of Physics, Chinese Academy of Sciences, Beijing
100190, China  University of Chinese Academy of Sciences, Beijing 100049,
China

###### Abstract

In the past decades, tremendous efforts have been made towards understanding
the exotic physics emerging from competition between various ordering
tendencies in strongly correlated systems. Employing state-of-the-art quantum
Monte-Carlo simulation, we investigate an interacting SU($N$) fermionic model
with varying interaction strength and value of $N$, and unveil the ground-
state phase diagram of the model exhibiting a plethora of exotic phases. For
small value of $N$, namely $N=2,3$, the ground state is antiferromagnetic
(AFM) phase, whereas in the large-$N$ limit, valence bound solid (VBS) order
is dominant. For the intermediate value of $N$ such as $N=4$, remarkably, our
study reveals distinct VBS orders appear in the weak and strong coupling
regimes. More fantastically, the competition between staggered and columnar
VBS ordering tendencies gives rise to a Mott insulating phase without
spontaneously symmetry breaking (SSB), existing in a large interacting
parameter regime, which is consistent with a gapped quantum spin liquid. Our
study not only provides a platform to investigate the fundamental physics of
quantum many-body systems, but also offers a novel route towards searching for
exotic states of matter such as quantum spin liquid in realistic quantum
materials.

Figure 1: Schematic ground-state phase diagram of SU($N$) fermions with SSH
interaction on square lattice at half filling: for $N=2,3$, the ground state
is AFM ordered phase. In the large $N$ limit, the dominant VBS ordering is
sVBS with momentum $(\pi,\pi)$. For the intermediate $N$, such as $N=4$, the
insulating phases in weakly and strongly interacting regime possess VBS orders
with distinct patterns of SSB. The sVBS order carrying momentum $(\pi,\pi)$ is
dominate in the weakly interacting regime owing to the presence of Fermi-
surface nesting, while in strongly coupling regime, cVBS order with momentum
$(\pi,0)/(0,\pi)$ appears. Most appealingly, in a large regime between two
different VBS phases ($1.5<J<3.5$), the ground state is a Mott insulator
without SSB consistent with gapped quantum spin liquid.

Introduction.- Understanding competing orders arising from strong electronic
interactions plays vital roles in fathoming fundamental theory of quantum
many-body physics [1, 2], as well as many appealing phenomena in quantum
materials including high $T_{c}$ superconductivity [3, 4, 5]. More
intriguingly, the competitions between various ordering tendencies offer a
route towards yielding exotic states of matter. A typical example is quantum
spin liquid (QSL) [6, 7, 8, 9, 10, 11], a Mott insulating phase featuring
deconfined fractionalization, whereas without symmetry spontaneously breaking
even at zero temperature. Despite vigorous studies in past decades [12, 13,
14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32,
33, 34], unambiguous demonstration of realizing quantum spin liquid in
realistic microscopic models by unbiased theoretical approach remains elusive.

On a different front, another promising mechanism to trigger exotic state of
matter is assigning more degrees of freedom to the microscopic models and
extending symmetries. By extending spin symmetry group SU($2$) in electronic
spin system to SU($N$), it is theoretically predicted that quantum fluctuation
of spin is enhanced with increasing $N$[35, 36], rendering the possibility of
exotic quantum paramagnetic phases such as valence bond solid (VBS)[37, 38,
39, 40] and quantum spin liquid[41, 42, 27]. Fruitful fascinating physics has
been revealed in the interacting SU($N$) systems by virtue of the cooperative
effects of enlarged symmetry and strong correlation [43, 44, 45, 46, 47, 48,
49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,
68, 69, 70, 71, 72, 73, 74]. Benefiting from celebrated progresses in cold-
atom experiments, various SU($N$) symmetric systems are realized in optical
lattice [75, 76, 77, 78, 79, 80, 81]. More recently, enlarged SU($N$) symmetry
is proposed to emerge in the twisted bilayer graphene system [82, 83, 84] as a
low-energy effective description combining electronic spin and valley degrees
of freedom[85, 86, 87, 88]. Hence, designing concrete model with extended
symmetry and investigating the effect of strong correlation in a theoretically
controlled approach is vastly desired.

In this letter, we construct a SU($N$) interacting fermionic model on square
lattice. Remarkably, the model is sign problem free at half filling for any
$N$[89, 90, 91, 92], thus amendable to approximation-free Quantum Monte Carlo
(QMC) simulation [93, 94, 95] at low temperature (or zero temperature) with
large system sizes. We perform large-scale projector QMC simulation and
investigate numerically accurate ground-state properties of the model with
varying interaction strength and value of $N$. At small value of $N$, the
ground state is AFM spin ordered state, like SU(2) Heisenberg model on square
lattice. With increasing value of $N$, quantum fluctuation of spin destroys
spin long-range order, rendering VBS ordering breaking translational symmetry,
whereas preserving spin rotational symmetry. The most intriguing phenomena
occur in the intermediate value of $N$, wherein the insulating phases in
weakly coupling regime and strongly coupling regime feature distinct patterns
of VBS. The competition between two ordering tendencies gives rise to an
intermediate insulating phase without any SSB even at zero temperature. Such
symmetric Mott insulator is a possible gapped quantum spin liquid phase
featuring topological order [96, 97, 98, 99].

Model and Method.- We consider the following interacting microscopic model of
SU($N$) fermions on the square lattice:[100](SSH interaction model):

$H=-t\sum_{\langle
ij\rangle,\alpha}(c^{\dagger}_{i\alpha}c_{j\alpha}+H.c.)-\frac{J}{2N}\sum_{\langle
ij\rangle}(\sum_{\alpha}c^{\dagger}_{i\alpha}c_{j\alpha}+H.c.)^{2}$ (1)

where $\langle ij\rangle$ refers to the bond between nearest-neighbor (NN)
sites $i$ and $j$, $c^{\dagger}_{i\alpha}$ creates a fermion on-site $i$ with
fermion flavor label $\alpha=1,...N$. Here $t$ is the fermion hopping
amplitude and is set $t=1$ as an energy unit. $J$ represents the strength of
the SSH interaction, which can be induced by SSH electron-phonon couplings in
the fast phonon limit[101, 102]. Remarkably, model in Eq. (1) is free from the
notorious sign problem for any $N$ at half filling and even $N$ at generic
filling, and thus a promising platform to investigate the exotic physics
emerging in the SU($N$) fermions. In this work, we perform un-biased, large-
scale projector QMC to study the ground-state properties of the model. The
details of projector QMC are introduced in the Supplementary Materials (SM).
We focus on the model at half filling with perfectly nested Fermi surface, and
explore the insulating phase with competing ordering arising from the
interplay between strong correlation and extended symmetry of multi-component
fermions.

Notice that for $N=1$, the interaction in Eq. (1) is reduced to the nearest-
neighbor density interaction, which favors charge-density-wave order on
bipartite lattice, as intensively studied previously[103, 104, 105, 106]. For
$N=2$, the model is SSH electron-phonon coupling model at anti-adiabatic
limit, the ground-state of which is revealed to be AFM ordered state in recent
works [107, 108, 109, 110]. In the present work, we focus on understanding the
nature of the multi-component fermions with enlarged symmetry, namely $N\geq
3$, in the presence of strong correlation.

Quantum phase diagram.- Before presenting the details of QMC results, we
summarize our main findings and the ground-state phase diagram of the model in
Eq. (1). The schematic phase diagram with varying fermion-flavor number $N$ is
shown in Fig. 1. For small $N$, more explicitly $N=2,3$, the ground state is
an insulating phase with AFM spin long-range order. As discussed in previous
literature, with increasing fermion-flavor number $N$, the quantum fluctuation
of spin is enhanced, resulting in the emergence of valence bond ordering. The
emergence of VBS order for large-$N$ case is revealed by our systematic
numerical calculation. In the large-$N$ limit, where the saddle-point
approximation can render the accurate results, we perform mean-field
calculation and unambiguously show that the ground state is VBS ordered phase
with momentum $(\pi,\pi)$, namely a staggered VBS (sVBS) phase. In the case of
a large but finite $N$, we perform QMC simulation to decipher the rich ground-
state phase diagram. For example, for $N=4$, the AFM spin order is suppressed
and VBS order becomes dominant. More strikingly, our state-of-the-art QMC
simulation shows that the patterns of VBS orders in the weak and strong
coupling regimes are completely different. In the weakly interacting regime,
the VBS order in the ground state carries momentum $(\pi,\pi)$ as a
consequence of Fermi surface nesting. In the strongly coupling regime, the
dominant pattern of VBS ordering is the colomnar VBS (cVBS) with momentum
$(\pi,0)/(0,\pi)$. Appealingly, in a large regime with intermediate
interacting strength ($1.5<J<3.5$), a phase without SSB arises from the
competition between two VBS ordering tendencies. This phase without any SSB in
the intermediate strength of $J$ possesses finite single-particle and spin
gaps. Therefore, it is a Mott insulating phase, and possibly, a spin liquid
phase with exotic topological order.

Figure 2: Numerical evidence of AFM order in $N=3$: (a) The color plot of
momentum distribution of AFM structure factor for $J=2,L=16$, clearly showing
that peak occurring at momentum $(\pi,\pi)$. (b) Finite-size scaling of spin
and VBS structure factor at peaked momentum $(\pi,\pi)$. The second-order
polynomial is used to fit the structure factor as a function of $1/L$, and the
intercept of extrapolation indicates the structure factor at thermodynamic
limit. (c) RG-invariant AFM correlation-length ratio as a function of
interacting strength $J$ for system size $L=10,12,14,16$. Increasing
correlation ratio with system size is a strong evidence of long-range order.
(d) The extrapolation of spin gap as a function of $1/L$ for $J=2,3,4$.
Second-order pronominal is used to fit the structure factors. The gapless spin
excitation is consistent with AFM long-range ordering.

Staggered VBS order in the large $N$ limit.- In the large-$N$ limit for Eq.
(1), the mean-field saddle-point approximation is capable of capturing the
exact ground-state properties. Hence, before presenting QMC results on
finite-$N$ system, we perform self-consistent mean-field calculation to
understand the ground-state properties of Eq. (1) in the large-$N$ limit. We
decouple the interaction in Eq. (1) in bond channel:
$\varphi_{x/y,i}=\langle(c^{\dagger}_{i}c_{i+\hat{x}/\hat{y}}+{\rm{h.c}})\rangle$
and solve the mean-field equation self-consistently. Our calculation indicates
sVBS with momentum $(\pi,\pi)$ and cVBS with momentum $(\pi,0)/(0,\pi)$ are
two dominant ordered states with close energies. We perform systematic self-
consistent mean-field calculation to compare energies of the two states,
unequivocally revealing that sVBS ordered state is energetically favored in
the whole interaction parameter regime under consideration. The details and
results of mean-field calculation are included in the section V of SM. To
conclude, in the large-$N$ limit, the ground state of Eq. (1) is sVBS ordered
phase with momentum $(\pi,\pi)$.

AFM order for small $N$.- At small or intermediate value of $N$, where most
physical realizations reside, quantum fluctuation beyond saddle-point
approximation becomes pronounced. Consequently, an approximation-free approach
incorporating the effect of quantum fluctuation is strongly desired to
investigate the possible SSB at small or intermediate value of $N$. We employ
numerically exact QMC simulation to explore the accurate ground-state
properties of Eq. (1). To demystify various possible symmetry-breaking orders,
we compute the static structure factor
$S(\vec{q},L)=\frac{1}{L^{4}}\sum_{i,j}e^{i\vec{q}\cdot(\vec{r}_{i}-\vec{r}_{j})}\langle\hat{O}_{i}\hat{O}_{j}\rangle$
of the corresponding order $\hat{O}$ and $L$ refers to the linear system size
of square lattice. The order parameter characterizing SSB is given by the
peaked-momentum structure factor at thermodynamic limit:
$\Delta^{2}=\lim_{L\rightarrow\infty}S(\vec{Q}_{\rm peak},L)$. In our study,
we focus on VBS and spin orderings, the definition of which are given in the
SM I. Firstly, we present the momentum distribution of spin structure factor
in Fig. 2(a), which clearly displays a sharp peak at momentum $(\pi,\pi)$,
implying the existence of AFM ordering at $N=3$. The scaling analysis of AFM
structure factor with $1/L$ unambiguously reveals the existence of AFM long-
range order, as shown in Fig. 2(b), while the VBS structure factor vanishes at
thermodynamic limit. To further confirm the AFM long-range ordered ground
state at $N=3$, we perform a more sophisticated finite-size scaling procedure
by calculating RG-invariant correlation-length ratio,
$R(L)_{{\rm{AFM/VBS}}}=\frac{S_{{\rm{AFM/VBS}}}(\vec{Q},L)}{S_{{\rm{AFM/VBS}}}(\vec{Q}-\delta\vec{q},L)}-1$,
where $\vec{Q}$ labels the momentum at which the structure factor is maximum,
$\delta\vec{q}=(\frac{2\pi}{L},\frac{2\pi}{L})$ is a minimal momentum shift
from $\vec{Q}$. As shown in Fig. 2(c), the correlation-length ratio of AFM
order increases with system size, corroborating the existence of AFM long-
range order. Conversely, the system-size dependence of VBS correlation-length
ratio, as shown in SM III, is a strong manifestation of short-range property
of VBS ordering at $N=3$.

In addition to the static structure factor, we investigate the spectral
properties of Eq. (1) by virtue of imaginary-time Green’s function. At $N=3$,
we extract spectral gap from imaginary-time Green’s function, with the details
illustrated in SM II. The results of single-particle gaps at $N=3$ with
varying $J$ are shown in SM III, explicitly demonstrating the existence of
finite single-particle gap at thermodynamic limit. Furthermore, in the AFM
ordered phase, the spin gap is expected to vanish at momentum $(\pi,\pi)$ as a
consequence of Goldstone modes. Fig. 2(d) depicts the scaling analysis of spin
gap as a function of $1/L$, explicitly indicating the absence of spin gap, and
hence corroborating the existence of AFM long-range order in the ground state
at $N=3$.

Figure 3: The results of correlation-length ratio and structure factors for
$N=4$: (a) RG-invariant correlation length ratio of sVBS order. The results
increase with system size in the weakly interacting regime $J<1.5$, indicating
sVBS ordered phase. (b) Correlation-length ratio of cVBS order. The results
increase with system size in the strongly interacting regime $J>3.5$,
indicating a cVBS ordered phase. (c) Finite-size scaling of sVBS structure
factor in the regime of intermediate interaction strength $1.5\leq J\leq 3.0$.
(d) Finite-size scaling of sVBS structure factor in the regime of intermediate
interaction strength $1.5\leq J\leq 3.0$. Both sVBS and cVBS structure factors
vanish in thermodynamic limit.

Intermediate Mott insulator without SSB.- In this section, we perform
projector QMC to investigate ground-state properties of the model Eq. (1) in
the intermediate values of $N$, wherein the competition between various
ordering tendencies is expected to result in exotic phases. We focus our study
on $N=4$. For $N=4$, the results clearly show that AFM order is destroyed by
the quantum fluctuation of spin by evaluating structure factor and RG-
invariant correlation ratio, as shown in section IV of SM. The dominant
instability for $N=4$ is VBS ordering. We present the results of correlation-
length ratio for sVBS and cVBS orders in Fig. 3(a) and (b), respectively. At
weak coupling regime, the dominant VBS ordering is sVBS with momentum
$(\pi,\pi)$ owing to the Fermi surface nesting on square lattice at half
filling. Remarkably, in the presence of strong interaction, the ground state
displays cVBS order with momentum $(\pi,0)/(0,\pi)$, distinct from the VBS
pattern at weak coupling regime. Different from sVBS at weak coupling regime
arising from Fermi surface instablity, the cVBS order is a consequence of
strong electronic interaction. The results of correlation-length ratio
unambiguously establish the existence of sVBS and cVBS long-range order in the
region $J<1.5$ and $J>3.5$, respectively. Hence, for the intermediate fermion-
component number $N=4$, our simulation unequivocally reveals two different VBS
ordered phase emerging in the weak and strong coupling regimes, with distinct
patterns and underlying mechanisms.

Figure 4: Momentum distribution of VBS structure factors for $N=4$: (a)
$J=0.5$. (b) $J=2.0$. (c) $J=3.0$. (d) $J=5.0$. The linear system size is
fixed as $L=20$.

The competition between various ordered phases with incompatible broken
symmetries offers a promising route towards yielding exotic QSL phase. In Eq.
(1) for $N=4$, the results of correlation-length ratio reveal an intermediate
phase in the region $1.5<J<3.5$, in which both sVBS and cVBS orders are short-
range. The finite-size scaling of VBS structure factors with momentum
$(\pi,\pi)$ and $(\pi,0)$, as depicted in Fig. 3(c) and (d), also confirm the
absence of VBS long-range ordering. We exclude the possibility of other SSB
orderings occurring in this intermediate phase by means of finite-scaling
analysis of structure factor, as shown in the section IV of SM, corroborating
the existence of an intermediate phase without any SSB. The occurrence of
intermediate phase in the absence of VBS long-range order is also witnessed by
the momentum distribution of VBS structure factors, as shown in Fig. 4. In the
weak coupling ($J=0.5$) and strong coupling ($J=5.0$) regimes, as expected,
the VBS structure factors are sharply peaked at momentum $(\pi,\pi)$ and
$(\pi,0)$, respectively. In the intermediate phase ($J=2.0,3.0$), the peaks of
VBS structure factors are obviously broadened, implying the suppression of
long-range ordering. To further scrutinize the nature of intermediate
symmetric phase, we present the results of correlation length normalized by
the linear-system size $L$,estimated from the correlation-length ratio
$\xi=\frac{1}{\delta q}\sqrt{\frac{S(Q)}{S(Q+\delta q)}-1}$, in Fig. 5(a). In
the intermediate regime $1.5\textless J\textless 3.5$, the correlation length
is much smaller compared with linear system size $L$, confirming that VBS
order is short-range and our computational system size is sufficiently large
to capture the physical properties at thermodynamic limit. Taken together, our
simulation provides compelling evidences supporting the emergence of an
intermediate phase without SSB for $N=4$ in model Eq. (1).

Figure 5: (a) Correlation length normalized by linear system size $L$
extracted from RG-invariant correlation ratio for $N=4$. In the intermediate
symmetric Mott insulating phase, the correlation length of VBS order is much
smaller than $L$. (b) Finite-size scaling of spin gap in the intermediate
symmetric Mott insulating phase for $J=2.0,2.5,3.0$ for $N=4$. (c) Finite-size
scaling of single-particle gap in the intermediate symmetric Mott insulating
phase for $J=2.0,2.5,3.0$ for $N=4$.

To further investigate the nature of intermediate symmetric phase, we
calculate the single-particle and spin gaps. The finite-size scaling results
of spin and single-particle gaps for various values of $J$ are shown in Fig.
5(b) and (c), respectively. The extrapolated values unambiguously reveal
finite single-particle and spin gaps, demonstrating the nature of a gapped
Mott insulating phase, and possibly, a gapped quantum spin liquid featuring
topological order. Moreover, the intermediate gapped quantum spin liquid phase
is also consistent with the Lieb-Schultz-Mattis theorem [111, 112, 113, 114].
In our simulation, the average filling number is two electrons per unit cell,
hence the possibility of gapped SU(4) symmetric phase without ground-state
degeneracy is excluded owing to the LSM theorem for SU($N$) symmetric
system[114]. In conclusion, our systematic simulations establish a possible
gapped quantum spin liquid phase existing in the intermediate interacting
regime for $N=4$.

Discussions and concluding remarks: The SU($N$) interacting fermionic model
offers a promising platform to explore the intriguing physics arising from the
interplay between strong correlation and enlarged symmetry of multi-component
fermions. Here, we perform approximation-free numerical simulation to decipher
the ground-state properties of SU($N$) fermionic model with SSH interaction on
square lattice. Various competing orderings emerge at varying interaction
strength and the number of fermion-component $N$. At small value of $N$,
namely $N=2,3$, the ground state is AFM spin ordered state. However, in the
large-$N$ limit, the dominant order is sVBS carrying momentum $(\pi,\pi)$
while preserving spin rotational symmetry. The results at intermediate value
of $N$ are particularly intriguing. At $N=4$, in weak coupling regime the
dominant VBS ordering is sVBS with momentum $(\pi,\pi)$, persisting up to
$J=1.5$, whereas the strong interaction between electrons destroy sVBS
ordering and establish the cVBS ordered state carrying momentum
$(\pi,0)/(0,\pi)$ in region $J>3.5$. More remarkably, in a large intermediate
interacting regime ($1.5<J<3.5$), an exotic Mott insulating phase without SSB
develops before the cVBS order is established. The intermediate symmetric Mott
insulating phase possesses finite single-particle and spin gaps. All these
features are consistent with a gapped quantum spin liquid phase.

Notably, although we focus on the simulation on model in Eq. (1) at half
filling, the model is sign-problem-free at generic filling for even $N$. The
possible superconducting pairing arising from doping various symmetry
spontaneously breaking phases or quantum spin liquid is of particular
interest. Furthermore, Eq. (1) is sign-problem-free at half filling for even
$N$ in the presence of repulsive Hubbard interaction, making it feasible to
investigate the effect of Hubbard interactions on the various VBS and QSL
order unveiled in our simulation. Another interesting direction to explore in
future is the nature of quantum phase transition between quantum spin liquid
and different VBS ordered phases. Therefore, our simulation paves a promising
way to constructing concrete theoretically tractable model to investigate
these crucial ingredients of quantum many-body physics in a unified framework.

Acknowledgement.- We thanks Hong Yao, Yi Zhou, and Linhao Li for very helpful
discussions. We thank the computational resources provided by the TianHe-1A
supercomputer, the High Performance Computing Platform of Peking University,
China. The QMC simulations are partly carried out with ALF Library[115].
X.-J.Y. and L.X. are supported by the National Natural Science Foundation of
China under Grant No.11935002, and the National 973 project under Grant No.
2021YF1400501. Z.X.L acknowledges support from the start-up grant of IOP-CAS.

## References

  * Sachdev [2011] S. Sachdev,  _Quantum Phase Transitions_, 2nd ed. (Cambridge University Press, 2011).
  * Fradkin [2013] E. Fradkin,  _Field theories of condensed matter physics_ (Cambridge University Press, 2013).
  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006).
  * Fradkin  _et al._ [2015] E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Rev. Mod. Phys. 87, 457 (2015).
  * Sachdev [2003] S. Sachdev, Rev. Mod. Phys. 75, 913 (2003).
  * Balents [2010] L. Balents, Nature 464, 199 (2010).
  * Lee [2007] P. A. Lee, Reports on Progress in Physics 71, 012501 (2007).
  * Zhou  _et al._ [2017] Y. Zhou, K. Kanoda, and T.-K. Ng, Rev. Mod. Phys. 89, 025003 (2017).
  * Savary and Balents [2016] L. Savary and L. Balents, Reports on Progress in Physics 80, 016502 (2016).
  * Norman [2016] M. R. Norman, Rev. Mod. Phys. 88, 041002 (2016).
  * Broholm  _et al._ [2020] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Science 367, eaay0668 (2020), https://www.science.org/doi/pdf/10.1126/science.aay0668 .
  * Hermele  _et al._ [2004] M. Hermele, T. Senthil, M. P. A. Fisher, P. A. Lee, N. Nagaosa, and X.-G. Wen, Phys. Rev. B 70, 214437 (2004).
  * Hermele  _et al._ [2005] M. Hermele, T. Senthil, and M. P. A. Fisher, Phys. Rev. B 72, 104404 (2005).
  * Yan  _et al._ [2011] S. Yan, D. A. Huse, and S. R. White, Science 332, 1173 (2011).
  * Jiang  _et al._ [2012] H.-C. Jiang, H. Yao, and L. Balents, Phys. Rev. B 86, 024424 (2012).
  * Han  _et al._ [2012] T.-H. Han, J. S. Helton, S. Chu, D. G. Nocera, J. A. Rodriguez-Rivera, C. Broholm, and Y. S. Lee, Nature 492, 406 (2012).
  * Pratt  _et al._ [2011] F. L. Pratt, P. J. Baker, S. J. Blundell, T. Lancaster, S. Ohira-Kawamura, C. Baines, Y. Shimizu, K. Kanoda, I. Watanabe, and G. Saito, Nature 471, 612 (2011).
  * Meng  _et al._ [2010] Z. Meng, T. Lang, S. Wessel, F. Assaad, and A. Muramatsu, Nature 464, 847 (2010).
  * Sorella  _et al._ [2012] S. Sorella, Y. Otsuka, and S. Yunoki, Scientific reports 2, 1 (2012).
  * Liao  _et al._ [2017] H. J. Liao, Z. Y. Xie, J. Chen, Z. Y. Liu, H. D. Xie, R. Z. Huang, B. Normand, and T. Xiang, Phys. Rev. Lett. 118, 137202 (2017).
  * Yao and Kivelson [2007] H. Yao and S. A. Kivelson, Phys. Rev. Lett. 99, 247203 (2007).
  * Yao and Lee [2011] H. Yao and D.-H. Lee, Phys. Rev. Lett. 107, 087205 (2011).
  * Jiang  _et al._ [2019] Y.-F. Jiang, T. P. Devereaux, and H.-C. Jiang, Phys. Rev. B 100, 165123 (2019).
  * Gong  _et al._ [2015] S.-S. Gong, W. Zhu, L. Balents, and D. N. Sheng, Phys. Rev. B 91, 075112 (2015).
  * Zhou  _et al._ [2008] Y. Zhou, P. A. Lee, T.-K. Ng, and F.-C. Zhang, Phys. Rev. Lett. 101, 197201 (2008).
  * He  _et al._ [2017] Y.-C. He, M. P. Zaletel, M. Oshikawa, and F. Pollmann, Phys. Rev. X 7, 031020 (2017).
  * Jin  _et al._ [2022] H.-K. Jin, R.-Y. Sun, H.-H. Tu, and Y. Zhou, Science Bulletin 67, 918 (2022).
  * Wang  _et al._ [2021] Y.-C. Wang, Z. Yan, C. Wang, Y. Qi, and Z. Y. Meng, Phys. Rev. B 103, 014408 (2021).
  * Yan  _et al._ [2021] Z. Yan, Y.-C. Wang, N. Ma, Y. Qi, and Z. Y. Meng, npj Quantum Materials 6, 1 (2021).
  * Wan and Tchernyshyov [2013] Y. Wan and O. Tchernyshyov, Phys. Rev. B 87, 104408 (2013).
  * Song  _et al._ [2019] X.-Y. Song, C. Wang, A. Vishwanath, and Y.-C. He, Nature Communications 10, 4254 (2019).
  * Wang  _et al._ [2018] Y.-C. Wang, X.-F. Zhang, F. Pollmann, M. Cheng, and Z. Y. Meng, Phys. Rev. Lett. 121, 057202 (2018).
  * Szasz  _et al._ [2020] A. Szasz, J. Motruk, M. P. Zaletel, and J. E. Moore, Phys. Rev. X 10, 021042 (2020).
  * Chen  _et al._ [2022] B.-B. Chen, Z. Chen, S.-S. Gong, D. N. Sheng, W. Li, and A. Weichselbaum, Phys. Rev. B 106, 094420 (2022).
  * Read and Sachdev [1989a] N. Read and S. Sachdev, Nuclear Physics B 316, 609 (1989a).
  * Read and Sachdev [1989b] N. Read and S. Sachdev, Phys. Rev. Lett. 62, 1694 (1989b).
  * Lang  _et al._ [2013] T. C. Lang, Z. Y. Meng, A. Muramatsu, S. Wessel, and F. F. Assaad, Phys. Rev. Lett. 111, 066401 (2013).
  * Assaad [2005] F. F. Assaad, Phys. Rev. B 71, 075103 (2005).
  * Kaul and Sandvik [2012] R. K. Kaul and A. W. Sandvik, Phys. Rev. Lett. 108, 137201 (2012).
  * Zhou  _et al._ [2016] Z. Zhou, D. Wang, Z. Y. Meng, Y. Wang, and C. Wu, Phys. Rev. B 93, 245157 (2016).
  * Corboz  _et al._ [2012] P. Corboz, M. Lajkó, A. M. Läuchli, K. Penc, and F. Mila, Phys. Rev. X 2, 041013 (2012).
  * Koga  _et al._ [2018] A. Koga, S. Nakauchi, and J. Nasu, Phys. Rev. B 97, 094427 (2018).
  * Honerkamp and Hofstetter [2004] C. Honerkamp and W. Hofstetter, Phys. Rev. Lett. 92, 170403 (2004).
  * Chen and Wu [2021] G. Chen and C. Wu, arXiv preprint arXiv:2112.02630 (2021).
  * Xu  _et al._ [2019] H. Xu, Y. Wang, Z. Zhou, and C. Wu, Mott insulating states of the anisotropic $\mathrm{SU}(4)$ dirac fermions (2019), arXiv:1912.11791 [cond-mat.quant-gas] .
  * Cai and Wu [2020] Z. Cai and C. Wu, arXiv preprint arXiv:2003.08566 (2020).
  * Wang  _et al._ [2019] D. Wang, L. Wang, and C. Wu, Phys. Rev. B 100, 115155 (2019).
  * Zhou  _et al._ [2018] Z. Zhou, C. Wu, and Y. Wang, Phys. Rev. B 97, 195122 (2018).
  * Xu  _et al._ [2018a] S. Xu, J. T. Barreiro, Y. Wang, and C. Wu, Phys. Rev. Lett. 121, 167205 (2018a).
  * Kaul  _et al._ [2013] R. K. Kaul, R. G. Melko, and A. W. Sandvik, Annual Review of Condensed Matter Physics 4, 179 (2013), https://doi.org/10.1146/annurev-conmatphys-030212-184215 .
  * Nataf  _et al._ [2016] P. Nataf, M. Lajkó, A. Wietek, K. Penc, F. Mila, and A. M. Läuchli, Phys. Rev. Lett. 117, 167202 (2016).
  * Zhou  _et al._ [2014] Z. Zhou, Z. Cai, C. Wu, and Y. Wang, Phys. Rev. B 90, 235139 (2014).
  * Wang  _et al._ [2014a] D. Wang, Y. Li, Z. Cai, Z. Zhou, Y. Wang, and C. Wu, Phys. Rev. Lett. 112, 156403 (2014a).
  * Cai  _et al._ [2013a] Z. Cai, H.-H. Hung, L. Wang, and C. Wu, Phys. Rev. B 88, 125108 (2013a).
  * Cai  _et al._ [2013b] Z. Cai, H.-h. Hung, L. Wang, D. Zheng, and C. Wu, Phys. Rev. Lett. 110, 220401 (2013b).
  * Nataf and Mila [2014] P. Nataf and F. Mila, Phys. Rev. Lett. 113, 127204 (2014).
  * Hung  _et al._ [2011] H.-H. Hung, Y. Wang, and C. Wu, Phys. Rev. B 84, 054406 (2011).
  * Raczkowski and Assaad [2020] M. Raczkowski and F. F. Assaad, Phys. Rev. Research 2, 013276 (2020).
  * Kim  _et al._ [2019] F. H. Kim, F. F. Assaad, K. Penc, and F. Mila, Phys. Rev. B 100, 085103 (2019).
  * Da Liao  _et al._ [2022a] Y. Da Liao, X. Y. Xu, Z. Y. Meng, and Y. Qi, Phys. Rev. B 106, 075111 (2022a).
  * Da Liao  _et al._ [2022b] Y. Da Liao, X. Y. Xu, Z. Y. Meng, and Y. Qi, Phys. Rev. B 106, 115149 (2022b).
  * Zhu  _et al._ [2022] X. Zhu, Y. Huang, H. Guo, and S. Feng, Phys. Rev. B 106, 075109 (2022).
  * Beach  _et al._ [2009] K. S. D. Beach, F. Alet, M. Mambrini, and S. Capponi, Phys. Rev. B 80, 184401 (2009).
  * Harada  _et al._ [2003] K. Harada, N. Kawashima, and M. Troyer, Phys. Rev. Lett. 90, 117203 (2003).
  * Yao  _et al._ [2022] X.-P. Yao, R. L. Luo, and G. Chen, Phys. Rev. B 105, 024401 (2022).
  * Assaad and Grover [2016] F. F. Assaad and T. Grover, Phys. Rev. X 6, 041049 (2016).
  * Janssen and Seifert [2022] L. Janssen and U. F. P. Seifert, Phys. Rev. B 105, 045120 (2022).
  * Yamamoto and Kawakami [2022] K. Yamamoto and N. Kawakami, Universal description of dissipative tomonaga-luttinger liquids with $\mathrm{SU}$($\mathrm{N}$) spin symmetry: Exact spectrum and critical exponents (2022).
  * Wang  _et al._ [2022] S.-Y. Wang, D. Wang, and Q.-H. Wang, Transition from band insulator to mott insulator and formation of local moment in half-filled ionic $\mathrm{SU}$($\mathrm{N}$) hubbard model (2022).
  * Ibarra-García-Padilla  _et al._ [2021] E. Ibarra-García-Padilla, S. Dasgupta, H.-T. Wei, S. Taie, Y. Takahashi, R. T. Scalettar, and K. R. A. Hazzard, Phys. Rev. A 104, 043316 (2021).
  * Yoshida and Katsura [2021] H. Yoshida and H. Katsura, Phys. Rev. Lett. 126, 100201 (2021).
  * Romen and Läuchli [2020] C. Romen and A. M. Läuchli, Phys. Rev. Research 2, 043009 (2020).
  * Fromholz and Lecheminant [2020] P. Fromholz and P. Lecheminant, Phys. Rev. B 102, 094410 (2020).
  * Da Liao  _et al._ [2022c] Y. Da Liao, X. Y. Xu, Z. Y. Meng, and Y. Qi, Phys. Rev. B 106, 155159 (2022c).
  * DeSalvo  _et al._ [2010] B. J. DeSalvo, M. Yan, P. G. Mickelson, Y. N. Martinez de Escobar, and T. C. Killian, Phys. Rev. Lett. 105, 030402 (2010).
  * Gorshkov  _et al._ [2010] A. V. Gorshkov, M. Hermele, V. Gurarie, C. Xu, P. S. Julienne, J. Ye, P. Zoller, E. Demler, M. D. Lukin, and A. M. Rey, Nature Physics 6, 289 (2010).
  * Taie  _et al._ [2010] S. Taie, Y. Takasu, S. Sugawa, R. Yamazaki, T. Tsujimoto, R. Murakami, and Y. Takahashi, Phys. Rev. Lett. 105, 190401 (2010).
  * Krauser  _et al._ [2012] J. S. Krauser, J. Heinze, N. Fläschner, S. Götze, O. Jürgensen, D.-S. Lühmann, C. Becker, and K. Sengstock, Nature Physics 8, 813 (2012).
  * Taie  _et al._ [2012] S. Taie, R. Yamazaki, S. Sugawa, and Y. Takahashi, Nature Physics 8, 825 (2012).
  * Pagano  _et al._ [2014] G. Pagano, M. Mancini, G. Cappellini, P. Lombardi, F. Schäfer, H. Hu, X.-J. Liu, J. Catani, C. Sias, M. Inguscio, and L. Fallani, Nature Physics 10, 198 (2014).
  * Hofrichter  _et al._ [2016] C. Hofrichter, L. Riegger, F. Scazza, M. Höfer, D. R. Fernandes, I. Bloch, and S. Fölling, Phys. Rev. X 6, 021030 (2016).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011), https://www.pnas.org/doi/pdf/10.1073/pnas.1108174108 .
  * Chichinadze  _et al._ [2022] D. V. Chichinadze, L. Classen, Y. Wang, and A. V. Chubukov, Phys. Rev. Lett. 128, 227601 (2022).
  * Xu  _et al._ [2018b] X. Y. Xu, K. T. Law, and P. A. Lee, Phys. Rev. B 98, 121406 (2018b).
  * Da Liao  _et al._ [2019] Y. Da Liao, Z. Y. Meng, and X. Y. Xu, Phys. Rev. Lett. 123, 157601 (2019).
  * You and Vishwanath [2019] Y.-Z. You and A. Vishwanath, npj Quantum Materials 4, 16 (2019).
  * Li  _et al._ [2015a] Z.-X. Li, Y.-F. Jiang, and H. Yao, Phys. Rev. B 91, 241117 (2015a).
  * Li  _et al._ [2016] Z.-X. Li, Y.-F. Jiang, and H. Yao, Phys. Rev. Lett. 117, 267002 (2016).
  * Wang  _et al._ [2015] L. Wang, Y.-H. Liu, M. Iazzi, M. Troyer, and G. Harcos, Phys. Rev. Lett. 115, 250601 (2015).
  * Wei  _et al._ [2016] Z. C. Wei, C. Wu, Y. Li, S. Zhang, and T. Xiang, Phys. Rev. Lett. 116, 250601 (2016).
  * Assaad and Evertz [2008] F. Assaad and H. Evertz, in  _Computational Many-Particle Physics_ (Springer, 2008) pp. 277–356.
  * Blankenbecler  _et al._ [1981] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Phys. Rev. D 24, 2278 (1981).
  * Li and Yao [2019] Z.-X. Li and H. Yao, Annual Review of Condensed Matter Physics 10, 337 (2019), https://doi.org/10.1146/annurev-conmatphys-033117-054307 .
  * Read and Sachdev [1991] N. Read and S. Sachdev, Phys. Rev. Lett. 66, 1773 (1991).
  * Wen [1991] X. G. Wen, Phys. Rev. B 44, 2664 (1991).
  * Chen  _et al._ [2010] X. Chen, Z.-C. Gu, and X.-G. Wen, Phys. Rev. B 82, 155138 (2010).
  * Wen [1990] X.-G. Wen, International Journal of Modern Physics B 4, 239 (1990).
  * Li  _et al._ [2017] Z.-X. Li, Y.-F. Jiang, S.-K. Jian, and H. Yao, Nature communications 8, 1 (2017).
  * Su  _et al._ [1979] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. 42, 1698 (1979).
  * Fradkin and Hirsch [1983] E. Fradkin and J. E. Hirsch, Phys. Rev. B 27, 1680 (1983).
  * Wang  _et al._ [2014b] L. Wang, P. Corboz, and M. Troyer, New Journal of Physics 16, 103008 (2014b).
  * Li  _et al._ [2015b] Z.-X. Li, Y.-F. Jiang, and H. Yao, New Journal of Physics 17, 085003 (2015b).
  * Li  _et al._ [2020] B.-H. Li, Z.-X. Li, and H. Yao, Phys. Rev. B 101, 085105 (2020).
  * Zhang  _et al._ [2020] Y.-X. Zhang, H.-M. Guo, and R. T. Scalettar, Phys. Rev. B 101, 205139 (2020).
  * Cai  _et al._ [2021] X. Cai, Z.-X. Li, and H. Yao, Phys. Rev. Lett. 127, 247203 (2021).
  * Cai  _et al._ [2022] X. Cai, Z.-X. Li, and H. Yao, Phys. Rev. B 106, L081115 (2022).
  * Götz  _et al._ [2022] A. Götz, S. Beyl, M. Hohenadler, and F. F. Assaad, Phys. Rev. B 105, 085151 (2022).
  * Feng  _et al._ [2022] C. Feng, B. Xing, D. Poletti, R. Scalettar, and G. Batrouni, Phys. Rev. B 106, L081114 (2022).
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Annals of Physics 16, 407 (1961).
  * Oshikawa [2000] M. Oshikawa, Phys. Rev. Lett. 84, 3370 (2000).
  * Hastings [2004] M. B. Hastings, Phys. Rev. B 69, 104431 (2004).
  * Yao and Oshikawa [2020] Y. Yao and M. Oshikawa, Phys. Rev. X 10, 031008 (2020).
  * Assaad  _et al._ [2022] F. F. Assaad, M. Bercx, F. Goth, A. Götz, J. S. Hofmann, E. Huffman, Z. Liu, F. P. Toldin, J. S. E. Portela, and J. Schwab, SciPost Phys. Codebases , 1 (2022).

## Supplemental Material

###  .1 I. The details of projector Quantum Monte Carlo simulation

We employ the method of projector QMC to investigate the ground-state
properties of the SU($N$) SSH model at the anti-adiabatic limit, as described
in Eq. (1) in the main text. In Projector QMC, the ground-state expectation
value of the observable is evaluated by projecting a trial wave function
$\left|\psi_{T}\right>$ along the imaginary time direction:

$\langle{\hat{O}}\rangle=\frac{\left<\psi_{0}\right|\hat{O}\left|\psi_{0}\right>}{\langle\psi_{0}\mid\psi_{0}\rangle}=\lim_{\Theta\rightarrow\infty}\frac{\left<\psi_{T}\right|e^{-\frac{\Theta
H}{2}}\hat{O}e^{-\frac{\theta
H}{2}}\left|\psi_{T}\right>}{\left<\psi_{T}\right|e^{-\Theta
H}\left|\psi_{T}\right>}$ (S1)

Here $\left|\psi_{T}\right>$ is the trial wave function and $\Theta$ is the
projection parameter. The key point of the algorithm is intrinsically unbiased
against the choice of $\left|\psi_{T}\right>$ assuming that the trivial wave
function is not orthogonal to the exact ground state, namely
$\langle\psi_{T}\mid\psi_{0}\rangle\neq 0$, which is generically satisfied for
quantum many-body models in finite systems. In our simulation, we choose
$\left|\psi_{T}\right>$ as the ground-state wave function of the non-
interacting part of the model under consideration. The projection parameter
$\Theta$ in Eq. (S1) is chosen to be $\Theta=40$ for most simulations in our
study. For the simulations of model Eq. (1) for $N=4$, we choose a larger
projection parameter $\Theta=60$ to determine the ground-state properties of
competing ordered phases and the possible quantum spin liquid phase. We have
checked that the value of $\Theta$ is sufficient large to access the ground-
state properties of the model.

Similar to the finite-temperature determinant QMC algorithm, in projector QMC
Trotter decomposition is implemented to discretize the imaginary time $\Theta$
in spacing $\Delta_{\tau}=\Theta/L_{\tau}$, introducing a Trotter error
scaling as $\Delta_{\tau}^{2}$. In our simulation, we choose
$\Delta_{\tau}=0.1$ for most cases and $\Delta_{\tau}=0.05$ for the simulation
of model Eq. (1) for $N=4$. Of course, We have checked the convergence of the
results against further decreasing the value of $\Delta_{\tau}$. To decouple
the four-fermion interacting terms in the SU($N$) model Eq. (1), we perform
discrete Hubbard-Stratonovich transformation via introducing a classical
space-time dependent auxiliary field after the procedure of Trotter
decomposition:

$\displaystyle
e^{\frac{\Delta_{\tau}J}{2N}(\sum_{\alpha}c^{\dagger}_{i\alpha}c_{j\alpha}+H.c.)^{2}}=\sum_{l=\pm
1,\pm
2}\gamma(l)e^{\sqrt{\frac{\Delta_{\tau}J}{2N}}\eta(l)(\sum_{\alpha}c^{\dagger}_{i\alpha}c_{j\alpha}+H.c.)}+O(\Delta_{\tau}^{4})$
(S2)

with the four-valued parameters introduced as: $\gamma(\pm 1)=1+\sqrt{6}/3$,
$\gamma(\pm 2)=1-\sqrt{6}/3$, $\eta(\pm 1)=\pm\sqrt{2(3-\sqrt{6})}$, $\eta(\pm
2)=\pm\sqrt{2(3+\sqrt{6})}$. This scheme of Hubbard-Stratonovich
transformation keeps the SU($N$) symmetry, and is utilized to decouple the
SU($N$) four-fermion interactions.

To characterize various long-range orderings, we evaluate structure factors
for the corresponding order parameters. For model Eq. (1), the dominant
orderings are AFM spin order breaking spin rotational symmetry and VBS order
breaking lattice transnational symmetry. The structure factors for spin and
VBS orders with system size $N_{s}=L\times L$ are defined as:

$S_{\rm{spin}}(L,\vec{k})=\frac{1}{N_{s}^{2}N^{2}}\sum_{i,j,\alpha,\beta}\langle\hat{S}_{\alpha\beta}(i)\hat{S}_{\beta\alpha}(j)\rangle
e^{i\vec{k}\cdot(\vec{r_{i}}-\vec{r_{j}})}$ (S3)
$\begin{split}&S_{{\rm{VBS}}}(L,\vec{k},\delta)=\frac{1}{N_{s}^{2}N^{2}}\sum_{i,j}\langle\hat{K}_{\delta}(i)\hat{K}_{\delta}(j)\rangle
e^{i\vec{k}\cdot(\vec{r_{i}}-\vec{r_{j}})},\end{split}$ (S4)

where $\vec{k}$ is the momentum of structure factor in consideration. For spin
order,
$\hat{S}_{\alpha\beta}(i)=c^{\dagger}_{i\alpha}c_{i\beta}-\frac{\delta_{\alpha\beta}}{N}\sum^{N}_{\gamma=1}c^{\dagger}_{i\gamma}c_{i\gamma}$($\alpha,\beta=1,\cdots,N$)
is SU($N$) spin operator on site $i$, representing the generators of SU($N$)
group. Because of the relation $\sum_{\gamma}S_{\gamma\gamma}=0$, there are
$N^{2}-1$ independent generators for the SU($N$) group. AFM order on square
lattice denotes the spin density wave ordering with momentum
$\vec{k}=(\pi,\pi)$, hence the definition of structure factor for AFM order
reads: $S_{\rm{AFM}}(L)=S_{\rm{spin}}(L,\vec{k}=(\pi,\pi))$. For VBS structure
factor,
$\hat{K}_{\delta}(i)=\sum_{\alpha}c^{\dagger}_{i\alpha}c_{i+\delta\alpha}+H.c.$
is the kinetic operator on bond in $\delta=x,y$ direction. For Eq. (1), two
distinct VBS orderings are dominant with varying interaction strength $J$ and
fermion flavor number $N$, namely sVBS and cVBS. The sVBS order carries
momentum $\vec{k}=(\pi,\pi)$. In our simulation, we evaluate structure factor
for sVBS order averaging over $\hat{x}$ and $\hat{y}$-direction, which is
defined as
$S_{\rm{sVBS}}=\frac{1}{2}(S_{\rm{VBS}}(L,\vec{k}=(\pi,\pi),\hat{x})+S_{\rm{VBS}}(L,\vec{k}=(\pi,\pi),\hat{y}))$.
The cVBS order carries momentum $\vec{k}=(\pi,0)$ for kinetic term in
$\hat{x}$-direction and momentum $\vec{k}=(0,\pi)$ for kinetic term in
$\hat{y}$-direction. The structure factor for the cVBS order is defined as
$S_{\rm{cVBS}}=\frac{1}{2}(S_{\rm{VBS}}(L,\vec{k}=(\pi,0),\hat{x})+S_{\rm{VBS}}(L,\vec{k}=(0,\pi),\hat{y}))$.
VBS order breaks the lattice $\mathbb{Z}_{4}$ symmetry, while AFM breaks spin
$SU(2)$ rotational symmetry. The order parameter characterizing SSB is given
by the corresponding structure factor at thermodynamic limit:
$\Delta^{2}={\rm{lim}}_{L\rightarrow\infty}S(L)$.

###  .2 II. Derivations of single-particle and spin gap from time-displaced
correlation function

To understand the spectral properties of the model, we evaluate the single-
particle excitation gap extracted from the imaginary-time displaced Green’s
function. The single-particle gap denotes the minimum energy cost to extract
one electron from the system, corresponding to the spectral gap detected in
the photoemission experiment. Moreover, we extract the spin gap
$\Delta_{\rm{s}}(\textbf{k})$ from imaginary-time displaced spin-spin
correlation function. The spin gap is expected to vanish in the AFM ordered
phase owing to the existence of Goldstone mode.

To extract single-particle gap, we compute SU($N$) imaginary-time displaced
single-particle Green’s function:

$\displaystyle G(\vec{k},\tau)=\sum_{\alpha}\left\langle
c_{\vec{k}\alpha}(\tau)c^{\dagger}_{\vec{k}\alpha}(0)\right\rangle=\frac{1}{N^{2}_{s}}\sum_{\alpha}\sum_{ij}\left\langle
c_{i\alpha}(\tau)c^{\dagger}_{j\alpha}(0)\right\rangle
e^{i(\vec{R}_{i}-\vec{R}_{j})\cdot\vec{k}}$ (S5)

where $\vec{k}$ is the momentum under consideration and $N_{s}$ is the total
number of lattice site. The asymptotic behaviour at large $\tau$ gives rise to
the estimation of single-particle gap at momentum $\vec{k}$:
$G(\vec{k},\tau)\sim e^{-\Delta_{\rm{sp}}(\vec{k})\tau}$. Similarly, spin gap
is extracted from the imaginary-time displaced correlation function of spin
operator:

$\displaystyle G_{s}(\vec{k},\tau)=\sum_{\alpha\beta}\left\langle
S_{\alpha\beta}(\vec{k},\tau)S_{\beta\alpha}(\vec{k},0)\right\rangle$ (S6)

where $S_{\alpha\beta}(i)$ is SU($N$) spin operator defined as
$S_{\alpha\beta}(i)=c^{\dagger}_{i\alpha}c_{i\beta}-\frac{\delta_{\alpha\beta}}{N}\sum^{N}_{\gamma=1}c^{\dagger}_{i\gamma}c_{i\gamma}$
and $S_{\alpha\beta}(i\tau)=e^{\tau H}S_{\alpha\beta}(i)e^{-\tau H}$.
$\vec{k}$ is the momentum under consideration and $N_{s}$ is the total number
of lattice sites. At sufficient large $\tau$, the spin correlation function
scales as $G_{s}(\vec{k},\tau)\sim e^{-\Delta_{\rm{spin}}(\vec{k})\tau}$.
Hence, the single-particle gap and spin gap at a given momentum are available
to extract from the single-particle Green’s function and spin correlation
function at large imaginary time $\tau$ respectively. In our simulation, we
present the results of single-particle gap and spin gap with the minimum value
of $\Delta_{\rm{sp}}(\vec{k})$ and $\Delta_{\rm{spin}}(\vec{k})$ in the whole
Brillouin zone. For the results of spin gaps for $N=3,4$ as shown in main
text, the minimum spin gap $\Delta_{\rm{spin}}(\vec{k})$ occurs at
$\vec{k}=(\pi,\pi)$ in the whole interaction regime under consideration. For
single-particle gap, we compute local real-space imaginary-time displaced
Green’s function equivalent to the average of $G(\vec{k},\tau)$ in the whole
Brillouin zone: $G(i,\tau)=\frac{1}{N_{s}}\sum_{i,\alpha}\left\langle
c_{i\alpha}(\tau)c^{\dagger}_{i\alpha}(0)\right\rangle=\frac{1}{N_{s}}\sum_{\vec{k}}G(\vec{k},\tau)$.
At sufficient long imaginary time, $G(i,\tau)$ is determined by the minimum
value of $\Delta_{\rm{sp}}(\vec{k})$ in the whole Brillouin zone, scaling as
$G(i,\tau)\sim e^{-\Delta_{\rm{sp}}\tau}$ where $\Delta_{\rm{sp}}$ is the spin
gap defined as the minimum gap in the whole Brillouin zone. Hence, we extract
single-particle gap from the scaling behaviour of $G(i,\tau)$ at long
imaginary time.

###  .3 III. Additional QMC results for $N=3$

In this section, we present additional QMC results of Eq. (1) for $N=3$,
substantiating the conclusion that the ground state is AFM long-ranged ordered
phase for $N=3$.

In the maintext, we provide convincing evidence that the ground state
possesses AFM long-range order. Here we present the results of VBS ordering,
confirming the absence of VBS long-range order in the ground state for $N=3$.
For $N=3$, the peaked momentum of VBS structure factor is $(\pi,\pi)$, hence
we evaluate correlation-length ratio of staggered VBS order versus $J$ for
various system sizes, as shown in Fig. S1(a). Additionally, we compute the
single-particle gap at several values of $J$ for $L=8,10,12,14,16$, and fit
the results using second-polynomial function of $1/L$, the intercepts of which
give rise to the single-particle gap at thermodynamic limit. The results of
single-particle gap are depicted in Fig. S1(b), clearly demonstrating the
single-particle gaps are opened in the ground state for $N=3$.

###  .4 IV. Additional QMC results for $N=4$

In the maintext, we present QMC numerical results for VBS order and establish
the ground-state phase diagram for $N=4$. In the weakly coupling regime, the
ground state is a sVBS ordered phase, whereas in the strongly coupling regime
the dominant VBS ordering is columnar order carrying momentum
$(\pi,0)/(0,\pi)$. More remarkably, there exists an intermediate disordered
phase preserving all the lattice and Hamiltonian symmetries between the
staggered and columnar VBS ordered phases in the weakly and strongly coupling
regime. In this section, we present additional data of QMC simulation for
$N=4$, to exclude the possibility for the existence of long-range ordered
phases. The results of AFM correlation-length ratio and structure factors are
shown in Fig. S2(a) and (b), respectively. The correlation-length ratio of AFM
order decreases with linear system size $L$ in the whole interacting parameter
regime under consideration, indicating the AFM order is short-ranged for
$N=4$. In addition, the fitting results of AFM structure factors further
confirm the absence of AFM long-range order in the ground state for $N=4$.

For even value of $N$, the model Eq. (1) on bipartite lattice features an
additional partial particle-hole $Z_{2}$ symmetry exists at half filling:
$c_{i\alpha}\rightarrow c^{\dagger}_{i\alpha}(-1)^{i_{x}+i_{y}}$ for even
$\alpha$. Under such partial particle-hole transformation, the diagonal spin
operator $\sum_{\alpha}c^{\dagger}_{i\alpha}c_{i\alpha}(-1)^{\alpha}$ is
transformed to the density operator
$n_{i}=\sum_{\alpha}c^{\dagger}_{i\alpha}c_{i\alpha}$. The off-diagonal spin
operator $c^{\dagger}_{i\alpha}c_{i\alpha+1}$ ($\alpha$ is odd) is transformed
to staggered pairing operator
$c^{\dagger}_{i\alpha}c^{\dagger}_{i\alpha+1}(-1)^{(i_{x}+i_{y})}$.
Consequently, the AFM order parameters are degenerate with CDW order parameter
$O_{\rm{CDW}}=\frac{1}{N_{s}}\sum_{i}n_{i}(-1)^{i_{x}+i_{y}}$ and on-site
superconducting order parameter
$O_{\rm{SC}}=\frac{1}{N_{s}}\sum_{i}c^{\dagger}_{i\alpha}c^{\dagger}_{i\alpha+1}$.
Because we have unambiguously shown the absence of AFM order in the ground
state for $N=4$, the CDW and on-site singlet SC pairing order are also short-
ranged in the ground state of model Eq. (1) for $N=4$.

Figure S1: (a) Correlation-length ratio of sVBS order parameters
$R_{\rm{sVBS}}$ as a function of interacting strength $J$ for $N=3$.
Correlation-length ratio $R_{\rm{sVBS}}$ decreases with linear system size
$L$, indicating the staggered VBS long-range order is absent. (b)The results
of single-particle gap as a function of $1/L$ for $J=2,3,4$ and $N=3$. The
extrapolation of single-particle gap as a function of $1/L$ to thermodynamic
limit is performed, which unambiguously reveal the finite single-particle gap
exists in the AFM ordered phase for $N=3$.  Figure S2: (a) Correlation-length
ratio of AFM order parameters $R_{\rm{AFM}}$ as a function of interacting
strength $J$ for $N=4$. Correlation-length ratio $R_{\rm{AFM}}$ decreases with
linear system size $L$, indicating the absence of AFM long-range order. (b)
The finite-size scaling analysis for AFM structure factor as a function of
$1/L$ for $J=2.0,2.5,3.0$ and $N=4$. The extrapolated value at thermodynamic
limit $1/L\rightarrow\infty$ within error bar, further confirming the absence
of AFM long-range order in the ground state for $N=4$.

###  .5 V. Staggered VBS ordered phase in the large $N$ limit

For model Eq. (1), in the large-$N$ limit, the saddle-point approximation,
namely mean-field approximation, is able to give rise to the exact solution
since the fluctuation effect beyond saddle-point approximation is suppressed
by $\frac{1}{N}$. Here we perform self-consistent mean-field calculation to
investigate the ground-state properties of the model. The interaction term
reads:

$\displaystyle H_{J}=-\frac{J}{2N}\sum_{\langle
ij\rangle}(\sum_{\alpha}c^{\dagger}_{i\alpha}c_{j\alpha}+h.c.)^{2}$ (S7)

We consider SU($N$) symmetric mean-field decoupling of the interaction in
terms of bond hopping mean-field value:

$\displaystyle\varphi(i,\hat{\delta})=\langle
c^{\dagger}_{i\alpha}c_{i+\hat{\delta}\alpha}+h.c.\rangle,$ (S8)

resulting in the mean-field Hamiltonian:

$H^{\rm{MF}}=-\sum_{i\alpha}(t+J\varphi(i,\hat{x}))(c^{\dagger}_{i\alpha}c_{i+\hat{x}\alpha}+h.c.)-\sum_{i\alpha}(t+J\varphi(i,\hat{y}))(c^{\dagger}_{i\alpha}c_{i+\hat{y}\alpha}+h.c)+\sum_{i}\frac{J}{2}(\varphi^{2}(i,\hat{x})+\varphi^{2}(i,\hat{y}))$
(S9)

Since the mean-field Hamiltonian $H^{\rm{MF}}$ is quadratic and can be
diagonalized straightforwardly. Thus, we can solve the mean-field equation Eq.
(S8) self-consistently. Notice that given the $SU(N)$ symmetric ansatz Eq.
(S8), the results of mean-field solution is independent on the value of $N$.
We consider two mean-field ansatz corresponding to the staggered and columnar
VBS orders. For the sVBS order:

$\displaystyle\langle c^{\dagger}_{i\alpha}c_{i+\hat{x}\alpha}+h.c.\rangle$
$\displaystyle=$ $\displaystyle\phi^{s}_{0}+\phi_{\rm{stag}}(-1)^{x+y}$
$\displaystyle\langle c^{\dagger}_{i\alpha}c_{i+\hat{y}\alpha}+h.c.\rangle$
$\displaystyle=$ $\displaystyle\phi^{s}_{1}$ (S10)

For the cVBS order:

$\displaystyle\langle c^{\dagger}_{i\alpha}c_{i+\hat{x}\alpha}+h.c.\rangle$
$\displaystyle=$ $\displaystyle\phi^{c}_{0}+\phi_{\rm{colu}}(-1)^{x}$
$\displaystyle\langle c^{\dagger}_{i\alpha}c_{i+\hat{y}\alpha}+h.c.\rangle$
$\displaystyle=$ $\displaystyle\phi^{c}_{1}$ (S11)

where $\phi_{\rm{stag}}$ and $\phi_{\rm{colu}}$ are order parameters
characterizing staggered and columnar VBS orders respectively. $\phi^{s}_{0}$,
$\phi^{s}_{1}$, $\phi^{c}_{0}$ and $\phi^{c}_{1}$ represent uniform parts in
the mean-field ansatz. The values of these parameters are solved self-
consistently. The results of sVBS and cVBS order parameters as a function of
interaction strength $J$ are shown in Fig. S3(a). The sVBS order parameter is
finite in the whole interacting parameter regime under consideration,
increasing with interacting strength $J$. On the other hand, for cVBS order,
the self-consistent mean-field calculations show that order parameter remains
zero in a large interacting parameter regime $J\leq 4$. Given by the mean-
field solutions, the ground-state energy of the staggered and columnar VBS
phases for mean-field Hamiltonian are easily achieved. We present the results
of ground-state energies for staggered and columnar VBS phases in Fig. S3(b),
unambiguously indicating staggered VBS phase is energetically favored in the
whole interaction parameter regime under consideration. Hence, the mean-field
calculations provide convincing evidence that the ground state of model Eq.
(1) is sVBS ordered phase in the large-$N$ limit.

Figure S3: Large $N$ mean field results for SU(N) SSH model. (a) Staggered and
columnar VBS order parameters as a function of interaction strength
$J\in[0,10]$. (b) Results of ground-state energies for staggered and columnar
VBS phase as a function of interaction strength $J\in[0,10]$.


# Frequency Hopping Joint Radar-Communications with Hybrid Sub-pulse Frequency
and Duration Modulation

Linh Manh Hoang, , J. Andrew Zhang, , Diep N. Nguyen, , and Dinh Thai Hoang
The authors are with the School of Electrical and Data Engineering, University
of Technology Sydney, NSW 2007. e-mail<EMAIL_ADDRESS>(Andrew.Zhang, Diep. Nguyen<EMAIL_ADDRESS>
###### Abstract

Frequency-hopping (FH) joint radar-communications (JRC) can offer excellent
security for integrated sensing and communication systems. However, existing
JRC schemes mainly embed information using only the sub-pulse frequencies and
hence the data rate is limited. In this paper, we propose to use both sub-
pulse frequencies and durations for information modulation, leading to higher
communication data rates. For information demodulation, we propose a novel
scheme by using the time-frequency analysis (TFA) technique and a ‘you only
look once’ (YOLO)-based detection system. As such, our system does not require
channel estimation, simplifying the transmission signal frame design.
Simulation results demonstrate the effectiveness of our scheme, and show that
it is robust against the Doppler shift and timing offset between the
transceiver and the communication receiver.

###### Index Terms:

Joint radar-communications, waveform design, YOLO.

##  I Introduction

The frequency-hopping (FH) joint radar-communications systems (JRC) [1, 2, 3,
4, 5] have become more and more prevalent in both military and civil
applications, e.g., airborne, shipborne, ground-based combat systems, the
connected autonomous vehicle (CAV). The most popular FH JRC systems include
multiple-input multiple-output orthogonal frequency division multiplexing
(MIMO-OFDM), multi-Carrier AgilE phaSed Array Radar (CAESAR), and FH-MIMO[2].

A key challenge in FH JRC is to embed/modulate data (onto frequency sub-
pulses) at the transceiver and demodulate the signal at the receiver. For data
embedding, the most dominant technique is the index modulation (IM), which
utilizes frequency selection/combination and/or antenna selection/permutation
for data representation. Specifically, the data bits are embedded by selecting
different sets of sub-pulse frequency (i.e., frequency combination) and
allocating them to different antennas (i.e., antennas permutation). For data
demodulation, the optimal demodulator can be based on the maximum likelihood
principle [3], while sub-optimal methods, which require lower computation
complexities, are based on compressive sensing (CS) or discrete Fourier
transform (DFT) [2]. However, the performance of these methods depends on the
accuracy of channel estimation. Unfortunately, a long training sequence
necessary for accurate channel estimation is not always feasible. This is
because a long training sequence reduces the communication fraction over the
whole time frame, thus decreasing the data bit rate. In particular, the data
embedding schemes using both sub-pulse frequency combinations and antenna
permutations require complex demodulation techniques that are prone to
demodulation error. On the other hand, the data embedding scheme using only
the sub-pulse frequency combination has a limited data transmission rate.

This paper proposes novel techniques to embed and demodulate data in an FH JRC
system to increase the data rate and reduce the demodulation error. For data
embedding, we use both sub-pulse frequency and duration, therefore increasing
the data transmission rate compared to only using the sub-pulse frequency. For
data demodulation, we propose a novel scheme based on the signal’s time-
frequency image (TFI) and a YOLO-based detection system. This demodulation
scheme does not require channel estimation and is robust to the Doppler shift
and the timing offset between the transceiver and the communication receiver.
Moreover, the proposed data embedding and demodulation schemes are spatially
flexible and not limited to the sidelobe of the transmit beampattern, since
the data is not embedded by utilizing the phase or amplitude of the
beampattern sidelobe.

##  II System and Signals

Figure 1: FH joint radar-communications system.

###  II-A System Overview

We consider a ground-to-air JRC system, where the JRC transceiver is located
on the ground and the communication receiver is mounted on an aircraft or a
high-altitude platform (HAP), as illustrated in Fig. 1. For sensing, the JRC
transceiver transmits the radio frequency (RF) signal and processes the
reflected RF signal to determine target’s position and velocity. For data
transmission, the communication receiver obtains the RF signal from the JRC
transceiver and demodulate the embedded data. In such a system, the line-of-
sight (LoS) path dominates, and non-LoS paths are negligible. At the
transceiver or communication receiver, the received RF signal is down-
converted into the intermediate frequency (IF) signal and sampled at sampling
frequency $f_{\rm s}$ to generate a discrete-time signal as

$\displaystyle r[k]$ $\displaystyle=s[k]+w[k].$ (1)

Here, $k$ is the sample index which increases every $T_{\rm s}=1/f_{\rm s}$,
$s[k]$ and $r[k]$ are the transmitted and received signal samples,
respectively, and $w[k]$ is the complex additive white Gaussian noise (AWGN).
The JRC system uses a pulse wave (PW) waveform with a duty cycle $D=\triangle
t_{\rm w}/\triangle t_{\rm r},D\ll 1$, where $\triangle t_{\rm w}$ and
$\triangle t_{\rm r}$ are the pulse width and pulse repetition intervals,
respectively. The received signal at the communication receiver is then used
to generate the Choi-Williams distribution time-frequency image (CWD-TFI).
Next, the CWD-TFI is preprocessed in the Preprocessing block to generate the
training and testing data for the YOLO detection system. The training process
of the YOLO detection system is demonstrated by the dashed line in Fig. 1
using the training data. The trained YOLO detection system and the testing
CWD-TFIs are used to generate the annotated CWD-TFIs, which are used to
demodulate the embedded data.

###  II-B Signals Model

The JRC system uses FH signals with each signal pulse being divided into
$N_{\rm f}$ sub-pulses (i.e., hops) with different frequencies and durations.
Let ${\bf{F}}=(f_{1},f_{2},...,f_{N_{\rm f}})$ denote the FH sequence and
${\bf{T}}=(\triangle t_{1},\triangle t_{2},...,\triangle t_{{N_{\rm f}}})$
denote the sub-pulse duration sequence. The JRC signals in each pulse is
represented by

$\displaystyle s[k]=\sum_{n=1}^{N_{\rm f}}$ $\displaystyle Ae^{j2\pi
f_{n}kT_{\rm s}}{\rm{rect}}\Big{[}\frac{kT_{\rm s}-\sum_{i=0}^{n-1}\triangle
t_{i}}{\triangle t_{n}}\Big{]},$ (2) $\displaystyle\text{where }A\text{ is the
}$ signal amplitude, $\displaystyle{\text{rect}}[k]$
$\displaystyle\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}}$
$\displaystyle\text{and }\triangle t_{0}$
$\displaystyle\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}}0.$

##  III Data Embedding Schemes

In this work, we use the Costas arrays to generate the FH sequence of each
signal pulse in one of the data embedding schemes. Using the Costas array, an
FH signal possesses a narrow peak at the origin of its ambiguity function (AF)
and low sidelobes elsewhere, which is desirable for the radar measurement
(i.e., range and speed) accuracy. Details of the Costas arrays, their
construction, and characteristics can be found in [6, Ch. 5].

We propose two data embedding techniques, namely Random and Costas-based
schemes. Unlike existing studies that only utilize sub-pulse frequencies to
convey information, both sub-pulse frequencies and durations are used in our
proposed schemes. We first design the codebooks $\mathcal{S}$ and
$\mathcal{S}_{\rm C}$ of the Random and Costas-based schemes, respectively.
Then, for each embedding scheme, data bits are embedded by mapping each
element of the corresponding codebook to a symbol. Let
${\bf{\mathcal{F}}}\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}}[(1,2,...,N_{\rm
f})\times f_{\rm f}]$ be the selection set of the sub-pulse frequency, where
$f_{\rm f}$ is the fundamental frequency. Let
${\bf{\mathcal{D}}}\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}}(\triangle_{1},\triangle_{2},...,\triangle_{N_{{\rm{T}}}})$
be the selection set of the sub-pulse duration, where $N_{{\rm{T}}}$ is the
number of sub-pulse duration selections. The codebooks $\mathcal{S}$ and
$\mathcal{S}_{\rm C}$ are designed as follows.

  * •

Random scheme: $\mathcal{S}$ is designed by selecting each sub-pulse frequency
from the frequency set ${\bf{\mathcal{F}}}$, and selecting each sub-pulse
duration from the duration set ${\bf{\mathcal{D}}}$ as

$\displaystyle f_{i}$ $\displaystyle\in{\bf{\mathcal{F}}},\triangle
t_{i}\in{\bf{\mathcal{D}}},\forall i\in(1,2,...,N_{\rm f}),$ (3)
$\displaystyle f_{i}$ $\displaystyle\neq f_{i+1},\forall i\in(1,2,...,N_{\rm
f}-1).$ (4)

Here, condition (4) means every two consecutive sub-pulse frequencies are not
equal, which is necessary for the data demodulation technique that will be
presented in Section V.

  * •

Costas-based scheme: Unlike the Random scheme, for designing $\mathcal{S}_{\rm
C}$, each sub-pulse frequency is not individually selected. Instead, the whole
FH sequence of each signal pulse is selected from the Costas-based FH sequence
defined as follows.

Let ${\bf{\mathcal{A}}_{\rm
C}}\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}}({\bf{A}}_{1},{\bf{A}}_{2},...,{\bf{A}}_{N_{\rm
C}})$ be the set of Costas array of length $N_{\rm f}$, where
${\bf{A}}_{i},\forall i\in(1,2,...,N_{\rm C})$ is an $\mathbb{N}^{N_{\rm
f}\times 1}$ Costas array, and $N_{\rm C}$ is the number of Costas sequence
with length of $N_{\rm f}$. Then, let ${\bf{\mathcal{F}}_{\rm
C}}\mathrel{\ensurestackMath{\stackon[1pt]{=}{\scriptstyle\Delta}}}[({\bf{A}}_{1},{\bf{A}}_{2},...,{\bf{A}}_{N_{\rm
C}})\times f_{\rm f}]$ be the selection set of the Costas-based FH sequence.
Then, the FH sequence of each signal pulse is selected from
${\bf{\mathcal{F}}_{\rm C}}$. On the other hand, similar to the Random scheme,
each sub-pulse duration is selected from the sub-pulse duration set
${\bf{\mathcal{D}}}$ as

$\displaystyle{\bf{F}}$ $\displaystyle\in{\bf{\mathcal{F}}_{\rm C}},\triangle
t_{i}\in{\bf{\mathcal{D}}},\forall i\in(1,2,...,N_{\rm f}).$ (5)

Note that for the Costas-based scheme, because the FH sequence is generated
using the Costas array, condition (4) is satisfied. The dimensions of
$\mathcal{S}$ and $\mathcal{S}_{\rm C}$ are given by

$\displaystyle|\mathcal{S}|$ $\displaystyle=\underbrace{N_{{\rm{T}}}^{N_{\rm
f}}}_{\text{sub-pulse duration selections}}\times\underbrace{N_{\rm f}(N_{\rm
f}-1)^{(N_{\rm f}-1)}}_{\text{sub-pulse frequency selections}},$ (6)
$\displaystyle|\mathcal{S}_{\rm C}|$
$\displaystyle=\underbrace{N_{{\rm{T}}}^{N_{\rm f}}}_{\text{sub-pulse duration
selections}}\times\underbrace{N_{\rm C}}_{\text{sub-pulse frequency
selections}}.$ (7)

As shown in (6), the number of selections of sub-pulse frequencies is $N_{\rm
f}(N_{\rm f}-1)^{(N_{\rm f}-1)}$ instead of $N_{\rm f}^{N_{\rm f}}$ due to the
requirement in (4).

It follows from (6) and (7) that the maximum number of bits $C$ and $C_{\rm
C}$ that can be embedded in each pulse of the two schemes are given by

$\displaystyle C=$
$\displaystyle\lfloor\log_{2}(|\mathcal{S}|)\rfloor=\lfloor\log_{2}(N_{{\rm{T}}}^{N_{\rm
f}}N_{\rm f}(N_{\rm f}-1)^{(N_{\rm f}-1)})\rfloor,$ $\displaystyle C_{\rm C}=$
$\displaystyle\lfloor\log_{2}(|\mathcal{S}_{\rm
C}|)\rfloor=\lfloor\log_{2}(N_{{\rm{T}}}^{N_{\rm f}}N_{\rm C})\rfloor,$ (8)

where $\lfloor.\rfloor$ denotes the floor function, determining the closest
smaller integer. As can be seen, the Random scheme achieves a higher
transmission rate than the Costas-based scheme.

##  IV Sensing at the JRC Transceiver

This section demonstrates that by using the sub-pulse frequencies and duration
for data embedding as presented in the previous section, the detection
probability of the sensing function is not affected. Let $P_{\rm D}$ and
$P_{\rm FA}$ denote the detection and false alarm probabilities of the JRC
transceiver, respectively. Let $E_{\rm{p}}$ be the signal energy per each
pulse of the JRC signal and ${\sigma_{\rm w}^{2}}$ be the noise variance. The
energy-to-noise ratio (ENR) of the JRC signal is given by

ENR $\displaystyle=10\log_{10}(\frac{E_{\rm{p}}}{{\sigma_{\rm w}^{2}}}).$ (9)

To detect the transmitted signal $s[k]$ buried in the reflected RF signal
$r[k]$, the JRC transceiver performs the matched filtering. Let $R_{s}$
denotes the output of the matched filter. Based on Neyman-Pearson criterion,
the transceiver decides a detection if $R_{s}$ exceeds a threshold $\gamma$

$\displaystyle R_{s}$ $\displaystyle=\sum_{k=0}^{N_{\rm
p}-1}r[k]s[k]>\gamma,\qquad\qquad\qquad\qquad$ (10) where $N_{\rm p}$ is the
number of signal sample per pulse, $\displaystyle\gamma$
$\displaystyle=\sqrt{\sigma_{\rm w}^{2}E_{\rm{p}}}Q^{-1}(P_{\rm FA}),$ (11)
$\displaystyle Q(x)$
$\displaystyle=\int_{x}^{\infty}\frac{1}{\sqrt{2\pi}}\exp(-\frac{1}{2}t^{2})\;dt\quad$
(12)

is the complement of the cumulative distribution function (CDF) of the
standard Gaussian distribution, and $Q^{-1}(x)$ denotes the inverse function
of $Q(x)$. The characteristic of $P_{\rm D}$ is given in Theorem 1 below.

Theorem 1: The value of $P_{\rm D}$ is independent from the JRC signal
waveform and only depends on $P_{\rm FA}$ and ENR.

Proof: From [7, Ch. 4], $P_{\rm D}$ can be expressed by $P_{\rm FA}$ and ENR
as

$\displaystyle P_{\rm D}$ $\displaystyle=Q\Big{[}Q^{-1}(P_{\rm
FA})-\sqrt{10^{{\rm ENR}/10}}\;\Big{]}.$ (13)

Accordingly, $P_{\rm D}$ only depends on $P_{\rm FA}$ and ENR, and is
independent from the JRC signal waveform. $\blacksquare$

Therefore, varying the sub-pulse frequencies and durations does not affect the
detection performance of the sensing function, as long as the signal energy
per each pulse is fixed. Note that, however, varying the sub-pulse durations
may cause some increase in the peak-to-average power ratio (PAPR), which
reduces the efficiency of the power amplifier at the front end of the JRC
transceiver and communication receiver. Therefore, the sub-pulse selection set
${\bf{\mathcal{D}}}$ should be carefully designed to avoid an excessive PAPR
value.

##  V Data Demodulation at the Communication Receiver

This section presents the detailed procedure of the proposed data demodulation
schemes and its advantages, including being unaffected by channel estimation
error (i.e., the channel estimation is not required), and robustness against
the Doppler shift or timing offset between the transceiver and the
communication receiver.

###  V-A Time-Frequency Analysis Technique

Figure 2: CWD-TFI of a Costas signal with ${\bf{F}}=[(4,2,5,1,3)\times f_{\rm
f}]$, ${\bf{T}}=(5,2,4,2,3)\times 1/(25f_{\rm s})$, $f_{\rm f}=f_{\rm s}/16$,
and a starting point of $520$.

The CWD is one of the bilinear TFA techniques that is more robust against
noise than the linear ones, thanks to the auto-correlation operation in their
formulas. Moreover, the CWD does not have the limitation of the cross-term
effect that can obscure essential features of the signal. Therefore, in this
work, we use the CWD to generate the CWD-TFI, which represents the time-
frequency characteristic of the signal. Specifically, the CWD is
mathematically given by[8]

$\displaystyle\text{CW}[k,n]$
$\displaystyle\\!=\\!2\sum_{\tau=-\infty}^{\infty}W_{N}(\tau)e^{-j2\pi
n\tau/N}\qquad\qquad\qquad\qquad\qquad$
$\displaystyle\times\Bigg{[}\sum_{\mu=-\infty}^{\infty}W_{M}(\mu)\frac{1}{\sqrt{4\pi\tau^{2}/\sigma}}e^{-\frac{\mu^{2}}{4\tau^{2}/\sigma}}$
$\displaystyle\times r[k+\mu+\tau]r^{*}[k+\mu-\tau]\Bigg{]},$ (14)

where $k$ and $n$ are the time index and frequency index, respectively,
$W_{N}(\tau)$ is a symmetrical window function with non-zero values for the
range of $-N/2\leq\tau\leq N/2$, while $W_{M}(\mu)$ is a rectangular window
function with a value of $1$ for the range of $-M/2\leq\mu\leq M/2$, and
$\sigma$ is the scaling factor. Similar to [9], the value of $N$ and $M$ are
often chosen as a power of $2$ to minimize the computational cost (by
exploiting the fast Fourier transform algorithm). Furthermore, $\sigma=1$ is
used to achieve a balance between frequency resolution and the cross-terms
reduction in the CWD.

###  V-B CWD-TFI Preprocessing

To generate the CWD-TFI, similar to [9], we capture $N_{\rm s}=2048$
consecutive signal samples, and assume these $N_{\rm s}$ samples contain a
complete signal pulse, and the signal pulse starts from an arbitrary point
within the $N_{\rm s}$ samples. Fig. 2 illustrates the CWD-TFI of a Costas
signal with a hopping sequence of ${\bf{F}}=[(4,2,5,1,3)\times f_{\rm f}]$, a
sub-pulse duration sequence of ${\bf{T}}=(5,2,4,2,3)\times 1/(25f_{\rm s})$, a
fundamental frequency of $f_{\rm f}=f_{\rm s}/16$, and a starting point of
$520$. For brevity, the signal uses a Costas array to generate the FH sequence
is referred to as a Costas signal. As can be seen, the CWD-TFI clearly
demonstrates the time-frequency characteristic of the Costas signal by
horizontal lines of different lengths and vertical positions.

Fig. 3 shows the CWD-TFI processing procedure used to generate the input to
the YOLO detection system. First, using the nearest neighbor interpolation,
the CWD-TFI of $N_{\rm s}\times N_{\rm s}$ pixels is resized to $L\times L$
pixels, which is the input size of the YOLO detection system, where $L=500$
pixels is the width of the resized CWD-TFIs. Then, the signal sub-pulse
locations in the CWD-TFI are specified by the bounding boxes (i.e., white-
color rectangles in Fig. 3), as a requirement for the training of the YOLO.
Note that the localization of the bounding boxes is not needed in the CWD-TFI
processing procedure in the testing phase, because the bounding boxes are
parts of the output of the trained YOLO detection system.

###  V-C “You only look once” (YOLO) Object Detection System

Among the object detection algorithm for determining the objects types and
their locations in the image, YOLO [10] is well-known for its real-time
detection speed, which is critical in tactical applications, such as military
JRC systems. The real-time detection speed is achieved by splitting the input
image into a grid of cells and predicting an object and its bounding box on
each cell, instead of performing region proposal as in other object detection
systems. Details of the YOLO detection system can be found in [10].

Fig. 4 shows the YOLO’s output for the CWD-TFI of a Costas signal with an FH
sequence of ${\bf{F}}=[(4,2,5,1,3)\times f_{\rm f}]$, a sub-pulse duration
sequence of ${\bf{T}}=(4,3,5,4,2)\times 1/(25f_{\rm s})$, a fundamental
frequency of $f_{\rm f}=f_{\rm s}/16$, and a starting point of $259$. As can
be seen, the output contains not only the signal type (i.e., “Costas”) but
also the locations of the signal objects (i.e., the white-color rectangular
bounding boxes) and the probabilities of the detection.

Figure 3: CWD-TFI preprocessing procedure for the Costas signal in Fig. 2.
Figure 4: YOLO’s output for a Costas signal with ${\bf{F}}=[(4,2,5,1,3)\times
f_{\rm f}]$, ${\bf{T}}=(4,3,5,4,2)\times 1/(25f_{\rm s})$, $f_{\rm f}=f_{\rm
s}/16$, and a starting point of $259$.

###  V-D Data Demodulation Technique

As aforementioned, the output from the YOLO detection system contains the
bounding boxes, which determine the position, horizontal, and vertical width
of the signal objects in the CWD-TFI. The $i$th left-most bounding box,
corresponds to the $i$th sub-pulse, is localized by a tuple of $(x^{\rm
min}_{i},x^{\rm max}_{i},y^{\rm min}_{i},y^{\rm max}_{i})$. Recall that for
the CWD-TFI, $L$ pixels in the vertical axis correspond to a frequency range
from $0$ to $f_{\rm s}/2$, while $L$ pixels in the horizontal axis correspond
to a time duration of $N_{\rm s}\times(1/fs)$, the frequency and duration of
the $i$th sub-pulse can be calculated by

$\displaystyle f_{i}$ $\displaystyle=\frac{f_{\rm s}(y^{\rm min}_{i}+y^{\rm
max}_{i})}{4L},$ $\displaystyle\triangle t_{i}$ $\displaystyle=\frac{N_{\rm
s}(x^{\rm max}_{i}-x^{\rm min}_{i})}{Lf_{\rm s}}.$ (15)

By comparing the combination of $f_{i}$ and $\triangle t_{i}$ values with the
selections in the codebook $\mathcal{S}$ or $\mathcal{S}_{\rm C}$, the
embedded symbol in each signal pulse can be determined.

It is worth noting that this demodulation technique does not require any
knowledge about the channel between the JRC transceiver and the communication
receiver, and hence the channel estimation is not necessary. Moreover, since
our data demodulation technique is based on the CWD-TFI, its performance is
not affected by the Doppler shift or timing offset between the transceiver and
the communication receiver. This is because the Doppler shift or timing offset
merely creates vertical or horizontal shifts of the whole set of lines,
respectively. On the other hand, the shapes of the lines in the CWD-TFI remain
unchanged.

##  VI Performance Evaluation

TABLE I: Parameters values. Parameter |  Description |  Value
---|---|---
$N_{\rm f}$ |  Number of sub-pulses per pulse |  5
$N_{\rm C}$ |  Number of Costas array |  40
${\bf{\mathcal{D}}}$ |  Sub-pulse duration selection set |  $(1,1.5,2,2.5,3)\times 1/(25f_{\rm s})$
$f_{\rm f}$ |  Fundamental frequency |  $f_{\rm s}/16$
$L$ |  Resized CWD-TFI width (pixels) |  500

Table. I shows the parameter values used to produce the simulation results.
Note that the number of Costas array $N_{\rm C}$ is not a freely set value,
but depends on $N_{\rm f}$, the number of sub-pulses per pulse. For instance,
$N_{\rm C}=40$ for $N_{\rm f}=5$. The list of available Costas arrays and
their construction can be found in [6, Ch. 5]. To train the YOLO detection
system, for each data embedding scheme, we generate $2700$ signals with SNR
ranging from $-6$ dB to $10$ dB and a step size of $2$ dB. The $2700$ signals
are divided into a training set of $2160$ signals ($80\%$ of the total) and a
validation set of $540$ signals (20% of the total). For the testing data, we
generate $3300$ signals for each embedding scheme, with SNR ranging from $-10$
dB to $10$ dB and a step size of $2$ dB.

Figure 5: Detection probability of the JRC sensing function for different
signal waveforms, $P_{\rm FA}$, and ENR values. Figure 6: Maximum number of
bits per pulse for different data embedding schemes.

Fig. 5 shows the detection probability of the JRC sensing function for
different signal waveforms (i.e., Random scheme, Costas-based scheme, and
fixed sub-pulse length), $P_{\rm FA}$, and ENR values. The results are
obtained by averaging $10,000$ independent trials. As can be seen, for each
$P_{\rm FA}$ value, the three curves almost overlap, meaning the detection
probability only depends on $P_{\rm FA}$ and ENR but not the signal waveform.
Therefore, varying the sub-pulse durations does not affect the detection
capability of the sensing function of the FH RJC system, as long as the signal
energy per each pulse is fixed.

Fig. 6 demonstrates the maximum number of bits that can be transmitted over a
signal pulse of the two proposed embedding schemes and those in [3, 5]. As can
be seen, the Random scheme can embed a significantly higher amount of data
than the Costas-based scheme, while the Costas-based scheme can convey a
slightly higher number of data bits than the scheme in [3, 5]. That is because
the schemes in this paper use both sub-pulse durations and frequencies to
embed signals. On the other hand, the technique in [3, 5] focuses on using
antenna permutations for data embedding. Therefore, the proposed techniques
are more suitable when the number of antennas is limited, while the technique
in [3, 5] is suitable in the case of an extensive number of antennas is
available.

Figure 7: Symbol error rate (SER) of different embedding and demodulation
schemes.

Fig. 7 illustrates the symbol error rate (SER) of the two proposed schemes and
those of the techniques in [5] as functions of the signal-to-noise ratio
(SNR). As can be seen, the two proposed schemes have lower SERs compared to
those of the techniques in [5] for all SNRs. Distinctively, at low SNRs (i.e.,
$\leq-2$ dB), the two proposed schemes have similar SERs, and are about $4$ dB
and $10$ dB better than those of the FH and BPSK techniques in [5],
respectively. This is because the proposed techniques use the CWD, which
contains the auto-correlation operation, to generate CWD-TFIs of the signals.
As such, the proposed techniques are more robust against noise compared to the
other techniques in [5].

##  VII Conclusions

We have proposed novel approaches for embedding and demodulating information
bits in FH JRC systems. For data embedding, we have used both sub-pulse
durations and frequencies to increase the data transmission rate. For data
demodulation, to reduce the error rate, we have proposed a CWD-TFI and YOLO-
based demodulation scheme that does not require channel estimation and is
robust against Doppler shift and timing offset between the JRC transceiver and
communication receiver. Simulation results have shown that our proposed
techniques achieve higher data transmission and lower symbol error rates than
the existing ones.

## References

  * [1] A. Hassanien, M. G. Amin, E. Aboutanios, and B. Himed, “Dual-function radar communication systems: A solution to the spectrum congestion problem,” _IEEE Signal Process. Mag._ , vol. 36, no. 5, pp. 115–126, Sep 2019.
  * [2] J. A. Zhang, F. Liu, C. Masouros, R. W. Heath, Z. Feng, L. Zheng, and A. Petropulu, “An overview of signal processing techniques for joint communication and radar sensing,” _IEEE J. Sel. Topics Signal Process._ , vol. 69, no. 2, pp. 1295–1315, Nov 2021.
  * [3] T. Huang, N. Shlezinger, X. Xu, Y. Liu, and Y. C. Eldar, “MAJoRCom: A dual-function radar communication system using index modulation,” _IEEE Trans. Signal Process_ , vol. 68, pp. 3423–3438, May 2020.
  * [4] K. Wu, J. A. Zhang, X. Huang, Y. J. Guo, and R. W. Heath, “Waveform design and accurate channel estimation for frequency-hopping MIMO radar-based communications,” _IEEE Trans. on Commun._ , vol. 69, no. 2, pp. 1244–1258, Oct 2020.
  * [5] W. Baxter, E. Aboutanios, and A. Hassanien, “Dual-function MIMO radar-communications via frequency-hopping code selection,” in _Proc. Asilomar Conf. Signals, Syst., Comput._ IEEE, 2018, pp. 1126–1130.
  * [6] N. Levanon and E. Mozeson, _Radar signals_. Hoboken, NJ, USA: John Wiley & Sons, 2004.
  * [7] S. M. Kay, _Fundamentals of statistical signal processing: detection theory, volume II_. Englewood Cliffs, NJ, USA: Prentice-Hall, 1998.
  * [8] H.-I. Choi and W. J. Williams, “Improved time-frequency representation of multicomponent signals using exponential kernels,” _IEEE Trans. Acoust., Speech, Signal Process._ , vol. 37, no. 6, pp. 862–871, Jun 1989.
  * [9] L. M. Hoang, M. Kim, and S.-H. Kong, “Automatic recognition of general LPI radar waveform using SSD and supplementary classifier,” _IEEE Trans. Signal Process._ , vol. 67, no. 13, pp. 3516–3530, May 2019.
  * [10] J. Redmon, S. Divvala, R. Girshick, and A. Farhadi, “You only look once: Unified, real-time object detection,” in _IEEE CVPR_ , 2016, pp. 779–788.