Measuring the topological phase transition via the single-particle density matrix

TL;DR

This paper presents a method to detect topological phase transitions in interacting fermionic systems using the single-particle density matrix, with a practical tomography scheme for cold atom experiments.

Contribution

It introduces a way to identify topological phases via the single-particle density matrix and proposes an experimental tomography scheme for cold atom systems.

Findings

Berry curvature and Chern number can be extracted from the density matrix.

Topological phase transition points are identifiable through this method.

A feasible tomography scheme for cold atom experiments is designed.

Abstract

We discuss the topological phase transition of the spin-$\frac{1}{2}$ fermionic Haldane model with repulsive on-site interaction. We show that the Berry curvature of the topological Hamiltonian, the first Chern number, and the topological phase transition point can be extracted from the single-particle density matrix for this interacting system. Furthermore, we design a tomography scheme for the single-particle density matrix of interacting fermionic two-band models in experimental realizations with cold atoms in optical lattices.

Figures

Equations

$$H_t^{-1}(\mathbf{k})\simeq \frac{{\rho }_{\mathbf{k}}^T}{{ϵ}_h(\mathbf{k})}+\frac{1-{\rho }_{\mathbf{k}}^T}{{ϵ}_p(\mathbf{k})},$$

$$\mathcal{G}^{\alpha\beta}_{{\bf k},i\omega}=\sum_{\eta}\Big{[}\frac{\langle 0|\hat{c}_{{\bf k}\alpha}|\eta\rangle\langle\eta|\hat{c}^{\dagger}_{{\bf k}\beta}|0\rangle}{i\omega-E_{\eta}}+\frac{\langle\bar{\eta}|\hat{c}_{{\bf k}\alpha}|0\rangle\langle 0|\hat{c}^{\dagger}_{{\bf k}\beta}|\bar{\eta}\rangle}{i\omega+E_{\bar{\eta}}}\Big{]},$$

$${\hat{H}}_0=-t\sum _{⟨ij⟩s}{\hat{c}}_{is}^{†}{\hat{c}}_{js}+\lambda \sum _{⟨⟨ij⟩⟩s}{e}^{i\varphi {\nu }_{ij}}{\hat{c}}_{is}^{†}{\hat{c}}_{js}+m\sum _{is}{\xi }_i{\hat{c}}_{is}^{†}{\hat{c}}_{is},$$

$${\hat{H}}_{\text{int}}=U\sum _i{\hat{n}}_{i\uparrow }{\hat{n}}_{i\downarrow }.$$

$${\gamma }_{\mathbf{k}}$$

$${\rho }_{\mathbf{k}}^T=\frac{1}{2}\sum _{i=0}^3{a}_{\mathbf{k},i}{\sigma }_i,$$

$$N_{\text{TOF}}(\mathbf{k})$$

$$N_{\text{TOF}}^{\text{I}}(\mathbf{k},\tau )\propto 1+\cos (\Omega \tau )a_{\mathbf{k},1}-\sin (\Omega \tau )a_{\mathbf{k},2}.$$

$$a_{\mathbf{k},1}(\tau )$$

$$N_{\text{TOF}}^{\text{II}}(\mathbf{k},\tau )\propto 1+\cos (\Omega \tau )a_{\mathbf{k},1}+\sin (\Omega \tau )a_{\mathbf{k},3}.$$

Full text

Topological insulators are a fascinating new phase without a local order parameter Hassan2010rmp ; Xiao2010rmp . They have been observed in solid-state materials Gehring2013nl , but have also been realized in quantum simulators such as photonic waveguides Rechtsman2013prl and ultracold atoms Aidelsburger2013prl ; Jotzu2014Nat ; Flaschner2016Sci . In two-dimensional systems, topology can be captured by the Chern number (ChN) as the topological index, which is given by the sum of the integral of the Berry curvature in the Brillouin zone over all occupied bands Xiao2010rmp . Topological insulators, which are characterized by a non-zero Chern number, possess robust conducting edge states at their boundaries. The number of edge states is equal to the Chern number (ChN) for noninteracting systems, according to the bulk-edge correspondence Hassan2010rmp . In solid-state systems and photonic systems, the topology is often revealed via the edge states Hassan2010rmp ; Rechtsman2013prl , while in quantum gas experiments, also the Berry curvature can be reconstructed from quench dynamics Duca2015Sci ; Flaschner2016Sci .

Generalized to interacting systems, the ChN is expressed by the Ishikawa-Matsuyama formula in terms of the single-particle Green’s function Ishikawa1986 . It still reflects the number of quasiparticle edge states when the interaction is weak or moderate, even though the bulk-edge correspondence breaks down in some situations with strong interactions Gurarie2011prb ; You2014prb ; He2016prb . On the other hand, it was proven that the ChN can be evaluated by mapping to a noninteracting topological Hamiltonian determined by the zero-frequency Green’s function, $H_{t}\equiv-1/G_{{\bf k},i\omega=0}$ Wang2012prx , or via the quasiparticle Berry curvature Shindou2008 ; Shindou2008prb ; Zheng2018 ; Wong2013 ; Sengupta2015 . Numerical simulations confirm that interaction could induce topologically nontrivial phases for specific systems Abanin2012prl ; Kumar2016prb ; Hofstetter2018 ; Rachel2018 ; Irsigler2018 ; Zheng2018a . However, these conclusions have so far not been confirmed experimentally. The main reason is that it is still unclear which observables correspond to the topological Hamiltonian and the quasiparticle Berry curvature.

In this letter, we consider the half-filled two-band model in a bipartite lattice with repulsive interaction. We illustrate that the Berry curvature of the topological Hamiltonian, the first Chern number, and the phase transition point can be extracted from the single-particle density matrix (SPDM) of the interacting system. The elements of the SPDM are $\rho_{{\bf k},\alpha\beta}=\langle\hat{c}^{\dagger}_{{\bf k}\alpha}\hat{c}_{{\bf k}\beta}\rangle$, where $\hat{c}^{\dagger}_{{\bf k}\alpha}$ and $\hat{c}_{{\bf k}\beta}$ are the fermionic creation and annihilation operators with momentum $\bf k$, and $\alpha,\beta$ represent the pseudospin from A-B sublattice. Furthermore, we develop a scheme of tomography for the SPDM of interacting fermions in two-dimensional optical lattices with a two-sublattice structure. This scheme involves time-of-flight (TOF) imaging of the momentum distribution following different sudden quenches, which can be implemented in cold atom experiments. Our method generalizes the scheme of tomographic measurement of pure or mixed states proposed in Refs. Hauke2014prl ; Ardila2018ar ; Tarnowski2017prl and realized in Ref. Flaschner2016Sci .

The topological Hamiltonian carries the full information on the topology of the interacting system and is theoretically important for understanding the topological phase transition via analogy with the noninteracting system (Wang2012prx, ), yet it is not a physical observable. The following statement builds a link between the topological Hamiltonian and the SPDM for half-filled fermionic two-band systems, which paves the way to probe it experimentally:

If the intrinsic quasiparticle linewidths $\gamma_{p}({\bf k})$ and $\gamma_{h}({\bf k})$ are much smaller than the quasiparticle energy $[\epsilon_{p}({\bf k})>0]$ and the quasihole energy $[\epsilon_{h}({\bf k})<0]$ respectively, i.e., $\gamma_{p}({\bf k})\ll|\epsilon_{p}({\bf k})|$ and $\gamma_{h}({\bf k})\ll|\epsilon_{h}({\bf k})|$, then the inverse of the topological Hamiltonian can be approximated as

[TABLE]

where $\rho^{T}_{\bf k}$ is the transpose of the SPDM, $\bf 1$ is the $2\times 2$ identity matrix.

In order to prove this, we start from the Lehmann representation of the Green’s function at zero temperature

[TABLE]

where $|0\rangle$ is the many-body ground state with zero energy. ${\eta}$ and $\bar{\eta}$ refer to the excitations ($E_{\eta},E_{\bar{\eta}}>0$). For each given momentum, the spectral density is given by the imaginary part of the trace of the retarded Green’s function, \varrho_{\bf k}(\omega)=\sum_{\eta\alpha}\big{[}{\left|\langle 0|\hat{c}_{{\bf k}\alpha}|\eta\rangle\right|^{2}}\delta({\omega-E_{\eta}})+{\left|\langle\bar{\eta}|\hat{c}_{{\bf k}\alpha}|0\rangle\right|^{2}}\delta({\omega+E_{\bar{\eta}}})]. The coefficient $|\langle 0|\hat{c}_{{\bf k}\alpha}|\eta\rangle|^{2}$ or $|\langle\bar{\eta}|\hat{c}_{{\bf k}\alpha}|0\rangle|^{2}$ becomes a nonnegligible contribution only when the energy of the many-body state is near to the quasiparticle energy, i.e., \big{|}E_{\eta}-\epsilon_{p}\big{|}\lesssim\mathcal{O}[\gamma_{p}] or \big{|}E_{\bar{\eta}}+\epsilon_{h}\big{|}\lesssim\mathcal{O}[\gamma_{h}]. When the linewidth is rather small compared to the quasiparticle energy, we have ${1}/{E_{\eta}}\simeq{1}/{\epsilon_{p}}$ and ${1}/{E_{\bar{\eta}}}\simeq-{1}/{\epsilon_{h}}$ for the contribution to $\varrho_{\bf k}$ and $\mathcal{G}^{\alpha\beta}_{{\bf k},i\omega=0}$. By using $\langle 0|\hat{c}^{\dagger}_{{\bf k}\beta}\hat{c}_{{\bf k}\alpha}|0\rangle=\rho_{{\bf k},\beta\alpha}$ and $[\hat{c}_{{\bf k}\alpha},\hat{c}^{\dagger}_{{\bf k}\beta}]=\delta_{\alpha\beta}$, we indeed obtain Eq. (1) from Eq. (2) at zero frequency. The error for this approximation is of order $\gamma_{p(h)}/\epsilon_{p(h)}$.

Eq. (1) shows that $H_{t}({\bf k})$ and $\rho^{T}_{\bf k}$ have exactly the same eigenvectors and the lower band of the former is mapped onto the higher band of the latter. This allows us to obtain the Berry curvature of the topological Hamiltonian through measuring the SPDM. In addition, Eq. (1) still holds when the temperature $T$ is finite but much smaller than the gap. The additional error is suppressed exponentially by $\exp[{-|\epsilon_{p(h)}|/k_{B}T}]$.

Let us consider the spin-$\frac{1}{2}$ Haldane model in a hexagonal optical lattice, which has been realized as a Floquet system in cold atom experiments Jotzu2014Nat ; Flaschner2016Sci . The Hamiltonian reads

[TABLE]

where the first and second terms are the nearest and the next-nearest neighbor hopping terms. $s$ refers to spin $\uparrow$ and $\downarrow$. $\nu_{ij}=\pm 1$, which is related to the hopping path. In the following, we restrict ourselves to the case of $\phi=\pi/2$, which maximally breaks time reversal symmetry. The third term is a staggered potential with $\xi_{i}=1$ for sublattice A and $\xi_{i}=-1$ for sublattice B. The system displays a transition into a normal insulator from the quantum Hall phase when $|m|$ becomes larger than $3\sqrt{3}|\lambda|$. The energy gap of the system is $2|m-3\sqrt{3}\lambda|$ for $m,\lambda>0$. The on-site interaction reads

[TABLE]

The system has SU(2) symmetry in spin space. Note that an interacting Floquet system contains additional subtleties such as micromotion corrections to the interaction Anisimovas2015prb . With our static effective model given by Eqs. (3) and (4), we focus on the high frequency regime, where these corrections are suppressed Eckardt2017rmp . Related interaction effects in static models can be found in Refs. Imriska2016prb ; Wu2016prb ; Vanhala2016prl ; Rubio-Garcia2018njp .

For weak interaction, using the Hartree-Fock (HF) approximation and HF plus the second-order perturbation correction (HF+2nd), respectively, we plot the phase diagram in Fig. 1a for the case $3\sqrt{3}\lambda=0.5t$. The HF approximation yields a renormalized staggered potential, m\rightarrow m+\frac{U}{2}\big{[}\langle\hat{n}_{\text{A}\downarrow}\rangle-\langle\hat{n}_{\text{B}\downarrow}\rangle\big{]}, where $\hat{n}_{\text{A(B)}\downarrow}$ is the number operator of spin down at each site of sublattice A(B). The repulsive interaction effectively smoothens the staggered potential, and induces the topological insulator phase, which is consistent with the result shown in Ref. Vanhala2016prl . For $m=0.6t$, we show the quasiparticle energy at the Dirac point K within HF and HF+2nd approximation in Fig. 1b. The gap of the system is exactly twice this energy due to particle-hole symmetry. The interaction closes the gap and inverts the bands at $U\simeq 0.6t$.

In the following, we confirm that the linewidth $\gamma_{\bf k}$ is rather small for the weak interaction regime. The linewidth of a HF quasiparticle excitation (corresponding to the HF approximation) can be obtained by considering all collision channels with one particle from the lower band (see Fig. 2a). Using Fermi’s Golden Rule, we obtain the linewidth for the quasiparticle state $|{\bf k}+\uparrow\rangle$,

[TABLE]

where $A_{r}$ is the area of the system and $|{\bf k}\pm\uparrow\rangle$ is the eigenstate of the higher (lower) band with spin up within the HF approximation. Each energy level is two-fold spin degenerate. The two outgoing particles occupy states in the higher band, since the lower band is filled. Momentum conservation demands $\bf\tilde{k}=\bf k^{\prime\prime}+\bf k-\bf k^{\prime}$. The $\delta$-function in Eq. (5) stems from energy conservation. The phase space of the final states is constrained by momentum and energy conservation. In particular, for the quasiparticle at the Dirac point K the linewidth vanishes for zero temperature, since all collision channels are forbidden. In comparison, a quasiparticle excitation with a higher energy has a larger linewidth and ratio $\gamma_{\bf k}/\epsilon_{\bf k}$ due to a larger phase space of the final states (see Fig. 2b). The linewidth as a function of interaction strength, Eq. (5), can formally be parameterized as $c_{1}U^{2}(1+c_{2}U+\cdots)$, where the first $U^{2}$ directly arises from $\hat{H}_{\text{int}}$ and the part $c_{2}U$ is due to the interaction dependent HF states. For weak interaction, the linewidth increases quadratically as a function of the interaction strength. In Fig. 2, we show the HF quasiparticle energy and the ratio of the linewidth to the energy for $U=t$. For different interaction strengths, the linewidth has a similar profile in momentum space but with an interaction-dependent rescaling. A large interaction enhances the linewidth, and thus the ratio $\gamma_{\bf k}/\epsilon_{\bf k}$. We find that up to $U=t$, the linewidth is still rather small compared to the energy ($<2.6\%$) for $m\in[0,t]$. A similar conclusion can be drawn for quasihole states. This confirms the validity of the approximation (1).

The ratio $\gamma_{\bf k}/\epsilon_{\bf k}$ also reflects how much the quasiparticle differs from a single-particle pure state. In principle, when the interaction becomes stronger, the deviation increases. On the other hand, also the temperature $T$ can mix states. For $T=0.01t$, we plot the eigenvalues of the SPDM within HF and HF+2nd approximation, respectively, for the K point in Fig. 1c. The position of the gap closing of the SPDM almost coincides with that of the energy in Fig. 1b. This means that the topological phase transition point can be obtained from the gap closing point of the SPDM as expected. The small deviation from the real phase transition point stems from finite $T$ and linewidth at the $\Gamma$ point, respectively.

We have shown that the higher band of $\rho_{\bf k}^{T}$ provides information on topological properties of the lower energy band of the system. In the following, we illustrate how to measure it in cold atom experiments. Including finite temperature and interaction effects, the many-body density matrix of an interacting system is $\mathcal{P}_{M}=\sum_{\eta}p_{\eta}|\eta\rangle\langle\eta|$, where $|\eta\rangle$ is a many-body energy eigenstate and $p_{\eta}$ is the thermal equilibrium probability distribution function with the constraint $\sum_{\eta}p_{\eta}=1$. The SPDM becomes ${\rho}_{{\bf k},\alpha\beta}=\text{Tr}[\hat{c}^{\dagger}_{{\bf k}\alpha}\hat{c}_{{\bf k}\beta}\mathcal{P}_{M}]$. Here and in the following, the spin index is dropped due to SU(2) symmetry. $\alpha$ and $\beta$ are the pseudospin sublattice index (A, B). The transpose of the SPDM can be represented as

[TABLE]

where $\sigma_{1(2,3)}$ is the Pauli matrix and $\sigma_{0}=\bf 1$. The coefficients are $a_{{\bf k},i}={\rm Tr}[\rho^{T}_{\bf{k}}\sigma_{i}]=\sum_{\eta}p_{\eta}\langle\eta|\hat{c}^{\dagger}_{\bf{k}}\sigma_{i}\hat{c}_{\bf{k}}|\eta\rangle$, where $\hat{c}_{\bf{k}}=(\hat{c}_{{\bf k}\text{A}},\hat{c}_{{\bf k}\text{B}})^{T}$. Note that $a_{{\bf k},0}>0$ is the total density $\rho_{{\bf k},\text{AA}}+\rho_{{\bf k},\text{BB}}$, and it equals 1 for the half filling case with particle-hole symmetry.

Quench dynamics can be used to reconstruct the vector ${\bm{a}}_{\bf k}\equiv(a_{{\bf k},1},a_{{\bf k},2},a_{{\bf k},3})$. Let us suppose that the system is suddenly quenched to a new noninteracting Hamiltonian $\hat{\mathcal{H}}=\frac{\Omega}{2}\sum_{\bf k}\hat{c}^{\dagger}_{\bf{k}}{\bm{\sigma}\cdot}{{\bm{d}}_{\bf k}}\hat{c}_{\bf{k}}$ at the time $\tau=0$, where ${{\bm{d}}_{\bf k}}$ is a momentum-dependent unit vector. The coefficients of $\rho^{T}_{\bf k}(\tau)$ become $a_{{\bf{k}},j}(\tau)=\sum_{\eta}p_{\eta}\langle\eta|\hat{c}^{\dagger}_{\bf{k}}V_{\bf k}^{\dagger}(\tau)\sigma_{j}{V_{\bf k}}(\tau)\hat{c}_{\bf{k}}|\eta\rangle$ after evolution to time $\tau>0$, where ${V_{\bf k}}(\tau)=e^{-i\tau(\Omega/2){\bm{\sigma}\cdot}{{\bm{d}}_{\bf k}}}$ is a $2\times 2$ matrix. Since $V_{\bf k}^{\dagger}(\tau)\sigma_{j}{V_{\bf k}}(\tau)$ is a linear combination of $\sigma_{i}$, this formula links the coefficients at time $\tau$ to those at time $\tau=0$. The evolution effectively rotates the vector ${\bm{a}}_{\bf k}$, and $a_{{\bf{k}},0}$ is time-independent after the quench. Thus, the initial SPDM can be deduced from the final coefficients. However, in TOF experiments, not all of the final coefficients can be recorded. The density operator of particles in momentum space observed in TOF experiments is $\hat{N}_{\text{TOF}}({\bf k})\simeq|\tilde{w}({\bf k})|^{2}\sum_{\bf{RR^{\prime}}}e^{-i{\bf{k\cdot(R-R^{\prime})}}}\hat{c}^{\dagger}_{{\bf R^{\prime}}}\hat{c}_{\bf R}$, where $\bf R$ and $\bf R^{\prime}$ denote lattice sites and $\tilde{w}(\bf k)$ is the Fourier transformation of the Wannier function Bloch2008rmp . So the particle density observed is

[TABLE]

where $a_{{\bf k},0}^{F}=a_{{\bf k},0}=1$ and $a^{F}_{{\bf k},1}$ are the components of the final SPDM at the time before the free expansion. Only the component $a_{{\bf k},1}^{F}$ can be detected.

Through rotating the initial vector ${\bm{a}}_{\bf k}$ during the dynamics after the quench, we can reconstruct ${\bm{a}}_{\bf k}$ by detecting its projection onto the first component. In the first protocol, the system is suddenly quenched to $\hat{\mathcal{H}}$ with $\bm{d}_{\bf k}=(0,0,1)$ at the time $\tau=0$, which can be realized by switching off all tunneling and interaction but with the staggered potential ${\Omega}/{2}$ remaining Hauke2014prl . The rotation couples $a_{{\bf k},1}$ and $a_{{\bf k},2}$, and we have $a_{{\bf{k}},1}(\tau)=\cos(\Omega\tau)a_{{\bf{k}},1}-\sin(\Omega\tau)a_{{\bf{k}},2}$. If the atoms are completely released at time $\tau$, then the particle density observed by TOF imaging is

[TABLE]

The protocol is the same as that for a single-particle pure state (SPPS) in noninteracting systems Hauke2014prl . By fitting to the experimental data, both $a_{{\bf{k}},1}$ and $a_{{\bf{k}},2}$ can be obtained from the oscillating behavior of ${N}^{\text{I}}_{\text{TOF}}({\bf k},\tau)$. For the SPPS, $|{\bm{a}}_{{\bf k}}|=1$ so that $|{a}_{{\bf k},3}|$ can be obtained from the known $a_{{\bf{k}},1}$ and $a_{{\bf{k}},2}$. This is not true for a general density matrix where $|{\bm{a}}_{{\bf k}}|\leq 1$. Additional experiments involving components $a_{{\bf{k}},1}$ and $a_{{\bf{k}},3}$ are needed for detecting $a_{{\bf{k}},3}$.

The second protocol uses the quench channel $\hat{\mathcal{H}}$ with ${\bm{d}}({\bf k})=(\cos(k_{y}l+\varphi),\sin(k_{y}l+\varphi),0)$, where $l$ is the lattice constant. This Hamiltonian can be realized by switching on (laser assisted) tunneling ${\Omega e^{-i\varphi}}/{2}$ between A-B sublattices only along $y$-direction, as shown in Fig. 3a. This Hamiltonian induces a similar precession dynamics on the Bloch sphere as the first protocol, but now along a vector, which lies in the xy-plane. The coefficient then becomes

[TABLE]

Using known $a_{{\bf{k}},1}$ and $a_{{\bf{k}},2}$, we can get $a_{{\bf{k}},3}$ by detecting the particle density ${N}^{\text{II}}_{\text{TOF}}({\bf k},\tau)\propto 1+a_{{\bf{k}},1}(\tau)$, except for the points $\bf k$ with $\sin(k_{y}l+\varphi)=0$. At these points, ${\bm{d}}({\bf k})=(1,0,0)$, which cannot generate an effective rotation that couples $a_{{\bf{k}},1}$ and $a_{{\bf{k}},3}$. To obtain $a_{{\bf{k}},3}$ in the whole Brillouin zone, a similar experiment with $\varphi\rightarrow\varphi+\pi/2$ can be implemented. One can choose the two experiments with $\varphi=0$ (normal tunneling) and $\varphi=\pi/2$ (laser assisted tunneling), respectively. Since the second protocol directly accesses $a_{{\bf{k}},3}$, there is no missing information on northern or southern hemisphere as it appears for the quench on flat bands Hauke2014prl .

Note that all the different Hamiltonians discussed above could be realized by starting with a static lattice with large AB-offset and shaking Flaschner2016Sci ; Weitenberg2017ar . Specifically, circular shaking is used for simulating the Haldane model, while asymmetric linear shaking along the y-direction can be used for realizing the situation in Fig. 3a Struck2012prl . For realizing the first protocol, the quench can simply be realized by switching off the shaking, which was used to realize the Haldane model before the quench. The interaction can be switched off by using a Feshbach resonance Chin2010rmp or by tuning the confinement strength along $z$-direction for a transverse confinement optical lattice. The time scale for ramping the interaction to zero should be much smaller than the time scale $1/\Omega$ and $1/U$, so that interaction effects during quench dynamics can be omitted.

The Berry curvature can then be extracted from the known ${\bm{a}}_{{\bf k}}$. Note that in the Fourier transformed basis, $\hat{c}_{{\bf k}\text{A(B)}}\propto\sum_{{\bf R}\in\text{A(B)}}\hat{c}_{\bf R}e^{-i{\bf k}\cdot{\bf R}}$, the Hamiltonian is not periodic but with an additional unitary transformation after translating by a reciprocal lattice vector $\bf Q$. We obtain $h({\bf k+Q})=U^{\dagger}_{\bf\bf Q}h({\bf k})U_{\bf Q}$, where $U_{\bf Q}=\text{diag}(1,e^{-iQ_{y}l})$ and $h({\bf k})$ is the matrix representation of the noninteracting Hamiltonian $\hat{H}_{0}$. Thus we introduce the unitary transformation $(\hat{\tilde{c}}_{{\bf k}\text{A}},\hat{\tilde{c}}_{{\bf k}\text{B}})=(\hat{c}_{{\bf k}\text{A}},\hat{c}_{{\bf k}\text{B}}e^{-ik_{y}l})$ to render the Hamiltonian periodic. The components for the periodic SPDM $\langle\hat{\tilde{c}}^{\dagger}_{{\bf k}\alpha}\hat{\tilde{c}}_{{\bf k}\beta}\rangle$ are ${\tilde{\bm{a}}}_{{\bf k}}=(\cos(k_{y}l)a_{{\bf k},1}+\sin(k_{y}l)a_{{\bf k},2},\cos(k_{y}l)a_{{\bf k},2}-\sin(k_{y}l)a_{{\bf k},1},a_{{\bf k},3})$. We plot the result of ${\bm{\tilde{}}{a}}_{\bf k}$ for different interaction strengths in Fig. 4. The two-dimensional vector $({\tilde{a}}_{{\bf k},1},{\tilde{a}}_{{\bf k},2})$ has an opposite winding behavior circuiting the Dirac points K and K*′*, as in noninteracting systems Hauke2014prl . The third component ${\tilde{a}}_{{\bf k},3}$ for the K point moves from the south pole to the north pole when increasing the interaction strength. It changes its sign when $U\simeq 0.6t$. The vector ${\bm{\tilde{}}{a}}_{\bf k}$ maps the Brillouin zone to a closed curved surface in three-dimensional space. For the noninteracting case, it looks like a deflated ball. Interaction inflates this ball to be round. The condition for a topological phase of the higher (and lower) band of $\rho^{T}_{\bf k}$ is that the origin is enclosed by that surface Niu2010rmp . This coincides with whether $a_{{\bf k},3}$ at the K point is positive. Recall that for a single-particle pure state $|{\bm{a}}_{\bf k}|=1$ and it lives on the surface of the sphere (see Fig. 4a). With interaction and finite temperature, ${\bm{a}}_{\bf k}$ can lie within the sphere and the topological phase transition occurs mildly. The Berry curvature of the higher band of $\rho^{T}_{\bf k}$ can be obtained by using the formula $-\frac{1}{4\pi}(\partial_{k_{x}}{\hat{\tilde{\bm{a}}}}\times\partial_{k_{y}}{\hat{\tilde{\bm{a}}}})\cdot{\hat{\tilde{\bm{a}}}}$, where ${\hat{\tilde{\bm{a}}}}={{\tilde{\bm{a}}}}/|{{\tilde{\bm{a}}}}|$.

To determine the phase transition point, we use the second protocol with $\varphi=\pi/2$. For the K point with momentum ${\bf k}=(4\pi/3\sqrt{3}l,0)$, the particle density observed becomes

[TABLE]

For the K point, $a_{{\bf{k}},1}$ is very small, and thus ${N}^{\text{II}}_{\text{TOF}}$ gets a $\pi$-phase shift when $a_{{\bf{k}},3}$ changes sign. This is shown in Fig. 3b. The point of sign change is exactly the phase transition point.

In conclusion, we have established a link between the SPDM and the topological Hamiltonian, and propose a scheme for detecting the SPDM in experiments. This opens up the possibility to experimentally measure the Berry curvature of the topological Hamiltonian, the first Chern number, and topological phase transitions in the interacting ultracold atom systems. The scheme for measuring the SPDM proposed here can be applied to other A-B sublattice structures. Without particle-hole symmetry, only the rescaled vector $\bm{a}_{\bf k}/a_{{\bf k},0}$ can be obtained by fitting to the experiment. However, this rescaled vector already contains the full topological information of the system. For very strong interaction, where the quasiparticle picture does not hold anymore, the connection between topological Hamiltonian and the SPDM is still an open question. A generalized scheme for systems with more bands (especially if more than one band is occupied) will be the subject of future research.