Quantum Monte Carlo study of topological phases on a spin analogue of Benalcazar-Bernevig-Hughes model

Recently, a new family of higher-order topological insulators (HOTIs) have attracted extensive attentions due to their novel and fundamental physics Benalcazar et al. (2017a, b); Song et al. (2017); Langbehn et al. (2017); Khalaf (2018); Schindler et al. (2018a); Slager et al. (2015); Călugăru et al. (2019). An $n$th-order HOTI, like its conventional cousins, has gapless excitations but at even lower $(d-n)$-dimensional boundaries that are protected by higher-order topological invariants defined in the $d$-dimensional bulk. For instance, the electric quadrupole insulator, proposed by Benalcazar et al., features zero-energy modes at the corners protected by a quantized electric quadrupole moment in the bulk Benalcazar et al. (2017a, b). This new concept of higher-order topological protection on states of matter has been immediately generalized to semimetalsEzawa (2018a); Lin and Hughes (2018), superconductorsYan et al. (2018); Wang et al. (2018), and even non-Hermitian systemsKunst et al. (2018); Liu et al. (2019a); Ghatak and Das (2019), constituting rich higher-order topological phases. From the experimental aspect, HOTIs have been reported in several platforms, such as the electric circuitsImhof et al. (2018); Serra-Garcia et al. (2019); Ezawa (2018b), microwave resonatorsPeterson et al. (2018), classical mechanical metamaterialsMa et al. (2019); Fan et al. (2019), photonic and phononic crystalsZilberberg et al. (2018); Noh et al. (2018); Xue et al. (2019); Ni et al. (2019); Xie et al. (2018); Zhang et al. (2019); Xie et al. (2019); El Hassan et al. (2019); Chen et al. (2019); Mittal et al. (2019), although they are difficult to be found in electric systems Schindler et al. (2018b); Liu et al. (2019b); Sheng et al. (2019); Chen et al. (2020).

Up to date, however, the main attention paid on the higher-order topological phases is within the framework of topological band theory, their closely related counterparts taking account the strong correlation effects are still less known. Introducing correlations will unexpectedly bring new and different properties to the system. For instance, it is found that the quantized electric quadrupole insulators are robust against weak interactions but are driven to antiferromagnetic insulators by strong enough correlations Peng et al. (2020). Besides, the electron correlations will generalize the bulk-boundary correspondence, i.e., there emerge gapless corner modes only in spin excitations while the single-particle excitations remain gapped Kudo et al. (2019). The bosonic counterparts of HOTIs have also been proposed in a two-dimensional superlattice Bose-Hubbard model and dimerized spin systemsBibo et al. (2020); Dubinkin and Hughes (2019); You et al. (2018), and the representative signatures such as fractional corner charges are explicitly demonstrated using the density matrix renormalization groupBibo et al. (2020).

To characterize the topological properties in system with strong correlations, there are several proposals for the possible many-body topological invariant. One recipe is the Green’s function formalism that is widely used in quantum Monte Carlo (QMC) simulations. Based on the zero-frequency Green’s function, the nested Wilson loop method can be directly applied to obtain the many-body topological invariants Peng et al. (2020); Wang and Zhang (2012); Wang and Yan (2013). Another method directly extends the charge polarization of one-dimensional (1D) systemsResta (1998) to many-body order parameters for bulk multipole moments in 2D and 3D systems, which has been verified to give the correct phase diagrams of topological multipole insulatorsWheeler et al. (2019); Kang et al. (2019) despite its finite defects Ono et al. (2019). Owing to its real space nature, it also allows the investigation of disorder effects in these electric multipole insulatorsLi et al. (2020); Li and Wu (2020); Agarwala et al. (2020). Furthermore, the $\mathbb{Z}_{Q}$ Berry phase may be an alternative topological invariant to characterize the HOTIs and the correlated onesAraki et al. (2020). We note the latter two may encounter finite-size difficulties when implementing many-body calculations using the exact diagonalization method that is available only on very limited lattice sizesWheeler et al. (2019); Araki et al. (2020).

As we know, numerical techniques, such as QMC method, contribute significantly to the understanding of exotic properties in the strongly correlated systems. The main reason is the absence of exact solutions and controlled analytic approximations for the strongly interacting systems. Although remarkable advances have been made in studying the correlated higher-order topological phases so far, a systematical understanding is still incomplete. Indeed, we note that the higher-order topological spin phases are even largely unexplored by unbiased large-scale numerical simulations.

In this work, we employ a sign-problem-free QMC method to study the higher-order topological properties in a spin-model analogue of the Benalcazar-Bernevig-Hughes (BBH) model. We first treat the spin BBH model with Heisenberg couplings using linear spin wave theory (LSWT). While the Néel-valence band solid (VBS)vbs transition can be qualitatively described by LSWT, the critical values are greatly underestimated. Then we determine the precise transition points using finite-size scaling of dimensionless quantities characterizing the antiferromagnetism (AFM) obtained by the QMC simulations. The phase diagrams in the $(g,\Delta)$ and $(J_{1x},J_{1y})$ planes are mapped out, in which the VBS phases with $J_{1x(y)}<J_{2}$ are identified to be spin higher-order topological phase (SHOTP), which are analogues of the HOTIs in the BBH model, and have gapless corner modes in the spin excitations. The gapless spin corner modes are demonstrated by applying a perpendicular magnetic fields to open lattice. In the presence of a spin gap, only four free corner spins can be aligned by the magnetic field, generating a quantized magnetization $M_{z}=2$. Besides, the magnetization mainly distributes near the corners, confirming the gapless corner spins cause the magnetization under the magnetic field. We finally discuss the connection between the SHOTPs and the electronic correlated HOTIs.

This paper is organized as follows. Section II introduces the precise model we will investigate, along with our computational methodology. Section III presents LSWT of the spin BBH model with Heisenberg couplings. Section IV uses QMC simulations to study the bulk properties of the spin models. Section V demonstrates the spin corner states. Section VI shows the results of higher-order topological invariants. Section VII discusses the connection to the fermionic BBH-Hubbard model, and is followed by some further discussion and conclusion in Section VIII.

We consider the following dimerized spin Hamiltonian on a square lattice with four spin-$1/2$ degrees of freedom per unit cell,

where $\hat{h}_{{\bf r},i;{\bf r^{\prime}},j}=[S^{x}_{{\bf r},i}S^{x}_{{\bf r^{\prime}},j}+S^{y}_{{\bf r},i}S^{y}_{{\bf r^{\prime}},j}+\Delta S^{z}_{{\bf r},i}S^{z}_{{\bf r^{\prime}},j}]$ with $i,j=1,2,3,4$ marking the site in the unit cell, and $\Delta$ a dimensionless parameter characterizing the anisotropy of the exchange couplings. $S^{\gamma}_{{\bf r},i}$ is spin-$1/2$ operator on the site $i$ of the unit cell at ${\bf r}$, which obeys commutation relations, $\left[S^{\alpha}_{{\bf r},i},S^{\beta}_{{\bf r},j}\right]=i\hbar\varepsilon_{\alpha\beta\gamma}S^{\gamma}_{{\bf r},i}\delta_{ij}\delta_{{\bf r},{\bf r^{\prime}}}$ with $\varepsilon_{\alpha\beta\gamma}$ the Levi-Civita symbol and $\alpha,\beta,\gamma=x,y,z$ representing spin direction. $J_{1,x},J_{1,y}(J_{2,x},J_{2,y})$ are intra(inter)-cell exchange antiferromagnetic (AF) couplings. The model with $\Delta=0(1)$ corresponds to a dimerized XY(Heisenberg) spin system. By tuning the anisotropy $\Delta$, the topological property of the XXZ spin model can also be explored. Throughout the manuscript, we set $J_{2,x}=J_{2,y}=J_{2}=1$ as the energy scale.

In the following discussions, we employ the approach of stochastic series expansion (SSE) QMC method Syljuåsen and Sandvik (2002); Syljuåsen (2003) with directed loop updates to study the model in Eq.(1). The SSE method expands the partition function in power series and the trace is written as a sum of diagonal matrix elements. The directed loop updates make the simulation very efficient et al. (2011); Alet et al. (2005); Pollet et al. (2004). Our simulations are on square lattices with the total number of sites $N=4L\times L$ with $L$ the linear size. There are no approximations causing systematic errors, and the discrete configuration space can be sampled without floating point operations. The temperature is set to be low enough to obtain the ground-state properties. For such spin systems on bipartite lattice, the notorious sign problem in the QMC approach can be avoided.

Let us first get some insights on the properties of the model Eq. (1) by LSWT analysis. The Heisenberg model with $\Delta=1$ described by Eq.(1) can be treated within LSWT by replacing the spin operators by bosonic ones via Holstein-Primakoff (HP) transformationHolstein and Primakoff (1940). The transformation on sublattice $1,2$ (the spin is in the positive $z$-direction) is defined as

On sublattice $3,4$ (the spin is in the negative $z$-direction), the spin operators are defined as

Then the bosonic tight binding Hamiltonian becomes $H=\sum_{\bf k}X^{\dagger}_{\bf k}{\cal{H}}({\bf k})X_{\bf k}$, where $X_{\bf k}=\left(a_{1\mathit{\mathbf{k}}}^{\dagger}\;\;a_{2\mathit{\mathbf{k}}}^{\dagger}\;\;b_{3\mathit{\mathbf{k}}}\;\;b_{4\mathit{\mathbf{k}}}\;\right)^{T}$ is the basis, and

with

The Hamiltonian ${\cal H}({\bf k})$ in Eq.(4) can be diagnalized using Bogliubov transformationToth and Lake (2015). The spin wave spectrum contains two branches, each of which is two-fold degenerate (see Appendix A for the analytical form). The lower energy spin-wave excitations in our model display a linear spectum for small $k$, which is consistent with their AFM nature. For the $C_{4}$ symmetric case, we let $J_{1x}=J_{1y}=J_{1}$, and define the ratio $g=J_{1}/J_{2}$. The magnon spectrums at $g=0.05,0.2$ are shown in Fig.2(a) and (b). Although the spectrums are gapless, there is a gap between the two branches, and its size monotonously decreases as $J_{1}$ is strengthened. The gap vanishes at $J_{1}=J_{2}$ when the system becomes 2D AF Heisenberg model. The resulting gapless spectrum with two branches is connected to the normal one under a two-site unit cell by a folding of the Brillouin zone.

The AFM staggered order parameter defined as

is obtained in LSWT, where $i\in(1,2)$ denotes those sites which are spin up sublattice sites, and $i\in(3,4)$ denotes those sites which are spin down sublattice sites. Writing $\left\langle S_{i}^{z}\right\rangle$ in terms of HP operators, we have

where $S$ is the spin.

We then use the Bogliubov transformation to write the basis $X$ in terms of the diagonal basis $Y=\left(\alpha_{1\mathit{\mathbf{k}}}\;\;\alpha_{2\mathit{\mathbf{k}}}\;\;\beta_{3\mathit{\mathbf{k}}}^{\dagger}\;\;\beta_{4\mathit{\mathbf{k}}}^{\dagger}\;\right)^{T}$: $X=BY$. The Bogliubov transformation matrix $B$ satisfies $Bs_{z}B^{\dagger}=s_{z}$, with

At zero temperature, only the terms containing $\left\langle\beta_{3\mathit{\mathbf{k}}}\beta_{3\mathit{\mathbf{k}}}^{\dagger}\right\rangle=\left\langle\beta_{4\mathit{\mathbf{k}}}\beta_{4\mathit{\mathbf{k}}}^{\dagger}\right\rangle$ are nonzero $(=1)$. Hence the staggered order parameter is

where $(B^{\dagger}B)_{ii}$ denotes the $i$-th element of the matrix $B^{\dagger}B$. Figure 2(c) shows $m_{s}^{z}$ in Eq.(8) as a function of $g$. LSWT greatly underestimates the transition point, and it predicts the Néel-VBS transition at $g_{c}=0.11$. The periodic lattice in Fig.1 is invariant under interchanging $J_{1}$ and $J_{2}$ due to the translation symmetry, indicating a duality between $g$ and $1/g$. Hence there is also a critical value in the large $g$ region, which is exactly $1/g_{c}$ with $g_{c}$ the small critical value. Although the gap size is not affected by the duality, the topological property of the VBS phase changes to be trivial when $J_{1}$ becomes stronger.

In this section, we investigate the bulk propertires of the higher-order topological spin models using QMC simulations. The properties of the system are controlled by the ratio $J_{1x}/J_{2}$ and $J_{1y}/J_{2}$, and there exist several scenarios. To obtain some insights, let us first discuss two special limits $J_{1x}=J_{1y}=1$ and $J_{1x}=J_{1y}=0$, at which the solutions are well known or exactly solvable. In the former limit, the Hamiltonian Eq. (1) is the XXZ spin model. The case with $\Delta=1$ corresponds to an isotropic AF Heisenberg model, which has a ground state with long-range AF order. In the latter limit, the lattice is decoupled into isolated $2\times 2$ plaquettes. The ground-state energy is $E_{0}=-J_{2}(\Delta+\sqrt{\Delta^{2}+8})/2$, which is separated from the first-excited state by a gap with the size $-J_{2}-E_{0}$ (see Appendix B). While the ground state is in the total $S_{z}=0$ sector, the first-excited state has $S_{z}=\pm 1$. Hence the low-energy excitation is a spin one, which is gapped in the four-site spin chain.

For general cases $J_{1x}\neq 0$ as well as $J_{1y}\neq 0$, a Néel-VBS quantum phase transition can happen by tuning the coupling ratio $J_{1x}/J_{2}$ or $J_{1y}/J_{2}$. The phase transition points can be precisely determined using the SSE QMC by measuring the staggered magnetization $m_{s}^{z}$, the spin stiffness $\rho_{s}$, and the uniform susceptibility $\chi$. Here, the staggered magnetization $m_{s}^{z}$ is defined asWenzel et al. (2008); Sandvik (2010); Ma et al. (2018); Ran et al. (2019); Jiang (2012)

and its Binder ratio is $Q_{2}=\frac{\langle(m_{s}^{z})^{4}\rangle}{\langle(m_{s}^{z})^{2}\rangle^{2}}$Binder (1981); Binder and Landau (1984). The uniform susceptibility writes asWenzel et al. (2008); Sandvik (2010); Ma et al. (2018); Ran et al. (2019)

In terms of SSE configurations, the spin stiffness can be obtained by the expressionWenzel et al. (2008); Sandvik (2010); Ma et al. (2018); Ran et al. (2019)

where $W_{\alpha}=(N_{\alpha}^{+}-N_{\alpha}^{-})/L_{\alpha}$ is the winding number in the $\alpha$ ($x$ or $y$) direction, which takes integer values, and represents the times of spin transporting across the system; $N_{\alpha}^{+}(N_{\alpha}^{-})$ denotes the total number of operators transporting spin in the positive (negative) directionnad . Due to the Lorentz symmetry, the correlation length is the same in the space and time directions, and hence the dynamic exponent is $z=1$Troyer et al. (1997). The temperature appears in the argument $L^{z}/\beta$ of the scaling function, thus we took $\beta=L$ to exclude the temperature dependence in the finite-size scaling. At the quantum critical point, the dimensionless quantities $Q_{2},L\chi,L\rho_{s}$ are size-independent, and should cross for different lattice sizes, from which we can extract the phase transition points without knowing the critical exponents.

Let us first consider the special case with $C_{4}$ symmetry, i.e., $J_{1x}=J_{1y}$ (the value is denoted as $J_{1}$). Figure 3 shows the critical point $g_{c}$ as a function of the anisotropy $\Delta$. For $g<g_{c}$, the system is in the VBS phase, characterized by $m_{s}^{z}=\rho_{s}=\chi=0$. The AF long-range order develops above $g_{c}$, where the values of the staggered magnetization, the spin stiffness and the uniform susceptibility become finite. Since the Néel-VBS phase transition is continuous, $g_{c}$ is determined by the crossings of the above dimensionless quantities. As examples, we show $L\rho_{s}$ at $\Delta=0,1$ for different lattice sizes in Fig.4. Using the leading scaling ansatz for a second-order phase transition,

the correlation length critical exponent $\nu$ is determined by the best data collapse. As shown in the inset of Fig.3, $\nu$ increases continuously with $\Delta$. The Heisenberg model with $\Delta=1$ belongs to the three-dimensional (3D) O(3) universality, and the obtained value of $\nu$ is consistent with the standard O(3) one of $\nu=0.7112(5)$Campostrini et al. (2002); Pelissetto and Vicari (2002) For the XY model with $\Delta=0$, the universality class is the 3D XY one. The critical exponent we obtain is $\nu=0.671$, consistent with the accurate results in the literaturePelissetto and Vicari (2002).

Next we map out the phase diagram for general cases in the $(J_{1x},J_{1y})$ plane. Here we only focus on the Heisenberg and XY exchange couplings, and the phase diagrams for $0<\Delta<1$ are similar. The critical points are determined more precisely using the best data collapses with the accurate $\nu$ for the O(3) and XY universality classes. The phase diagrams in Fig. 5 consist of the AFM and VBS phases. The region of VBS shrinks as $\Delta$ is decreased. In both figures, the VBS phases with $J_{1x(y)}/J_{2}<1$ are SHOTPs, which will be demonstrated in the following section. Here the SHOTPs survive only in part of the unit square. In contrast, the electronic HOTIs of the BBH model are in the entire square region $|t_{1x(y)}/t_{2}|<1$ [ $t_{1x(y)}(t_{2})$ denotes the hopping amplitude of the weak (strong) bonds].

Then we consider another limit $J_{1x}=0$ or $J_{1y}=0$, at which the system becomes one-dimensional dimerized spin chain on the boundaries. Under the Jordan-Wigner transformationColeman (2015)

the following 1D fermionic Hamiltonian is obtained,

where $\alpha_{i},\beta_{i}$ are annihilation operators of spinless fermions on the two sites of unit cell $i$, and $n_{i,\alpha}=\alpha^{\dagger}_{i}\alpha_{i},n_{i,\beta}=\beta^{\dagger}_{i}\beta_{i}$ are the corresponding number operators. In the hopping term connecting the boundary, there is an additional phase $e^{i\pi n_{t}}$ ($n_{t}$ the total number of fermions), which may be $1$ or $-1$ depending on even or odd $n_{t}$. Since such a sign has no effect on the 1D boundary mode, it does not affect the topological property of the system, and we omit it in the above Hamiltonian. For the XY model, the mapped fermionic Hamiltonian is the noninteracting Su-Schrieffer-Heeger (SSH) modelSu et al. (1979), which has a topological phase transition at $J_{1x(y)}=J_{2}$. The 1D dimerized Heisenberg model corresponds to an interacting SSH model with dimerized nearest-neighbor repulsions. It has been shown that an AFM transition occurs at $J_{1x(y)}=J_{2}$Cazalilla et al. (2011); Mishra et al. (2011). We calculate the 1D topological invariant, i.e., the Zak phase, which is $\gamma=1$ for $J_{1x(y)}<J_{2}$ and $\gamma=0$ for $J_{1x(y)}>J_{2}$ (see Appendix C). This unambiguously shows the topological phase transition happens exactly at the uniform point, which determines the on-axis points of the phase boundary in Fig.5.

It is noted that the transition line between VBS I and VBS II is a straight line. In this region, the dimensionless quantities, such as $L\rho_{s}$ $L\chi$, do not exhibit universal crossings as $J_{1x}(J_{1y})$ varies at fixed $J_{1y}(J_{1x})$, and instead their values continuously decrease with increasing $L$. These imply the absence of a Néel-VBS transition, and the system keeps in the VBS phase. Since there is a duality between $J_{1x(y)}<1$ and $J_{1x(y)}>1$, the transition point should be exactly at $J_{1x(y)}^{c}=1$. Otherwise it will generate a contradictory result. Let us explain this point more clear. Suppose $J_{1x}^{c}<1$ for fixed $J_{1y}$ in the straight-line region. Then there are two topological transition points due to the duality, and an intermediate phase exists in the region $J_{1x}^{c}<J_{1x}<1/J_{1x}^{c}$. Since the VBS with $J_{1x}<J_{1x}^{c}(J_{1x}>1/J_{1x}^{c})$ is adiabatically connected to the $J_{1x}=0(J_{1x}=\infty)$ topologically nontrivial (trivial) phase, the intermediate phase has contradictory topological properties from the topological transitions of the two different sides, suggesting it should not exist in real system and the transition point is exactly at $J_{1x}^{c}=1$.

In the following sections, we will show that the VBS for $J_{1x}<1$ or $J_{1y}<1$, which is denoted as VBS I, has nontrivial higher-order topological property, while the other one (denoted as VBS II) is topologically trivial.

In the above phase diagrams, the VBS phases are gapped in the spin excitation while the AF state is gapless. To see this explicitly, we add a magnetic field perpendicular to the lattice described by

The total Hamiltonian $\hat{H}_{0}+\hat{H}_{1}$ is then simulated by the QMC method. The magnetization density

is further measured. We find $m_{z}$ and $\rho_{s}$ remain to be zero up to a finite magnetic field. Afterwards, they continuously increase, and the system is in a ferromagnetic phase. To determine the transition values precisely, we fit the data with $\rho_{s}\sim[h-h_{c}(L)]^{\beta}$ around the critical point for each lattice sizeKrzakala et al. (2008); Capogrosso-Sansone et al. (2008). As shown in Fig.6(b), our data are compatible with the mean-field exponent $\beta=\frac{1}{2}$. Then $h_{c}$ is obtained by extrapolating $h_{c}(L)$ to the thermodynamic limit [see Fig.6(c)]. By this way, the phase diagram in the $(g,h_{c})$ plane is plotted in Fig.6(a). Below $h_{c}$, the system remains in the VBS phase with gapped spin excitations. However the spins are aligned along the direction of the applied magnetic field above $h_{c}$, displaying a ferromagnetic behavior with gapless spin excitations. The critical value is exactly $h_{c}=1$ at $g=0$, when the system is decoupled into isolated plaquettes, and can be understood analytically. The spin dimer can be in a spin-singlet ( or spin-triplet) state, with the energy $-\frac{3}{4}J_{2}$ (or $\frac{1}{4}J_{2}-h$). So the transition point is $h=J_{2}=1$, after which the ground state becomes the spin-triplet one where the spins are aligned to the same direction by the magnetic field.

It is interesting to note that the spin-gapped VBS phases are higher-order topological spin states. Next we will demonstrate that gapless spin corner modes will appear on open lattices in both $x$- and $y$-directions due to nontrivial higher-order bulk topology.

In this section, we demonstrate the spin corner modes, which consists of the characteristic signatures of nontrivial SHOTPs. Here we consider a square geometry with open boundary condition, which preserves the mirror symmetry $M_{x},M_{y}$, the inversion symmetry, and the rotation symmetry $C_{2}$. For other geometries, a general principle is that a corner state will appear as long as a domain wall is formed by the two edges crossing at the corner. The SHOTP within $J_{1x(y)}/J_{2}<1$ is adiabatically connected to the limit $J_{1x(y)}=0$, where there are four dangling spins at the corners. For finite $J_{1x(y)}$, the gapless spin corner modes should remain on open lattices in the SHOTP. A infinitesimal magnetic field $h$ can align these corner spins, thus the total magnetization becomes $M_{z}=Nm_{z}=2$, compared to $M_{z}=0$ in the bulk SHOTP. Besides, the magnetization is mainly localized near the corners, and the localization length increases with increasing $J_{1x(y)}$ due to the decreasing of the spin gap. Figure 7 shows the total magnetization as a function of $h$ at $J_{1}/J_{2}=0.1$ (the $C_{4}$ symmetric case) on a $L=40$ open lattice. It is noted that the origin of the $M_{z}=2$ plateau tends to be zero as $\beta$ is increased, implying the plateau is quantized as long as $h$ becomes nonzero in the limit of zero temperature. Beside, the $M_{z}=2$ plateau persists to a finite $h$, and only begins to vanish when the spin gap closes. We also show the local distribution of the magnetization for $J_{1}/J_{2}=0.1$ in Fig.7(b). Indeed the magnetization is localized near the corners, implying the gapless corner spins leads to the quantized magnetization.

Such results are also valid for larger $J_{1}/J_{2}$ as long as the ratio is less than the critical value $g_{c}=0.548$ (see Appendix D). Since the appearance of the spin corner states is a manifestation of the nontrivial higher-order topological property, the Néel-VBS transition is also a higher-order topological phase transition. Besides, we perform the calculations of $M_{z}-h$ curves with very small steps near the straight-line boundary of the phase diagram in Fig.5(b). The $M_{z}$ plateau characterizing the spin corner states immediately vanishes as $J_{1x}$ crosses the boundary located at $J^{c}_{1x}=1$, thus verifies the straight phase boundary with very high accuracy (see Appendix D).

The appearance of gapless spin corner states is due to the bulk higher-order topological protection, and the bulk topology is generally characterized by a topological invariant. A $\mathbb{Z}_{4}$ Berry phase is proposed as the topological invariant for the $C_{4}$-symmetric SHOTP. It is defined with respect to the local twist of the Hamiltonian, which is parameterized by three independent variables $\Theta=(\theta_{1},\theta_{2},\theta_{3})$ ranged in $[0,2\pi)$Araki et al. (2020). The Hamiltonian is decomposed into two types of plaquettes with weak and strong bonds, respectively. The twist parameters are introduced to the spin operators on any one of the plaquettes as: $\tilde{S}^{+}_{i}=e^{-i\varphi_{i}}S^{+}_{i},\tilde{S}^{-}_{i}=e^{i\varphi_{i}}S^{-}_{i}$, with $\varphi_{i}=\sum_{q=1}^{i}\theta_{i}$ for $i=1,2,3,4$ and $\varphi_{4}=0$. The $\mathbb{Z}_{4}$ Berry phase is defined as

where the path $L_{j}(j=1,2,3,4)$ is in the parameter space: $E_{j-1}\rightarrow G\rightarrow E_{j}$ with $E_{1}=(2\pi,0,0),E_{2}=(0,2\pi,0),E_{3}=(0,0,2\pi),E_{0}=E_{4}=(0,0,0)$, and $G=1/4\sum_{j=1}^{4}E_{j}$, and $\Psi_{\Theta}$ is the corresponding groundstate many-body wave function at half-filling. In our calculations, we choose the left-bottom plaquette to apply the local twist (see Fig.1). Figure 8 shows the $\mathbb{Z}_{4}$ Berry phase in the $(g,\Delta)$ plane. As $g$ increases at fixed $\Delta$, the value of the $\mathbb{Z}_{4}$ Berry phase changes abruptly from $\pi$ to $0$ at a critical $g_{c}$, marking a topological phase transition. The critical values are $g_{c}(\Delta=1)=1$ and $g_{c}(\Delta=1)=1.19$, respectively. For intermediate $\Delta$, the critical value continuously decreases with increasing $\Delta$. Compared to the QMC results, the $\mathbb{Z}_{4}$ Berry phase greatly overestimates the transition point, and the quantized value persists in the AF phase, which has also been found in other works, and has been attributed to the finite-size errorBibo et al. (2020). Since the Néel-VBS transition is continuous, the two phases coexist near the transition point on small finite lattices. Besides, a finite-size gap exists even in the AF phase. Hence although the $\mathbb{Z}_{4}$ Berry phase is quantized, it still fails to determine the phase boundary precisely. Another method of calculating the topological invariant in terms of the many-body quadrupole momentum also suffers significant finite-size errorsWheeler et al. (2019); Kang et al. (2019) (see Appendix E), and identifying a solid many-body topological invariant requires further studiesYou et al. (2020).

Here we further compare SHTOPs discussed in previous sections with correlated HOTIs for a better understanding of their properties. The Heisenberg model with $\Delta=1$ in Eq.(1) is the large-$U$ limit of the BBH electronic model subjected to an on-site Hubbard interaction $H_{U}=U\sum_{i}n_{i\uparrow}n_{i\downarrow}$, with $n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ and $c^{\dagger}_{i\sigma},c_{i\sigma}$ the electronic creation and annihilation operators with spin $\sigma=\uparrow,\downarrow$. At $t_{1}/t_{2}=1$, the BBH model becomes the $\pi$-flux lattice hosting 2D Dirac fermionsWang et al. (2014); Li et al. (2015); Parisen Toldin et al. (2015); Otsuka et al. (2016); Guo et al. (2018). The Hubbard interaction drives a semimetal-AFM quantum phase transition, with the most accurate value of the critical point $U_{c}=5.65\pm 0.05$Parisen Toldin et al. (2015); Otsuka et al. (2016). In the large-$U$ limit, the relation $J=4t^{2}/U$ gives the exchange constant in terms of $U$ and the hopping amplitude. Since the VBS-AFM transition happens at $g_{c}=0.584$, there is a critical ratio for $t_{1}/t_{2}$ at $f_{c}=\sqrt{g_{c}}=0.764$. Below $f_{c}$, the correlated HOTI continuously evolves into the SHOTP as $U$ increases, and meanwhile the spin excitation remains gapped in the bulk systems. However above the critical point $f_{c}$, there is a AFM transition, from which the bulk spin excitation becomes gapless, and hence the correlated HOTI vanishes. The above discussion is summarized as the phase diagram in the $(t_{1}/t_{2},U)$ plane, shown in Fig.9(a).

In the correlated HOTI, the spin excitation is gapless, but the charge excitation is gapped on open lattices. Figure 9(b) and (c) shows these gaps as a function of $U$ in the limit $t_{1}=0$. At half filling, each of the four isolated corner sites is occupied by one electron. Adding an opposite spin electron to any such site unavoidably induces the interacting energy $U$, and thus the charge excitation gap is $\Delta_{c}=U$. However the spin excitation by flipping the corner spin does not cost any energy, and the spin excitation gap is $\Delta_{s}=0$. We also calculated these gaps in the bulk system, which are defined as follows: $\Delta^{bulk}_{c}=(E_{N+1,S_{z}=1/2}+E_{N-1,S_{z}=-1/2}-2E_{N,S_{z}=0})/2$, and $\Delta^{bulk}_{s}=(E_{N,S_{z}=1}+E_{N,S_{z}=-1}-2E_{N,S_{z}=0})/2$, where $E_{N,S_{z}}$ is the ground-state energy with the total number of electrons $N$ and the total spin in the $z$-direction $S_{z}$. As shown in Fig.9(b) and (c), both of them are gapped. So the correlated HOTIs have similar bulk-boundary correspondence as the spin counterparts, i.e., the appearance of gapless spin corner states protected by the higher-order topology.

In summary, we investigate the higher-order topological properties of spin analogues of the BBH model using large-scale QMC simulations. Generally two spin quantum phases, i.e., VBS and AFM, may appear in the system. A continuous phase transition can be driven by tuning the ratio between the weak and strong couplings. We identify the VBS phase with $J_{1x(y)}/J_{2}<1$ to be a nontrivial SHOTP, which is characterized by a bulk spin excitation gap and gapless spin corner states on open lattice. Besides, in the large-$U$ limit, the SHOTP is continuously connected to the electronic correlated HOTI in the BBH-Hubbard model. They share the same bulk-boundary correspondence: the appearance of gapless corner spin excitations protected by nontrivial higher-order topology. The SHOTPs also exist in XXZ spin models, but their regions in the phase diagrams shrink as the anisotropy $\Delta$ decreases. Our results explicitly demonstrate the higher-order topological properties of spin systems, and contribute to further understandings of the many-body higher-order topological phenomena.

The authors thank Xiancong Lu, Zhixiong Li, Tianhe Li, Hiromu Araki, Julian Bibo, Rongqiang He, Hua Jiang, Nvsen Ma for helpful discussions. H.G. acknowledges support from the NSFC grant Nos. 11774019 and 12074022, the Fundamental Research Funds for the Central Universities and the HPC resources at Beihang University. X.Z. and S.F. are supported by the National Key Research and Development Program of China under Grant No. 2016YFA0300304, and NSFC under Grant Nos. 11974051 and 11734002.