Simulation of higher-order topological phases in 2D spin-phononic crystal networks

The periodic driving is known to render effective Hamiltonian in which specific terms can be adiabatically eliminated. In particular, the periodic driving can be used to trigger nonequilibrium topological behavior in a trivial setup, which offers an efficient tool to simulate topological phases in quantum systems Goldman and Dalibard (2014); Gómez-León and Platero (2013). We consider a periodic driving quantum system with $H(t)=H(t+T)$, characterized by time period $T=2\pi/\omega$. In this case we can introduce Floquet theorem to investigate long-time dynamics of the system, as developed in Ref. Eckardt and Anisimovas (2015). With the Floquet-Bloch ansatz, the time-dependent Schrodinger equation be given by

$$i\hbar d_{t}\left|\psi_{\alpha}(t)\right\rangle=H(t)\left|\psi_{\alpha}(t)\right\rangle,$$

(7)

where

$$\left|\psi_{\alpha}(t)\right\rangle=\left|\phi_{\alpha}(t)\right\rangle e^{-i\epsilon_{\alpha}t/\hbar}=e^{-i\epsilon_{\alpha}t/\hbar}\sum_{m}e^{-im\omega t}\phi_{m}.$$

(8)

$\left|\psi_{\alpha}(t)\right\rangle$ is the so-called Floquet eigenstate, and $\epsilon_{\alpha}$ is the quasienergy with band index $\alpha$. $\phi_{\alpha}(t)=\phi_{\alpha}(t+T)$ denotes the time-periodic Floquet eigenmode, which can be constructed by a complete set of orthonormal basis state $\phi_{m}$. With respect to the basis $\left|\psi_{\alpha}(t)\right\rangle$, the system can be effectively described by the Hamiltonian

$$H_{eff}^{mn}=\frac{1}{T}\int_{0}^{T}dte^{i(m-n)\omega t}H(t)\text{.}$$

(9)

The effective Hamiltonian is the time-average of the Hamiltonian $H(t)$ in a driving period, which is the core of Floquet theorem. Note that the Floquet state in time-periodically driven systems is analogous to the Bloch state in spatially periodic systems.

In a recent work Li et al. (2020b), we proposed a periodic driving protocol to simulate topological phases with a color center-phononic crystal system. By applying a standing wave field between the two lowest sublevels ($|g\rangle,|e\rangle$) of the SiV center, we get the Floquet engineering of the spin-spin interactions, resulting in the well-known SSH-type Hamiltonian. Here we consider a fundamentally different driving protocol, which allows us to selectively control the spin-spin interactions. What is more important, the resulting spin-spin interactions possess chiral symmetry and support rich quantum phases associated with topological invariants. The time-periodic driving has form Pérez-González et al. (2019)

$$H_{driv}(t)=\sum_{j}\hbar V_{j}f(t)\sigma_{j}^{z},$$

(10)

where $\sigma_{j}^{z}=|e\rangle_{j}\langle e|-|g\rangle_{j}\langle g|$ is the Pauli operator component. $f(t)$ denotes the standard square-wave function

	$\displaystyle f(t)$	$\displaystyle=$	$\displaystyle-1\text{ \ \ \ for }t\in[0,\frac{T}{2}],$
	$\displaystyle f(t)$	$\displaystyle=$	$\displaystyle 1\text{ \ \ \ \ for }t\in[\frac{T}{2},T].$		(11)

$V_{j}$ denotes the on-site potential

$$V_{j}=\left\{\begin{array}[]{l}b_{0}+\frac{(a_{0}+b_{0})}{2}(j-1)\ \ \ \ j=1,3,5,7,...\\ \frac{(a_{0}+b_{0})}{2}j\ \ \ \ \ \ \ \ \ \ \ \ \ \ j=2,4,6,8,...\end{array}\right..$$

(12)

This stair-like form offers alternating potential difference between two adjacent spins, i.e., $V_{j}-V_{j-1}=a_{0}$ and $V_{j+1}-V_{j}=b_{0}$ are staggered along the spin array.

We first consider the interaction of the periodic driving and the 1D spin array, but the case of 2D will be studied in the next section. Now we transform the total Hamiltonian $H_{1D}=H_{array}+H_{driv}(t)$ into the interaction picture, with the unitary operator $U(t)=e^{-i\int_{0}^{t}d\tau H_{driv}(\tau)/\hbar}$. After the unitary transformation, we obtain

	$\displaystyle\sigma_{eg}^{j}$	$\displaystyle\rightarrow$	$\displaystyle e^{i\Delta_{j}(t)\sigma_{j}^{z}}\sigma_{eg}^{j}e^{-i\Delta_{j}(t)\sigma_{j}^{z}}=\sigma_{eg}^{j}e^{2i\Delta_{j}(t)},$
	$\displaystyle\sigma_{ge}^{j}$	$\displaystyle\rightarrow$	$\displaystyle e^{i\Delta_{j}(t)\sigma_{j}^{z}}\sigma_{ge}^{j}e^{-i\Delta_{j}(t)\sigma_{j}^{z}}=\sigma_{ge}^{j}e^{-2i\Delta_{j}(t)},$		(13)

with

	$\displaystyle\Delta_{j}(t)$	$\displaystyle=$	$\displaystyle V_{j}\int_{0}^{t}d\tau f(\tau)$		(14)
		$\displaystyle=$	$\displaystyle V_{j}\int_{0}^{t}d\tau[\underset{n\neq 0}{\sum}\frac{1}{n\pi i}(e^{-in\pi}-1)e^{in\omega\tau}],$

where we expanded $f(t)$ into its Fourier series. In the interaction picture, the total Hamiltonian has the form as

$$H_{1D}=\underset{i,j}{\sum}\hbar J_{ij}(t)\sigma_{eg}^{i}\sigma_{ge}^{j},$$

(15)

where $J_{ij}(t)=J_{ij}e^{2i(\Delta_{i}(t)-\Delta_{j}(t))}$ is the hopping rate with a temporal periodicity, $J_{ij}(t)=J_{ij}(t+T)$. The Floquet components of the Hamiltonian ($14$) read

$$H_{1D}^{mn}=\underset{i,j}{\sum}\hbar J_{ij}^{mn}\sigma_{eg}^{i}\sigma_{ge}^{j},$$

(16)

$$J_{ij}^{mn}=\frac{1}{T}\int_{0}^{T}dtJ_{ij}(t)e^{i(m-n)\omega t}.$$

(17)

For the time-periodically driven system, $H_{1D}^{mn}$ can be expressed by the Floquet-Magnus expansion. In the high-frequency regime $\omega\gg J_{i,j}$, it is a good approximation to neglect the rapid oscillation of the external driving Rahav et al. (2003); Eckardt and Anisimovas (2015); Mikami et al. (2016). As a result, the spin-spin interaction can be given by the zeroth-order expansion term

$$\mathcal{J}_{ij}=J_{ij}\frac{i\omega}{2\pi(V_{i}-V_{j})}(e^{-i2\pi(V_{i}-V_{j})/\omega}-1).$$

(18)

As for the SSH model, the interparticle interaction is characterized by staggering hopping amplitudes. Thus, the two nearest-neighbor spins can be grouped into a unit cell and classified as odd and even spins, as we proposed in Ref. Li et al. (2020b). Likewise, here we consider the phononic-mediated spin-spin interaction as a bipartite lattice of the form $ABABAB$, with

	$\displaystyle A_{n}$	$\displaystyle=$	$\displaystyle\sigma_{ge}^{j}\ \ \ j=1,3,5,7,...,$
	$\displaystyle B_{n}$	$\displaystyle=$	$\displaystyle\sigma_{ge}^{j}\ \ \ j=2,4,6,8,....$		(19)

Based on this definition, we rewrite the renormalized hopping amplitude $\mathcal{J}_{ij}$. For simplicity, here we introduce $n_{l}$ to label the spins at the $l$ site of the $n$th cell, $l=A$ or $B$. In general, there are two types of interparticle hopping. For the even-neighbor hopping, which describes the spin-spin interaction of the same sublattice, as shown in Fig. $2(a)$. The potential difference are $V_{i}-V_{j}=\pm m(a_{0}+b_{0})$, with $m=n_{l}^{\prime}-n_{l}$. In consequence, the even-neighbor spin-spin hopping rate can be written as

$$\mathcal{J}_{n_{l},n_{l}^{\prime}}=\frac{iJ_{n_{l},n_{l}^{\prime}}}{\mp 2\pi qm}(e^{\pm i2\pi qm}-1),$$

(20)

where $q=(a_{0}+b_{0})/\omega$, and “$\pm$” correspond to  the coupling to the right and left spins, respectively. From Eq. (20), the even-neighbor hopping is always zero if we assign $q=1,2,3,...$. Thus we conclude that the even-neighbor hopping can be suppressed by tuning the parameters $a_{0}$ and $b_{0}$. Note that the even-neighbor hopping is a detrimental source for the chiral symmetry Pérez-González et al. (2019).

For the odd-neighbor hopping, which describes the spin-spin interaction of the different sublattice. To better describe the physical picture of the spin-spin interaction, we further classify two kinds of odd-neighbor hopping. We first discuss the case with $n_{A}\leq n_{B}$, for which the schematic diagram is shown in Fig. $2(b)$. If we define $n_{B}=n_{A}+r$ $(r=0,1,2,...)$, the spin-spin interaction can be described by

	$\displaystyle\mathcal{J}_{n_{A},n_{B}}$	$\displaystyle=$	$\displaystyle-\frac{iJ_{n_{A},n_{B}}}{2\pi(qr+\frac{a_{0}}{\omega})}[e^{2i\pi(qr+\frac{a_{0}}{\omega})}-1],$
	$\displaystyle\mathcal{J}_{n_{B},n_{A}}$	$\displaystyle=$	$\displaystyle\frac{iJ_{n_{A},n_{B}}}{2\pi(qr+\frac{a_{0}}{\omega})}[e^{-2i\pi(qr+\frac{a_{0}}{\omega})}-1],$		(21)

where $\mathcal{J}_{n_{A},n_{B}}$ and $\mathcal{J}_{n_{B},n_{A}}$ describe the forward $(A\rightarrow B)$ and backward $(B\rightarrow A)$ hoppings, respectively. For the case with $n_{A}>n_{B}$, the schematic diagram is shown in Fig. $2(c)$. If we define $n_{B}=n_{A}-r^{\prime}$ $(r^{\prime}=1,2,3,...)$, the spin-spin interaction can be described by

	$\displaystyle\mathcal{J}_{n_{A},n_{B}}^{\prime}$	$\displaystyle=$	$\displaystyle\frac{iJ_{n_{A},n_{B}}}{2\pi(qr^{\prime}-\frac{a_{0}}{\omega})}[e^{-2i\pi(qr^{\prime}-\frac{a_{0}}{\omega})}-1],$
	$\displaystyle\mathcal{J}_{n_{B},n_{A}}^{\prime}$	$\displaystyle=$	$\displaystyle-\frac{iJ_{n_{A},n_{B}}}{2\pi(qr^{\prime}-\frac{a_{0}}{\omega})}[e^{2i\pi(qr^{\prime}-\frac{a_{0}}{\omega})}-1].$		(22)

Likewise, $\mathcal{J}_{n_{A},n_{B}}^{\prime}$ and $\mathcal{J}_{n_{B},n_{A}}^{\prime}$ represent the backward $(A\rightarrow B)$ and forward $(B\rightarrow A)$ hopping, respectively. From Eqs. (20) and (21), we can conclude

$$\mathcal{J}_{n_{A},n_{B}}=(\mathcal{J}_{n_{B},n_{A}})^{\ast},\mathcal{J}_{n_{A},n_{B}}^{\prime}=(\mathcal{J}_{n_{B},n_{A}}^{\prime})^{\ast}.$$

(23)

Unlike the case of the SSH model, the backward and forward hoppings of the odd-neighbor spin-spin interaction are not equal.

According to this bipartite solution, the Hamiltonian $H_{1D}$ can be rewritten as

	$\displaystyle H_{1D}$	$\displaystyle=$	$\displaystyle\underset{n,r,r^{\prime}}{\sum}\hbar(\mathcal{J}_{n_{A},n_{B}}A_{n}B_{n+r}^{\dagger}+\mathcal{J}_{n_{B},n_{A}}A_{n}^{\dagger}B_{n+r}$		(24)
			$\displaystyle+\mathcal{J}_{n_{A},n_{B}}^{\prime}A_{n}B_{n-r^{\prime}}^{\dagger}+\mathcal{J}_{n_{B},n_{A}}^{\prime}A_{n}^{\dagger}B_{n-r^{\prime}}).$

Here we neglected the even-neighbor hopping terms. By applying a particular periodic driving field to the SiV centers, we obtain the Floquet engineering of the spin-spin interactions with unique properties. In this case, the even-neighbor hopping is suppressed by tuning the parameters of the driving field, while the odd-neighbor hopping can be enhanced as needed. This scheme enforces the chiral symmetry which provides topological protection for the spin-spin interaction. In Fig. $3$, we numerically calculate the quasienergy spectrum as a function of $a_{0}$. We can see that all the eigenmodes are grouped into chiral symmetric pairs with opposite energies. For simplicity, here we express the bare spin-spin interaction as

$$J_{i,j}=\frac{g_{c}^{2}}{2\Delta_{BE}}e^{-|x_{i}-x_{j}|/L_{c}}=J_{0}e^{-|x_{i}-x_{j}|/L_{c}}.$$

(25)

Given that the band-gap mediated spin-spin interaction decays exponentially with the spin spacing, here only the first- and third- neighbor interactions are included.

The periodic driving protocol offers an effective method to investigate the topological character of the spin-phononic crystal system. To explore topological features of the Floquet engineering spin-spin system, we convert $H_{1D}$ to the momentum space. Considering periodic boundary conditions, we can make the Fourier transformation

$$O_{n}=\frac{1}{\sqrt{N}}\sum_{k}e^{ink}O_{k},(O=A,B)$$

(26)

where $k=2\pi m/N(m=1,2,...,N)$ is the wavenumber in the first Brillouin zone, and $A_{k}$ and $B_{k}$ are the momentum space operators. Defining the unitary operator $\psi(k)=\left(\begin{array}[]{cc}A_{k}&B_{k}\end{array}\right)^{T}$, the Hamiltonian $H_{1D}$ be expressed as

$$H_{1D}=\sum_{k}\psi(k)^{\dagger}H(k)\psi(k).$$

(27)

Then we obtain $2\times 2$ matrix form of the Hamiltonian in the $k$-space

$$H(k)=\hbar\left(\begin{array}[]{cc}0&f(k)\\ f^{\ast}(k)&0\end{array}\right).$$

(28)

with

$$f(k)=\underset{r,r^{\prime}}{\sum}(\mathcal{J}_{n_{B},n_{A}}e^{ikr}+\mathcal{J}_{n_{B},n_{A}}^{{}^{\prime}}e^{-ikr^{\prime}}).$$

(29)

Here $f(k)$ describes the coupling between the $A$ and $B$ spins in momentum space.

The dispersion relation can be obtained by solving the eigenvalue equation

$$H(k)\psi(k)=E(k)\psi(k),$$

(30)

using the fact that $H^{2}(k)=E^{2}(k)I$, with $I$ being the identity operator in the Hilbert space. Then we obtain the energy band structure as

$$E(k)=\pm\hbar|f(k)|,$$

(31)

$$\psi(k)=\frac{1}{\sqrt{2}}\binom{1}{\pm e^{-i\vartheta(k)}}.$$

(32)

$\psi(k)$ corresponds to the eigenfunctions for the lower and upper band, and $\vartheta(k)$ is defined as the argument of $f(k)$. Figs. 4(a1)-(f1) show the energy spectra for different driving field parameters, which are split into two branches and there exists a band gap between the lower and upper branches. It should be noticed that the band gap will be vanished at the critical point of topological phases.

The band gap structures are generally associated with topological properties of bulk-boundary correspondence. For the 1D Floquet engineering spin-spin system, we introduce the topological Zak phase Shen et al. (2018)

	$\displaystyle\mathcal{\varphi}_{Zak}$	$\displaystyle=$	$\displaystyle-i\overset{occ.}{\underset{j=1}{\sum}}\int_{0}^{2\pi}dk\psi^{\dagger}(k)\partial_{k}\psi(k)$		(33)
		$\displaystyle=$	$\displaystyle N_{occ.}\frac{1}{2}\int_{0}^{2\pi}dk\frac{d}{dk}\vartheta(k)$
		$\displaystyle=$	$\displaystyle\mathcal{W}\pi$

where $\mathcal{W}$ is the topological winding number, $N_{occ.}$ describes the number of occupied energy bands. Now we need to investigate the winding number of the system. Alternatively, $f(k)$ can be expressed in the form

$$f(k)=d(k)\cdot\sigma,$$

(34)

where $\sigma=(\sigma_{x},\sigma_{y},\sigma_{z})$ is the Pauli matrix, and $d(k)$ denotes a three-dimensional vector field

	$\displaystyle d_{x}(k)$	$\displaystyle=$	$\displaystyle\frac{1}{2}\underset{r,r^{\prime}}{\sum}(\mathcal{J}_{n_{B},n_{A}}e^{ikr}+\mathcal{J}_{n_{B},n_{A}}^{{}^{\prime}}e^{-ikr^{\prime}}+c.c.),$
	$\displaystyle d_{y}(k)$	$\displaystyle=$	$\displaystyle\frac{i}{2}\underset{r,r^{\prime}}{\sum}(\mathcal{J}_{n_{B},n_{A}}e^{ikr}+\mathcal{J}_{n_{B},n_{A}}^{{}^{\prime}}e^{-ikr^{\prime}}-c.c.),$
	$\displaystyle d_{z}(k)$	$\displaystyle=$	$\displaystyle 0.$		(35)

For general $2$-band topological insulators, owing to the periodicity of the momentum-space Hamiltonian, the path of the endpoint of $d(k)$ is a closed loop in the auxiliary space $(d_{x},d_{y})$ Asbóth et al. (2016). The topology of this loop can be characterized by an integer, the winding number

$$\mathcal{W}=\frac{1}{2\pi}\int_{0}^{2\pi}\mathbf{n}\times\partial_{k}\mathbf{n}dk,$$

(36)

where $\mathbf{n}=(n_{x},n_{y})=(d_{x},d_{y})/\sqrt{d_{x}^{2}+d_{y}^{2}}$ is the normalized vector. Here the winding number $\mathcal{W}$ counts the number of times the loop winds around the origin of the $dx$-$dy$ plane. Figs. 4(a2)-(f2) present the path of the endpoint of $d(k)$ on the $dx$-$dy$ plane. For different values of $a_{0}$, the winding number of the system exhibits four possible values, $-1$, $0$, $1$, $2$. According to Eq. (33), we can derive the relevant topological Zak phases directly. Furthermore, one can implement the topological phase transition in this SiV-phononic system by modulating the periodic driving.

From the numerical simulation results, we show rich quantum phases related to topological invariants. As for the generalized SSH model, a prototypical example to investigate topological properties in a trivial system, the 1D Zak phase has only two possible values $0$ or $\pi$. This work offers an effective scheme for studying topological phases induced by periodic driving. The distinct feature is that it enables to simulate higher-order topological phases and related topological phase transitions in topological trivial systems.

The existence of edge states at the boundary is a distinguished feature for topological insulator states. In the following, we first simulate the edge states in a 1D spin-phononic system. The core step is to look for the zero-energy eigenstates. Here we introduce the single-excited state

$$\psi=\underset{n}{\sum}(a_{n}A^{\dagger}_{n}+b_{n}B^{\dagger}_{n})\left|0\right\rangle,$$

(37)

where $a_{n}$ and $b_{n}$ are the amplitudes of occupying probability in the $n$th cell. $\left|0\right\rangle=\left|ggg...\right\rangle$ is the vacuum state, which describes that all spins stay in the ground state $|g\rangle$. In the single-excited state subspace, we can get the the zero-energy eigenstates by sloving

$$H_{1D}\underset{n}{\sum}(a_{n}A^{\dagger}_{n}+b_{n}B^{\dagger}_{n})\left|0\right\rangle=0.$$

(38)

There will be $2N$ equations for the amplitudes $a_{n}$ and $b_{n}$. Considering the boundary conditions, $b_{0}=a_{N+1}=0$. We can analytically derive the left and right zero-energy edge states, respectively.

To verify the model, we numerically simulate the energy spectrum and zero-energy eigenstates of the system. Figs. 5(a1)-(d1) show the eigenvalues for various parameter settings of the periodic driving. As for the non-topological regime, $\mathcal{W}=0$, there will be an energy band gap but no gapless modes appear. For the cases with $\mathcal{W}=1$ and $\mathcal{W}=-1$, there are two zero-energy eigenvalues. For the case with $\mathcal{W}=2$, there are four zero-energy eigenvalues. Correspondingly, we plot the zero-energy edge states in Figs. 5(a2)-(d2). We see that the wavefunctions are located at the vicinity of the array boundaries, which are the so-called topological edge states. In addition, the edge states only distribute at certain (odd or even) sites, which is related to the chiral symmetry of the system.

Now we proceed to generalize the above 1D results to 2D spin-phononic crystal networks. Here we consider adding two mutually perpendicular microwave fields to the color center arrays Zou and Liu (2017). The first one is a time-dependent microwave driving of frequency $\omega_{x}$ in the $x$ direction. The other is a time-dependent driving of frequency $\omega_{y}$ in the $y$ direction. These two periodic driving fields have the form

	$\displaystyle H_{driv}^{(x)}$	$\displaystyle=$	$\displaystyle\overset{2N}{\sum_{l=1}}\overset{2N}{\sum_{j=1}}\hbar V_{j,l}f_{x}(t)\sigma_{j,l}^{z},$
	$\displaystyle H_{driv}^{(y)}$	$\displaystyle=$	$\displaystyle\overset{2N}{\sum_{j=1}}\overset{2N}{\sum_{l=1}}\hbar V_{j,l}f_{y}(t)\sigma_{j,l}^{z}.$		(39)

$V_{j,l}=(V_{j},V_{l})$ describes the on-site potential in the 2D phononic network, the two components of which are in the form of Eq. $(12)$. $f_{x}(t)$ and $f_{y}(t)$ denote the square-wave function in the $x$ and $y$ directions, respectively.

Let us discuss the two directions separately. For the periodic driving spin arrays along the $x$ direction, the total Hamiltonian can be written as

$$H_{2D}^{(x)}=H_{array}^{(x)}+H_{driv}^{(x)}.$$

(40)

In the interaction picture, we introduce the unitary operator $U_{x}(t)=e^{-i\int_{0}^{t}d\tau H_{driv}^{(x)}/\hbar}$. After the unitary transformation, we obtain

$$\sigma_{eg}^{(j,l)}\rightarrow\sigma_{eg}^{(j,l)}e^{2i\Delta_{j}(t)},\sigma_{ge}^{(j,l)}\rightarrow\sigma_{ge}^{(j,l)}e^{-2i\Delta_{j}(t)},$$

(41)

with

$$\Delta_{j}(t)=V_{j}\int_{0}^{t}d\tau f(\tau).$$

(42)

While for the spin arrays along the $y$ direction, the total Hamiltonian is given by

$$H_{2D}^{(y)}=H_{array}^{(y)}+H_{driv}^{(y)}.$$

(43)

Similary, we introduce the unitary operator $U_{y}(t)=e^{-i\int_{0}^{t}d\tau H_{driv}^{(y)}/\hbar}$. After the unitary transformation, we obtain

$$\sigma_{eg}^{(j,l)}\rightarrow\sigma_{eg}^{(j,l)}e^{2i\Delta_{l}(t)},\sigma_{ge}^{(j,l)}\rightarrow\sigma_{ge}^{(j,l)}e^{-2i\Delta_{l}(t)},$$

(44)

with

$$\Delta_{l}(t)=V_{l}\int_{0}^{t}d\tau f(\tau).$$

(45)

In the following, analogous to the 1D case, we consider the bipartite interaction in both $x$ and $y$ directions. Then we get a 2D system with $N\times N$ unit cells. As depicted in Fig. $6$, there are four spins in each unit cell, which are labeled as {$A,B,C,D$}, respectively. In the regime $\omega\gg J_{i,j},J_{k,l}$, we derive the Floquet engineering spin-spin interactions along the $x$ and $y$ directions

	$\displaystyle H_{2D}^{(x)}$	$\displaystyle=$	$\displaystyle\underset{m}{\sum}\underset{n,r,r^{\prime}}{\sum}\hbar[\mathcal{J}_{n,n+r}(A_{n,m}B_{n+r,m}^{\dagger}+C_{n,m}D_{n+r,m}^{\dagger})$
			$\displaystyle+\mathcal{J}_{n,n-r^{\prime}}^{\prime}(A_{n,m}B_{n-r^{\prime},m}^{\dagger}+C_{n,m}D_{n-r^{\prime},m}^{\dagger})+H.c.],$
	$\displaystyle H_{2D}^{(y)}$	$\displaystyle=$	$\displaystyle\underset{n}{\sum}\underset{m,r,r^{\prime}}{\sum}\hbar[\mathcal{J}_{m,m+r}(A_{n,m}C_{n,m+r}^{\dagger}+B_{n,m}D_{n,m+r}^{\dagger})$		(46)
			$\displaystyle+\mathcal{J}_{m,m-r^{\prime}}^{\prime}(A_{n,m}C_{n,m-r^{\prime}}^{\dagger}+B_{n,m}D_{n,m-r^{\prime}}^{\dagger})+H.c.].$

For simplicity, we introduce $(n,m)$ to describe the position of each unit cell in the 2D spin-spin networks, with $n,m=1,2,\ldots,N$.

To simplify the model, here we suppose that the spin spacing $d_{x}=d_{y}$ and the periodic driving frequencies $\omega_{x}=\omega_{y}$. In this case, we can derive

	$\displaystyle\mathcal{J}_{n,n+r}$	$\displaystyle=$	$\displaystyle\mathcal{J}_{m,m+r},$
	$\displaystyle\mathcal{J}_{n,n-r^{\prime}}^{\prime}$	$\displaystyle=$	$\displaystyle\mathcal{J}_{m,m-r^{\prime}}^{\prime}.$		(47)

Thus we can define $s=n$ or $m$, and the two-dimensional Hamiltonian can be further integrated as

	$\displaystyle H_{2D}$	$\displaystyle=$	$\displaystyle\underset{r,r^{\prime}}{\sum}\underset{n,m}{\sum}\hbar[\mathcal{J}_{s,s+r}(A_{n,m}B_{n+r,m}^{\dagger}+C_{n,m}D_{n+r,m}^{\dagger}$		(48)
			$\displaystyle+A_{n,m}C_{n,m+r}^{\dagger}+B_{n,m}D_{n,m+r}^{\dagger})$
			$\displaystyle+\mathcal{J}_{s,s-r^{\prime}}^{\prime}(A_{n,m}B_{n-r^{\prime},m}^{\dagger}+C_{n,m}D_{n-r^{\prime},m}^{\dagger}$
			$\displaystyle+A_{n,m}C_{n,m-r^{\prime}}^{\dagger}+B_{n,m}D_{n,m-r^{\prime}}^{\dagger})+H.c.].$

To investigate the topological features in the 2D Floquet engineering spin-spin system, we convert the Hamiltonian $H_{2D}$ to the momentum space. Here we consider periodic boundary conditions along both the $x$ and $y$ directions. Then we apply the Fourier transformation to the four spins in a unit cell

$$O_{n,m}=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}e^{i(k_{x}n+k_{y}m)}O_{\mathbf{k}},(O=A,B,C,D)$$

(49)

where $\mathbf{k}=(k_{x},k_{y})$ is the wavenumber in the first Brillouin zone. If we define the unitary operator $\psi(\mathbf{k})=\left(\begin{array}[]{cccc}A_{\mathbf{k}}&B_{\mathbf{k}}&C_{\mathbf{k}}&D_{\mathbf{k}}\end{array}\right)^{T}$, the two-dimensional Hamiltonian can be rewritten as

$$H_{2D}=\sum_{\mathbf{k}}\psi^{\dagger}(\mathbf{k})H(\mathbf{k})\psi(\mathbf{k}).$$

(50)

Along with that we get $4\times 4$ matrix form of the Hamiltonian in the $\mathbf{k}$-space

$$H(\mathbf{k})=\hbar\left(\begin{array}[]{cccc}0&f(k_{x})&f(k_{y})&0\\ f^{\ast}(k_{x})&0&0&f(k_{y})\\ f^{\ast}(k_{y})&0&0&f(k_{x})\\ 0&f^{\ast}(k_{y})&f^{\ast}(k_{x})&0\end{array}\right).$$

(51)

$f(k_{x})$ and $f(k_{y})$ have the same form as Eq. (29), which describe the spin-spin couplings in the $x$ and $y$ directions, respectively.

Let us study the 2D dispersion relation by solving the eigenvalue equation

$$H(\mathbf{k})\psi(\mathbf{k})=E(\mathbf{k})\psi(\mathbf{k}),$$

(52)

then obtain

$$E(\mathbf{k})=\epsilon_{x}\hbar\left|f(k_{x})\right|+\epsilon_{y}\hbar\left|f(k_{y})\right|,$$

(53)

$$\psi(\mathbf{k})=\frac{1}{2}\left(\begin{array}[]{c}1\\ \epsilon_{x}e^{-i\vartheta_{x}(k_{x})}\\ \epsilon_{y}e^{-i\vartheta_{y}(k_{y})}\\ \epsilon_{x}\epsilon_{y}e^{-i[\vartheta_{x}(k_{x})+\vartheta_{y}(k_{y})]}\end{array}\right),$$

(54)

where $\epsilon_{i}=\pm 1$, $\vartheta_{i}(k_{i})=\arg[f(k_{i})],$ $i=x,y$. In Figure. 7(a), we numerically calculate the 2D energy spectrum in the momentum space. There are four energy bands since there are four spins in a unit cell. The lowest and highest bands are isolated, while the two middle bands are jointed at the edges of the Brillouin zone ($0,0$), ($\pm\pi,\pm\pi$), ($\mp\pi,\pm\pi$). According to Eq. (53), there exist two equal energy band gaps. When assigning suitable values of $a_{0}$, these four bands will be jointed together, and the band gaps vanished. This is a signature of topological phase transition.