Many-body multipole indices revealed by the real-space dynamical mean-field theory

The polarization is a fundamental property of insulators, and it has been under intensive study for decades [1, 2, 3, 4, 5, 6, 7]. It is a ground-breaking progress to realize that the change of polarization is closely related to the integral of the Berry curvature in the parameter space [1, 2]. The emergence of topological insulators offers a new platform for the study of polarization, which works as a topological index for the topological phase transition. However, for a long time, a systematic method to reveal the polarization in a correlated insulator remains elusive, even if there are detailed discussion on this formalism [4, 8]. In Resta’s pioneering work [9], a formula for the polarization is proposed, and affords the possibility of exploring the polarization in the interacting case. Furthermore, with the concept of localization length [10], its connection with Kohn’s theory of insulating state [11] has been established. To explore the properties of multipole moments, which offers a distinctive point of view to uncover the nontrivial topological physics, is fundamentally important.

The topics on the polarization and its multipole generalization have become more active with the generalized concept of higher-order topological insulator. The higher-order topological phase, which is a key playground for multipole indices, is characterized by the quantized multipole indices, pioneered by Benalcazar et al. [12]. This novel concept has been identified in different non-crystalline lattices without periodicity [13, 14]. The generalization of Resta’s formula to the multipole case has opened the avenue to study the multipole [15], which can serve as an index to classify the higher-order topological insulators. Even though a couple of issues have been pointed out in this generalization [16], detailed calculations have been done for non-interacting models [17], and its effectiveness has been verified. Its relationship with the Wilson-loop formulation [12, 18] has been elaborated. Furthermore, several ways have been tried out to study multipole indices of correlated insulator. In Ref. [19], the polarization of the Hubbard model is explored with the dynamical mean-field theory (DMFT) [20]. However, it is hard to apply it to the multipole moments. In Ref. [21], correlation effects in quadrupole insulators are studied with the quantum Monte Carlo simulations and topological Hamiltonians. In Ref. [22] a many-body invariant is extracted from Resta’s operator for the polarization with methods of exact diagonalization and infinite density matrix renormalization group. Clearly, these are indirect ways to calculate the multipole indices in correlated systems.

However, a practical procedure for directly calculating the many-body multipole index of correlated insulators is still lacking. Here, we propose a Green’s function formalism in the real space for the multipole indices, which can be perfectly fitted into the scheme of the real-space dynamical mean-field theory (R-DMFT) [23, 24, 25]. We have carried out calculations of multipole indices of the correlated Benalcazar-Bernevig-Hughes (BBH) model with different spatial symmetries. The results are consistent with the symmetry analysis, and shed light on the exploration of multipole indices of the correlated insulator. Furthermore, we would like to point out that our formalism can be adopted to study the higher-order topological phase in crystalline systems with quenched disorder and in non-crystalline systems, and that it also affords a way to compute the localization length of the electron [10], the polarization fluctuation [26], and the quantum geometric tensor [27] in the correlated lattices.

The rest of the manuscript has been organized as follows. In Sec. II, we show our derivation details for the general formula of the multipoles, explain our protocol on how to combine them with the R-DMFT, and point out its further generalization. In Sec. III, we show our demonstrating calculation results for the correlated BBH model with different spatial symmetries. A summary and outlook is presented in Sec. IV.

We propose a Green’s function formalism for the electrical polarization for a crystal, which is ready to be studied in a correlated system,

This is the central result in this work. Here $\tilde{G}_{x}=\alpha_{x}G$. $\alpha_{x}=\mathrm{diag}\{e^{\frac{2\pi ix_{1}}{L_{x}}}-1,\cdots,e^{\frac{2\pi ix_{N}}{L_{x}}}-1\}$, and $G_{x}$ is a matrix of Green’s function with the $i,j$-th element as $\left\langle c_{i}^{\dagger}c_{j}\right\rangle=\frac{1}{\beta}\sum_{n}G_{ij}\left(i\omega_{n}\right)$, where $\omega_{n}=\frac{\left(2n+1\right)\pi}{\beta}$ is the Matsubara frequency, and $c_{i}^{\dagger}$ ($c_{j}$) is the creation (annihilation) operator on the lattice site $i$ ($j$) in a supercell with sizes as $L_{x}$, $L_{x}L_{y}$, and $L_{x}L_{y}L_{z}$ for one, two and three-dimensional cases, respectively. $I$ is the identity matrix. $\mathcal{S}$ is dimension-dependent as $\mathcal{S}=1$, $L_{y}$, $L_{y}L_{z}$ for different dimensions, which is the area perpendicular to the $x$ direction. $P_{bg}=\frac{n_{f}}{\mathcal{S}}\sum_{i}\frac{x_{i}}{L_{x}}$ is the contribution from the background charge, with $n_{f}$ the local density of electrical states. It can be combined with fixed boundary condition, periodic boundary condition and a hybrid boundary condition. The generalization to include the spin flavor is straightforward.

We outline the crucial derivation steps for Eq. (1). Starting from the formula [9, 15], $P_{x}=\frac{1}{2\pi\mathcal{S}}\mathrm{Im}\left[\ln\left\langle\hat{U}_{x}\right\rangle\right]-P_{bg}$, where $\left\langle\hat{U}_{x}\right\rangle=\left\langle GS\right|\hat{U}_{x}\left|GS\right\rangle=\left\langle GS\right|e^{2\pi i\sum_{l}\frac{x_{l}}{L_{x}}\hat{n}_{l}}\left|GS\right\rangle$, we have,

For the electron, we note the nilpotency property of the electron number operator $\hat{n}_{l}^{2}=\hat{n}_{l}$, which is a key observation here, and it follows that $e^{2\pi i\frac{x_{l}}{L_{x}}\hat{n}_{l}}=\sum_{m=0}^{\infty}\frac{1}{m!}\left(2\pi i\frac{x_{l}}{L_{x}}\hat{n}_{l}\right)^{m}=1+\alpha_{l}\hat{n}_{l}$ with $\alpha_{l}=e^{\frac{2\pi ix_{l}}{L_{x}}}-1$. Therefore, we come to

where the last step is obtained by a direct multiplication. We define a $N\times N$ matrix $\tilde{G}_{x}$, with $\left(\tilde{G}_{x}\right)_{ij}=\alpha_{i}\left\langle c_{i}^{\dagger}c_{j}\right\rangle$, ($i,j=1,2,\cdots,N$). A second key observation [28] is that different orders of principal minors for $\det\tilde{G}_{x}$ coincides with different order of terms for $\left\langle\hat{U}_{x}\right\rangle$ in Eq. (4) with Wick’s theorem. As an example, we can see that the second order of principal minors of $\det\tilde{G}_{x}$ are $\alpha_{i}\alpha_{j}\left|\begin{array}[]{cc}\left\langle\hat{n}_{i}\right\rangle&\left\langle c_{i}^{\dagger}c_{j}\right\rangle\\ \left\langle c_{j}^{\dagger}c_{i}\right\rangle&\left\langle\hat{n}_{j}\right\rangle\end{array}\right|$, with $1\leq i<j\leq N$, which coincide with $\sum_{1\leq i<j\leq N}\alpha_{i}\alpha_{j}\left\langle\hat{n}_{i}\hat{n}_{j}\right\rangle$ in Eq. (4) by using Wick’s theorem to the latter form. Furthermore, assuming eigenvalues for $\tilde{G}_{x}$ are $\lambda_{i}$ ($i=1,2,\cdots,N$), we have an identity for $\left\langle\hat{U}_{x}\right\rangle$ with eigenvalues as [29],

Compared with the characteristic polynomials $\det\left(\tilde{G}_{x}-\lambda I\right)=\left(-\lambda\right)^{N}+\left(-\lambda\right)^{N-1}\sum_{i=1}^{N}\lambda_{i}+\left(-\lambda\right)^{N-2}\sum_{1\leq i<j\leq N}\lambda_{i}\lambda_{j}+\cdots$, and setting $\lambda=-1$, we find

Therefore, finally we have the polarization in the formalism of Green’s functions,

With the identity $\mathrm{Tr}\ln A=\ln\det A$, we have obtained Eq. (1). We have a remark here for the numerical implementation. As is known that polarization is an intensive quantity. In practical calculations, one cannot adopt Eq. (7), of which the first term tends to go to zero in the thermodynamic limit in the more than one-dimensional space [30], because numerically the imaginary part of the logarithm function lies in the interval $\left[-\pi,\pi\right)$.

Similarly, we can write down the formula for the quadrupole moment in a crystal,

Here $\tilde{G}_{xy}=\alpha_{xy}G$, and $\alpha_{xy}=\mathrm{diag}\{e^{\frac{2\pi ix_{1}y_{1}}{L_{x}L_{y}}}-1,\cdots,e^{\frac{2\pi ix_{N}y_{N}}{L_{x}L_{y}}}-1\}$. $\mathcal{L}=1$, $L_{z}$ for two- and three-dimensional case, which is the length perpendicular to the $x-y$ plane. $Q_{bg}=\frac{n_{f}}{\mathcal{L}}\sum_{i}\frac{x_{i}y_{i}}{L_{x}L_{y}}$. A generalization to the case of the octupole moment is straightforward.

The key issue in the calculations in Eqs. (1) and (8) is to obtain the matrix of Green’s function $G$, and we are interested in correlated systems described by the fermionic Hubbard model. It is necessary to obtain $G$ of a correlated lattice model for calculating of many-body multipole moments. Real-space dynamical mean-field theory (R-DMFT) [20, 23, 24, 25] is perfectly suitable for this purpose. In this scheme, a strongly correlated system on a lattice is mapped into a group of impurities on lattice sites with the same spatial symmetry, along with hopping between different sites. In this way both the onsite interactions and the inhomogeneous due to hopping are taken into account. One then solves the DMFT problem on all the sites until it is fully converged. Therefore, the interaction effects are taken into account on the DMFT level. Once the R-DMFT loop is converged, we then adopt Eqs. (1) and (8) to compute multipole indices. The flow chart for the whole computation process is shown in Fig. 1, which can be mostly divided into two parts: (1) the R-DMFT part, and (2) the computation of many-body multipole indices. We comment on the efficiency of this algorithm. The most time-consuming part would be the R-DMFT routines, which can be accelerated with parallel implementation and optimized by the spatial symmetry consideration. With convergent Green’s function, the calculations of multipole indices are quite efficient.

To sum up, we have achieved two significant advancements. (1) We have a formula in the formalism of real-space Green’s functions for the multipole indices, which is ready to be combined with the DMFT method, and one can further introduce disorders to lattices to explore the rich physics of disorder effects, quasi-crystal lattice [31], amorphous lattices [32, 33] and so on. (2) With the multipole indices as a bridge, one can certainly study the localization length of the electron [10], the polarization fluctuation [26], and the quantum geometric tensor [27] in the correlated lattices.

In this section, we implement the multipole moment of the correlated BBH model based on Eqs. (1) and (8) with the real-space Green’s function calculated by the R-DMFT, and discuss the correlation effects explicitly.

To sum up, we propose a practical calculation scheme for the multipole moments in correlated materials on the level of R-DMFT. By calculations of BBH with different spatial symmetries due to staggered potentials, and different strengths of interactions we show the effectiveness of our method in exploring the properties of the multipole indices. However, its full potential is not fully unleashed. Our general method is ready to be implemented to explore a disordered correlated lattice. The close relation of polarization with the localization length [10], and the quantum geometric tensor [27] affords us a way to explore these interesting physical quantities. Our demonstrating calculations show that the spatial symmetry is a dominant factor for the quantization of the polarization and quadrupole moments even in the correlating regime, at least in the paramagnetic phase.