Correlation Effects in Quadrupole Insulators: a Quantum Monte Carlo Study

TL;DR

This study uses large-scale quantum Monte Carlo simulations to explore how interactions affect the topological properties of quadrupole insulators, revealing a transition to an antiferromagnetic phase with unique critical behavior.

Contribution

First application of Green's function formalism in interacting quadrupole insulators, demonstrating stability against weak interactions and identifying a new universality class for the phase transition.

Findings

Quadrupole insulator remains topologically stable under weak interactions.

Strong interactions induce a continuous transition to an antiferromagnetic insulator.

Critical exponents suggest a new universality class distinct from known AFM transitions.

Abstract

The quadrupole insulator, a high-order topological insulator, with on-site Hubbard interaction is numerically studied by large-scale projector quantum Monte Carlo (PQMC) simulations. The Green's function formalism is successfully used to characterize topological properties in interacting quadrupole insulators for the first time. We find that the topological quadrupole insulator is stable against weak interactions and turns into a trivial antiferromagnetic (AFM) insulator by a continuous topological phase transition (TPT) for strong interactions. The critical exponents related to the TPT are estimated to be $\nu=0.67(4)$, $\beta=0.40(2)$, which are distinct from those of the known AFM transitions and suggest a new universality class.

Figures

Equations

$$H_0=$$

$$H_I=$$

$$H_0=$$

$$H_I=$$

$$⟨\hat{O}⟩=\frac{⟨{\Psi }_T|e^{-\frac{\Theta }{2}H}\hat{O}{e}^{-\frac{\Theta }{2}H}|{\Psi }_T⟩}{⟨{\Psi }_T|e^{-\Theta H}|{\Psi }_T⟩}.$$

$$⟨{\Psi }_T|e^{-\Theta H}|{\Psi }_T⟩=⟨{\Psi }_T|{\left(e^{-{\Delta }_{\tau }{H}_0}{e}^{-{\Delta }_{\tau }{H}_I}\right)}^M|{\Psi }_T⟩+O\left({\Delta }_{\tau }^2\right).$$

$$e^{-\frac{{\Delta }_{\tau }U}{2}{\left(n_{i\uparrow }+n_{i\downarrow }-1\right)}^2}=\frac{1}{4}\sum _{s_i}\gamma \left(s_i\right)e^{\sqrt{-{\Delta }_{\tau }U}\eta \left(s_i\right)\left(n_{i\uparrow }+n_{i\downarrow }-1\right)},$$

$$H=\sum _{n,k}{\gamma }_{n,k}^{+}{ϵ}_{n,k}{\gamma }_{n,k},{\gamma }_{n,k}=\sum _{\alpha }{\left[u_k^{\ast n}\right]}^{\alpha }{c}_{k,\alpha },$$

$${\mathcal{W}}_x\left(\mathbf{k}\right)=G_{\mathbf{k}+\left(N_x-1\right)\Delta k_x}...G_{\mathbf{k}+\Delta k_x}{G}_{\mathbf{k}}.$$

$$|{\Psi }_R^j(k_y)⟩=\frac{1}{\sqrt{N_x}}\sum _{n=1}^{N_{\mathrm{occ}}}\sum _{k_x}{\left[{\nu }_{x,k_x}^j(k_y)\right]}^n{e}^{-i{k}_{x}R}{\gamma }_{n,k_x,k_y}^{+}|0⟩,$$

$${\mathcal{W}}_{y,\mathbf{k}}^{{\nu }_x^{+}}=$$

$$|w_{x,\mathbf{k}}^{+}⟩=\sum _{n=1}^{N_{occ}}|u_{\mathbf{k}}^n⟩{\left[{\nu }_{x,\mathbf{k}}^{+}\right]}^n.$$

$${\mathcal{W}}_{x,k_x}|{\nu }_{x,k_x}⟩=\exp \left[i2\pi {\nu }_x\right]|{\nu }_{x,k_x}⟩,{\nu }_x\in [-\frac{1}{2},+\frac{1}{2}].$$

$$\begin{array}{rl}{M}_x:p_x^{{\nu }_y^{\pm }}& =0\mathrm{or}1/2,\\ M_y:p_y^{{\nu }_x^{\pm }}& =0\mathrm{or}1/2.\end{array}$$

$$\mathcal{I}:p_x^{{\nu }_y^{+}}=-p_x^{{\nu }_y^{-}},p_y^{{\nu }_x^{+}}=-p_y^{{\nu }_x^{-}}.$$

Full text

††preprint: APS/123-QED

The bulk-edge correspondence Ryu and Hatsugai (2002) is one of the most fundamental concepts for understanding topological insulators (TIs) in $d$-dimensional systems which possess gapless modes in their ($d-1$)-dimensional boundaries Fu et al. (2007); Kane and Mele (2005); Fu and Kane (2007). Recently, topological crystalline insulators generalize this bulk-edge correspondence to high-order topological insulators (HOTIs) Benalcazar et al. (2017, 2018); Song et al. (2017a); Liu et al. (2019a). In addition, the concept of HOTIs is also generalized to high-order topological superconductors Bultinck et al. (2019); Wang et al. (2018); Yan et al. (2018); Langbehn et al. (2017a, b). The predicted Majorana corner/hinge states in such systems are promising building blocks (qubits) for realizing quantum computation Bomantara and Gong ; Pahomi et al. .

Attracted by the novel properties and potential applications of HOTIs, many theoretical works have been done to predict HOTIs in realistic systems, such as quasicrystals, sonic crystals, non-hermitian systems, and solid state systems Schindler et al. (2018a); Robredo et al. (2019); Olsen et al. (2019); Song et al. (2017b); Chen et al. (a); Lin et al. ; Yang et al. ; Liu et al. (2019b). And some experimental works have even shown the direct observation of HOTIs in photonic/phononic crystals, electrical circuits as well as crystalline solids, such as bismuth Serra-garcia et al. (2018); Chen et al. (b); Peterson et al. (2018); Serra-Garcia et al. (2019); Ni et al. (2019); Schindler et al. (2018b).

Unlike bosonic systems, fermionic systems with spin degrees of freedom are subjected to the Coulumb repulsion which may give rise to many exotic phenomena Fang et al. (2013); Zhang and Liu (2015); Jiang et al. (2018); Kivelson and Heim (1982); Guo and Shen (2011); He et al. (2016a, 2018); Hohenadler et al. (2011). However, there are only a few studies which discuss correlation effects in HOTIs Rasmussen and Lu ; Peng et al. (2017); You et al. (2018). Unbiased numerical research on interacting HOTIs is almost a blank and is the main motivation of this work.

In this paper, we focus on such a family of HOTIs, which are called quadrupole topological insulators. The model is constructed in a non-interacting spinless system with only nearest neighbor hopping terms in two dimensions Benalcazar et al. (2018, 2017). In conventional $d$-dimensional TIs, $(d-1)$-dimensional edge modes can be understood by electric dipoles and charge pumps, which are related to the Wilson loops and energy bands. However, the bulk topology of multipole insulators need to be understood by electric multipole moments and dipole pumps, which are respectively related to the nested Wilson loops and Wannier bands, respectively. The two-dimensional quadrupole insulator studied here has a vanishing dipole moment and a quantized quadrupole moment in its bulk. Then, the bulk quadrupole moment in a finite sample gives rise to dipole moments on its one-dimensional edges and zero-dimensional corner charges. In order to discuss correlation effects in quadrupole insulators, we introduce spin degrees of freedom and on-site Hubbard interactions to quadrupole insulators and simulate the system using sign-problem-free projector quantum Monte Carlo (PQMC) method.

Model and method.—Together with a concrete spinless lattice model with just nearest neighbor hopping terms defined in a square lattice, the concept of electric multipole moments in solid state systems is first established in Ref.Benalcazar et al. (2018). In order to study the correlation effects in quadrupole insulators, we stack up two copies of this model and introduce spin degrees of freedom ($\uparrow$ and $\downarrow$) to each copy, respectively. Mirror symmetries $M_{x}$ and $M_{y}$, which protect topological phases, still survive, since we do not introduce spin-orbital coupling terms. So, it is enough to characterize this spinful system by simply discussing topological properties carried by only one spin flavor. After adding the on-site Hubbard interaction, which accounts simply for the Coulumb repulsion, the model studied here is described by $H=H_{0}+H_{I}$:

[TABLE]

Wannier-sector polarizations $p_{x}^{\nu_{y}^{\pm}}(p_{y}^{\nu_{x}^{\pm}})$, which describe the electronic polarizations along $x(y)$ within Wannier bands $\nu_{y}^{\pm}(\nu_{x}^{\pm})$, are quantized to $1/2$ or [math] due to the mirror symmetries $M_{x}$ and $M_{y}$, where $\nu_{x}^{\pm}$($\nu_{y}^{\pm}$) stands for the upper ($+$) or lower ($-$) Wannier band along $k_{x}$($k_{y}$). (More details in Supplemental Material.) Then, different topological phases can be distinguished by the topological index ${\bf{P}}^{\nu}=(p_{x}^{\nu_{y}^{-}},p_{y}^{\nu_{x}^{-}})$ (Fig.1(b)). Wannier values and the locations of Wannier gap closing point at each phase boundary are also given in Fig.1(b).

Setting parameters $\lambda_{x}=\lambda_{y}=\lambda$ as energy unit, we investigate the model introduced above at half filling using large-scale PQMC simulation, which is sign-problem free as long as $H_{0}$ includes only nearest-neighbour hopping terms Assaad and Evertz (2008). Another issue we need to address is to distinguish different topological phases in a correlated system. In the conventional TIs, it is the Green’s function formalism that plays a key role in identifying different topological phases after introducing interactions Wang and Zhang (2012a); Wang and Yan (2013); Wang et al. (2010); Wang and Zhang (2012b); Trifunovic and Brouwer (2019). In this paper, We simply use the zero-frequency Green’s function $G^{-1}({\bf{k}},i\omega=0)$ to replace the non-interacting bloch Hamiltonian $H_{0}(\bf{k})$ in the calculation of the topological index $\bf{P}^{\nu}$. (cf. Supplemental Material.)

Single-particle level TPTs with interactions.—In order to investigate correlation effects in quadrupole insulators, we need to verify whether the zero-frequency Green’s function mentioned above can characterize the topological properties in quadrupole insulators. Unlike in the conventional topological insulators, topological properties in quadrupole insulators studied here are encoded in the Wannier bands. In other words, the topological phase is protected by the mirror symmetries and can not transform to another phase as long as the Wannier gaps, rather than the energy bands for the conventional topological insulators, keep open. So, we first investigate whether TPTs can still be detected by the Wannier gap closing events after introducing on-site Hubbard interaction as in the non-interacting case.

The TPTs at the single-particle level can be driven by a pair of hopping parameters $(\gamma_{x},\gamma_{y})/\lambda$ as shown in Fig.1(b). We choose the path in parameter space, $(0.75,0.75)\to(0.75,1)\to(0.75,1.25)\to(1,1.25)\to(1.25,1.25)$, which is discussed in Ref.Benalcazar et al. (2017). Starting from the quadrupole phase characterized by ${\bf{P}^{\nu}}=(1/2,1/2)$, we increase $\gamma_{y}/\lambda$ from 0.75 to 1.25 while keeping $\gamma_{x}/\lambda$ fixed at 0.75. As the quadrupole phase approaches the TPT, the Wannier bands $\nu_{x}(k_{y})$ become gapless at $k_{y}=\pi$ as they become twofold degenerate at $\nu_{x}(\pi)=1/2$. After this TPT, the Wannier gap reopens, and the ground state is still in a non-trivial topological phase (with ${\bf{P}^{\nu}}=(1/2,0)$) but not a quadrupole phase. Then with keeping $\gamma_{y}$ fixed at $\gamma_{y}/\lambda=1.25$ and increasing $\gamma_{x}$, the gap of another Wannier band $\nu_{y}(k_{x})$ is closing at $k_{x}=\pi$ at the transition point $(1,1.25)$. Finally, the electric polarization $p_{x}^{\nu_{y}^{-}}$ within Wannier band $\nu_{y}^{-}$ changes from $1/2$ to [math], which means the ground state is in a completely trivial phase with ${\bf{P}^{\nu}}=(0,0)$.

The evolution of Wannier bands through the path mentioned above are calculated by the PQMC using the Green’s function formalism. Numerical results are shown by square symbols in Fig.2 with system size $L=6$ and interaction strength $U=2\lambda$. We find that the Wannier bands $\nu_{x}$ and $\nu_{y}$ experience a gap closing and reopening process through the TPT similar to the non-interacting case. To make a concrete contrast, we also calculate Wannier bands and Wannier-sector polarizations with $U=0$ in the same path. Numerical results are shown in Fig.2 by pure lines. Since system the interacting case with $U=2\lambda$ is not far from the non-interacting case, we can find that the evolution of Wannier bands in the interacting system through the path in parameter space only slightly deviate from that of the non-interacting case. What’s more, we find that the Wannier gaps calculated by $G({\bf{k}},i\omega=0)$ with finite $U$ also experience a gap closing and reopening process through the phase boundaries. These results show that the TPTs happen at $(0.75,1.0)$ and $(1,1.25)$ and also mean that TPTs can be successfully detected in interacting quadrupole insulators by our method.

To further confirm the validity of the Green’s function formalism, we calculate topological index $\bf{P}$ by definition as a direct indication of TPTs (Fig.3). Still in the system with $U=2\lambda$, we find Wannier-sector polarizations $p_{x}^{\nu^{-}}$ and $p_{y}^{\nu^{+}}$ jump at $(1,1.25)$ and $(0.75,1)$ respectively, which is consistent with the results we derived from the Wannier bands argument. Thus, the topological properties in interacting quadrupole insulators can be correctly captured by $G({\bf{k}},i\omega=0)$.

In order to investigate how the phase boundaries vary with interactions, we increase the interaction strength up to $U=8\lambda$. We first investigate the transition between ${\bf{P}}^{\nu}=(1/2,1/2)$ and ${\bf{P}}^{\nu}=(1/2,0)$ by fixing $\gamma_{x}=0.8\lambda$ and increasing $\gamma_{y}$ from $0.9\lambda$ to $1.1\lambda$. Wannier bands $\nu_{x}$ close the gap at $k_{y}=\pi$ and Wannier-sector polarization $p_{y}^{\nu_{x}^{-}}$ jumps from $1/2$ to [math] at the same point $\gamma_{y}=\lambda$ as in the non-interacting case (Figs.4(a,b)). Similar results are also shown in the phase transition from ${\bf{P}}^{\nu}=(1/2,0)$ to ${\bf{P}}^{\nu}=(0,0)$ with increasing $\gamma_{x}$ and fixing $\gamma_{y}=1.2\lambda$ (Figs.4(c,d)). Unlike in the conventional topological insulators Hung et al. (2013, 2014), the phase transitions between different topological phases of Wannier bands are not affected by the Hubbard interaction. The stability of phase boundaries may be accounted for by the novel definition of quadrupole insulators that the topological properties are encoded in Wannier bands rather than energy bands.

Interaction induced TPTs.—Although we have studied correlation effects in quadrupole insulators, TPTs discussed above still stay at the single-particle level as the phase transition is driven just by tight-binding parameters. A direct route to construct an interaction driven TPT in our model is increasing the on-site Hubbard $U$ until an AFM long-range order is formed. Here, we study the transition from the quadrupole phase to AFM insulator (AFMI) and find that they are continuous. The Wannier-sector polarization $p_{x}^{\nu^{-}}$ is calculated and our results clearly reveal that the bulk quadrupole moment is gradually destroyed by the AFM order.

In order to characterize the quadrupole insulator to AFMI transition, we calculate the correlation ratio of AFM order parameter: ${R_{z}}\left({U,L}\right)=1-\frac{{{C_{z}}\left({{\bf{\Gamma}}+\delta{\bf{q}}}\right)}}{{{C_{z}}\left({\bf{\Gamma}}\right)}}$ where ${C_{z}}\left({\bf{\Gamma}}\right)$ is the structure factor ${C_{z}}\left({\bf{\Gamma}}\right)=\frac{1}{{4N^{2}}}\sum\nolimits_{ij\gamma}{\left\langle{S_{i\gamma}^{z}S_{j\gamma}^{z}}\right\rangle}$ where $i,j=1,2,...,N$ run over all unit cells and $\gamma$ stands for sub-lattice index in a unit cell. We choose parameters $(\gamma_{x},\gamma_{y})/\lambda=(0.8,0.8)$, which guarantee the ground state is a quadrupole insulator. Correlation ratios $R_{z}(U,L)$ with different system sizes cross at nearly the same point, which implies the transition is continuous since $R_{z}(U,L)$ is a renormalization invariant quantity in a continuous phase transition (Fig.5(a)). Further, we make a data collapse to find the critical exponents by assuming $C_{z}({\bf{\Gamma}},U,L)L^{2\beta/\nu}=f_{C}((U/U_{c}-1)L^{1/\nu})$. The data collapse process is to find $\nu$ and $\beta$ to make all data points collapse at one single unknown curve described by function $f_{C}$. Results with optimized exponents are shown in Fig.5(b). Critical exponents $\nu=0.67(4)$, $\beta=0.40(2)$ and critical point $U_{c}=13.15(2)\lambda$ is obtained, which reveals the TPT is actually continuous (Details to determine those values are given in Supplemental Material). These exponents should be compared with those of the 3D-Heisenberg universality class, the 2D dimerized Heisenberg model, and the deconfined quantum criticality in VBS-AFM transitions (found in $J_{1}$-$J_{2}$ and $J$-$Q$ models) Campostrini et al. (2002); Wenzel et al. (2008); Sandvik (2007); Gong et al. (2014). The critical exponents obtained here are inconsistent with those in these three cases, and suggest a new universality class for TPTs between quadrupole insulators and and AFM insulators.

Having known that the AFM phase is formed around $U_{c}\sim 13.15\lambda$, the interplay between the bulk topology and AFM order is still uncertain. The physics underlying the TPT here can be understood by spontaneous symmetry breaking. Different topological phases described by index ${\bf{P}}^{\nu}$ are well-defined only when preserving mirror symmetries $M_{x},M_{y}$ Benalcazar et al. (2018). In other words, Wannier-sector polarizations $p_{x}^{\nu_{y}^{\pm}},p_{y}^{\nu_{x}^{\pm}}$ are no longer quantized to $1/2$ or [math] when mirror symmetry is spontaneously broken by AFM order. Meanwhile, the Wannier gap closing mechanism, which is used to describe TPT, also fails since the topological properties encoded in different Wannier sectors $\nu^{+}$ and $\nu^{-}$ can mess up without closing the Wannier band gaps. However, in previous studies, topological index failed to describe the destruction of topology by long-range order in numerical simulations, as spontaneous symmetry breaking did not happen in a finite size system or due to the breakdown of the Green’s function formalism Lang et al. (2013); He et al. (2016b).

Hence, in order to visualize the competition between topology and AFM order directly, we simulate the system with manually introducing a very small staggered magnetic field, which can break the mirror symmetries. Specifically, we add positive or negative magnetic fields with strength $\mu=0.001\lambda$ to our system corresponding to the sites marked with red or blue in Fig.1(a). Since the AFM order plays the same role on $p_{x}^{\nu_{y}}$, $p_{y}^{\nu_{x}}$ and the constraint due to the inversion symmetry $\mathcal{I}$ are still valid Benalcazar et al. (2018) (details in Supplemental Material), we only use one of them to characterize the topology in the quadrupole phase in Fig.6(a). Far from the discontinuous behavior in TPTs at the single-particle level, electric polarization $p_{x}^{\nu_{y}^{-}}$ shows a smooth decrease from 1/2 to 0, which means the topological order is suppressed by the AFM order gradually as we expect. By extrapolating structure factor $C_{z}(\Gamma)$ to the thermodynamic limit with different sizes $L=4,6,8,10,12$, we find that the AFM order appears around $U=9\lambda$, which is smaller than $U_{c}\sim 13.15\lambda$, since spins are polarized toward $z$-direction under the small external staggered magnetic field $\mu$ (Fig.6(b)).

From directly monitoring the topological index around the TPT with introducing a tiny staggered magnetic field, we can clearly find the competition between AFM and topological order. The suppression of non-trivial topology happens simultaneously with the formation of AFM order, which means the AFM phase obtained here is a simple AFMI.

Summary.—We apply the Green’s function formalism to study interacting topological quadrupole insulators numerically for the first time. Through our numerical results using PQMC, we find that this formalism can successfully capture all the topological properties encoded in Wannier bands even when interactions are introduced. For the single-particle level TPTs, due to the novel machenism of TPTs in quadrupole insulators, the phase boundaries between different topological phases marked by ${\bf{P}}^{\nu}$ keep unchanged with increasing the strength of on-site Hubbard interaction up to $U=8\lambda$.

Besides the single-particle level TPTs, which have non-interacting counterpart, we also study a truly interaction driven TPT quantitatively. A continuous phase transition is found between the quadrupole insulator and AFMI. Critical exponents $\nu,\beta$ are obtained by data collapse, which suggest a new universality class for HOTIs. By adding a tiny staggered magnetic field, the competition between topology and AFM order is revealed through the smooth decreasing of topological index $p_{x}^{\nu_{y}^{-}}$ from 1/2 to 0. This indicates that the topological order is destroyed by AFM simultaneously and the AFM phase is a trivial Mott insulator.

This work was supported by the National Natural Science Foundation of China (Grants No. 11774422 and No. 11874421). R. Q. H. was supported by the Fundamental Research Funds for the Central Universities, and the Research Funds of Renmin University of China (Grant No. 18XNLG11). Computational resources were provided by the Physical Laboratory of High Performance Computing at Renmin University of China and National Supercomputer Center in Guangzhou with Tianhe-2 Supercomputer.

S1 solving the model with projector quantum Monte Carlo method

Our model writes as $H=H_{0}+H_{I}$:

[TABLE]

Using the projection operator we can obtain the wave function of ground state $\left|{{\Psi_{0}}}\right\rangle={e^{-\frac{\Theta}{2}H}}\left|{{\Psi_{T}}}\right\rangle$ and we choose as the trial wave function $\left|{{\Psi_{T}}}\right\rangle$ the Slater determinant describing the many-body ground state of $H_{0}$. The physical observables are measured by:

[TABLE]

To treat the interaction part, we divide the projection time into $M$ slices ($\Theta=M\Delta_{\tau}$) and make the Trotter decomposition to separate interaction term $H_{I}$ from $H$ in exponential:

[TABLE]

And after the Hubbard-Stratonovich (HS) transformation:

[TABLE]

with $\gamma\left({\pm 1}\right)=1+\sqrt{6}/3,\gamma\left({\pm 2}\right)=1-\sqrt{6}/3,\eta\left({\pm 1}\right)=\pm\sqrt{2(3-\sqrt{6})},\eta\left({\pm 2}\right)=\pm\sqrt{2(3+\sqrt{6})}$, we can finally obtain the interaction terms in the formula with only fermion bilinears coupled to auxiliary fields. After performing a particle-hole transformation $c_{i\downarrow}^{+}\to{\left({-1}\right)^{i}}{c_{i\downarrow}}$, it is easy to check that the Monte Carlo weight is always semipositive. In the simulations, we set parameters $\lambda_{x}=\lambda_{y}=\lambda$ as energy unit and $\Delta_{\tau}=0.05/\lambda$. We choose $\Theta=35/\lambda,40/\lambda,45/\lambda,50/\lambda,60/\lambda$ to simulate system with size $L=6,8,10,12,14$, respectively, which guarantees that the projected wave functions converge to ground states.

S2 detailed implementation to calculate $p_{y}^{\nu_{x}}$, $p_{x}^{\nu_{y}}$

Since the model described in Eq.(1) consists of four orbitals with spin degeneracy on square lattice, we take a 4-bands system as an example here to show the detailed implementation that we use to calculate Wannier-sector polarizations in the main text.

A two-dimensional non-interacting system $H=\sum\limits_{k}{c_{k,\alpha}^{+}{{\left[{{h_{k}}}\right]}^{\alpha\beta}}{c_{k,\beta}}}$ can be diagonalized by a linear combination of operator $c_{k,\alpha}$:

[TABLE]

where $n$ is band index. If we set periodic boundary condition with $N_{x}$ sites in $x$-direction, Wilson loop along $k_{x}$ can be defined within the subspace of occupied energy bands as:

[TABLE]

In the 4-band gapped system at half filling, $G_{k}$ is $2\times 2$ matrix constructed by the eigenvectors of occupied bands obtained in Eq.(S5): ${\left[{{G_{\bf{k}}}}\right]^{mn}}=\left\langle{{u_{{\bf{k}}+\Delta{k_{x}}}^{m}}}\mathrel{\left|{\vphantom{{u_{{\bf{k}}+\Delta{k_{x}}}^{m}}{u_{\bf{k}}^{n}}}}\right.\kern-1.2pt}{{u_{\bf{k}}^{n}}}\right\rangle$, where $\left\langle{{u_{{\bf{k}}+\Delta{k_{x}}}^{m}}}\mathrel{\left|{\vphantom{{u_{{\bf{k}}+\Delta{k_{x}}}^{m}}{u_{\bf{k}}^{n}}}}\right.\kern-1.2pt}{{u_{\bf{k}}^{n}}}\right\rangle=\sum\limits_{\alpha}{{{\left[{u_{{\bf{k}}+\Delta{k_{x}}}^{*m}}\right]}^{\alpha}}{{\left[{u_{\bf{k}}^{n}}\right]}^{\alpha}}}$. In the thermodynamic limit, eigenvalues of each Wilson loop are just a phase $E_{x}^{j}(k_{y})=e^{i2\pi\nu_{x}^{j}(k_{y})}$. The phase factor $\nu_{x}^{j}(k_{y})$ depends on the starting point of the Wilson loop is called Wannier center which interpreted as the center of electron distribution in the $x$ direction within a unit cell. After summing over Wannier centor with different occupied bands marked by $j$, we obtain dipole moments in each $k_{y}$: ${\bf{p}}_{x}(k_{y})=\sum_{j}{\nu_{x}^{j}(k_{y})}$. And then, the total dipole moment in the bulk of this 2D system is given by ${\bf{p}}_{x}=1/N_{y}\sum_{k_{y}}{\bf{p}}_{x}(k_{y})$

We should notice that, the definition of the Wilson loop in Eq.(S6) also depends on the starting point of the loop $k_{x}$. Given a definite $k_{y}$, we can define $N_{x}$ Wilson loops $\mathcal{W}_{x}(n\Delta_{k_{x}},k_{y})$ where $n=0,1,2...,N_{x}-1$. Although eigenvalues $\nu_{x}^{j}(k_{y})$ are independent of the starting point $k_{x}$, Wannier functions $\left|{\Psi_{R}^{j}(k_{y})}\right\rangle$ can be constructed by the eigenvectors of those Wilson loops with different starting points:

[TABLE]

where the $k_{x}$ dependent ${\left[{\nu_{x,k_{x}}^{j}}\right]}^{n}$ is the $n$th component of the eigenvector of $2\times 2$ Wilson loop $\mathcal{W}_{x}(k_{x},k_{y})$ with eigenvalue $\nu_{x}^{j}(k_{y})$.

In our 4-band example at half filling, we simply use $\pm$ to mark two values of index $j$. From a series of Wilson loops $\mathcal{W}_{x}(n\Delta_{k_{x}},k_{y})$, we can derive a one-to-one mapping between Wannier functions $\left|{\Psi_{R}^{\pm}(k_{y})}\right\rangle$ and Wannier centers $\nu^{\pm}(k_{y})$ for each $k_{y}$. Then we can make the definition of Wannier bands which describe the dispersion relation $\nu^{\pm}(k_{y})$ versus $k_{y}$. After making the interchange of variables ${k_{x}}\leftrightarrow{k_{y}}$, we can derive another two Wannier bands which stand for $y$-direction polarizations $\nu^{\pm}_{y}(k_{x})$. A vivid sketch of Wannier bands is shown in Ref.Benalcazar et al. (2017).

In the above statement, we first introduce the concept of Wilson loops defined in the subspace of occupied energy bands. Then, electronic polarizations of these bands can be derived from the eigenvalues of them. With the same procedure, we can further define the so-called nested Wilson loops based on each Wannier bands and the corresponding Wannier functions. For example, without providing tedious details, we directly give the expression of nested Wilson loop along $y$-direction within the subspace of Wannier band $\nu_{x}^{+}$ :

[TABLE]

where the vector $|{w_{x,{\bf{k}}}^{+}}\rangle$ is defined as:

[TABLE]

In the thermodynamic limit, The nested Wilson loop we obtained here is just a complex number with module one. We can directly extract observable $p_{y}^{\nu_{x}^{+}}(k_{x})$ from the phase factor which describes the polarization in $y$-direction within the Wannier band $\nu_{x}^{+}$. And the total Wannier-sector polarization $p_{y}^{\nu_{x}^{+}}$ is obtained by summing over all the $p_{y}^{\nu_{x}^{+}}(k_{x})$ with different $k_{x}$ and divided by $N_{x}$. With the same procedure, we can derive all the Wannier-sector poloarizations $p_{x}^{\nu_{y}^{\pm}}$, $p_{y}^{\nu_{x}^{\pm}}$ which are able to distinguish different phases in our model.

S3 quantized topological index due to symmetry

As we have known, in the 1D Su-Schrieffer-Heeger (SSH) model, we can calculate Wannier states and Wannier centers of occupied bands from Wilson loops:

[TABLE]

The bulk charge polarization ${\bf{p}}_{x}$ simply equals to the value of the Wannier center $\nu_{x}$ , which is quantized to [math] or $1/2$ due to the chiral symmetry. Without symmetry breaking, topological phases with ${\bf{p}}_{x}=1/2$ could be distinguished from trivial phases ${\bf{p}}_{x}=0$, which means you can not adiabatically evolve from one phase to the other without energy band gap closing. Quadrupole insulators studied here can be also understood by a similar way.

Back to our model, $M_{x}$, $M_{y}$ impose constraints on allowed values of four Wannier-sector polarizations defined in S2:

[TABLE]

Recall that in 2D systems, inversion $\mathcal{I}$ and $C_{2}$ transform the coordinates in the same way. Another constraint on the Wannier-sector polarizations due to $\mathcal{I}$ (since $C_{2}$ symmetry can be obtained by combining two mirror operators $M_{x}M_{y}$ in 2D systems) is given as follows:

[TABLE]

We can learn from the above realtions that our system possesses zero bulk polarizations in both $x$, $y$ directions, which is the premise of defining quadrupole phases. And we can simply use a vector ${{\bf{P}}^{\nu}}=({p_{x}^{\nu_{y}^{-}},p_{y}^{\nu_{x}^{-}}})$ to indicate different topological phases, since we can obtain $p_{x}^{\nu_{y}^{+}},p_{y}^{\nu_{x}^{+}}$ using the relations shown above.

S4 details to determine the critical point and critical exponents

As we known, it is a tough work to determine the precise value of critical exponents through numerical simulations. In order to derive reliable results, we first make data collapse by assumption $R_{z}(U,L)=f_{R}((U/U_{c}-1)L^{1/\nu})$. Although, the result of regression analysis is insensitive to the variable $\nu$, a precise value of critical point $U_{c}=13.15(2)$ can be convincingly determined. Then, an optimized value of $\nu=0.711$ can be also obtained by minimizing the error of regression analysis by sweeping $\nu$ (Fig.S1(a)). With a definite value $U_{c}=13.15$, we can also implement the same idea to find a proper $\beta$ in $C_{z}(U,L)L^{2\beta/\nu}=f_{C}((U/U_{c}-1)L^{1/\nu})$ with different $\nu$. At this time, a pair of value $(\nu=0.637,\beta=0.389)$ can be determined by minimizing the fitting error. Related results are shown in Fig.S1(b) and an error bar of $\nu$ can be given by the deviation between these two minimum points of $\nu$ (inset of Fig.S1(a) and Fig.S1(b)).

Furthermore, since we can get an optimized $\beta$ by a arbitrary $\nu$ through the regression analysis, we also plot the optimized values of $\beta$ versus $\nu$ in the Fig.S2(a). The relation between $\beta$ and $\nu$ can be described by a linear function $\beta=0.0169+0.572*\nu$. According to our numerical results, we can’t give a certain pair of critical exponents $(\nu,\beta)$ since all the points that fall on the line can give a fair fitting error. However, we can compared the quantum criticality obtained here with these three cases: (a) 3D Heisenberg universality class, which describes transitions between the AF and a featureless gapped state, (b) the disorder-AFM transition in dimerized quantum Heisenberg model and (c) the VBS-AFM deconfined quantum criticality found in $J1$-$J2$ and $J$-$Q$ models Campostrini et al. (2002); Gong et al. (2014); Wenzel et al. (2008); Sandvik (2007). The critical exponents of these three universality classes deviate from the line which stands for the optimized exponents of HOTI-AFMI transitions. The distinct deviation is further shown in Fig.S2(b-d) by making data collapses with the exponents corresponding to each of these universality classes, that strongly suggest the TPT investigated in our work is highly non-trivial.