Phase transitions in the diagonal ensemble of two-band Chern insulators

The discovery of the quantum Hall effect [1, 2], i.e. a quantized Hall conductance in the ground state which jumps from one plateau to another, inspired the study of characterizing different topological states of matter outside the conventional framework of spontaneous symmetry breaking [3, 4]. Considerable effort has been devoted to understanding topological phases in the ground state. More recently, a lot of attention was devoted to the properties of the topological states driven out of equilibrium, in particular for quantum quenches of the Hamiltonian [5–13].

Consider an isolated system initially in the ground state of a Hamiltonian ${\hat{H}}_{i}$ and suddenly changing the Hamiltonian to ${\hat{H}}_{f}$. The wave function follows a unitary time evolution, while the local observables in the long time limit settle to the prediction of the diagonal ensemble [14], which in some cases can be reduced to a thermal ensemble for nonintegrable systems or a generalized Gibbs ensemble [15, 16] for integrable systems. Topological order or topological invariant cannot be expressed as a local observable. Therefore, its identification in a nonequilibrium state is far from trivial. In the toric code model, the topological entropy which measures the long-range entanglement in the long time limit is found to be the same as its initial value independent from whether the ground states of ${\hat{H}}_{i}$ or ${\hat{H}}_{f}$ are topologically trivial or not [5–7]. This result agrees with a universal argument for gapped spin liquids [17], which states that local unitary transformations do not change the long-range entanglement. Similarly, for the Fermi gas on a honeycomb lattice which essentially simulates the Haldane model, the Chern number is proved to be conserved under unitary evolution [18, 19]. However, in the two-dimensional topological superfluid, the winding number of the retarded Green's function after a quench shows a strong dependence on the post-quench Hamiltonian ${\hat{H}}_{f}$ [8, 9], even if the winding of the Anderson pseudo spin texture is conserved [10]. Here the Chern number and the winding numbers are all topological invariants which equal each other in ground states. Also in the one-dimensional topological superfluid, an analysis of tunneling spectroscopy by coupling the system to an auxiliary thermal bath shows that the topological properties are mostly determined by ${\hat{H}}_{f}$ [11]. But in topological superconductors with proximity-induced superconductivity, the Majorana order parameter [12] or the entanglement spectrum [13] indicate that the quenched state is topologically trivial if ${\hat{H}}_{i}$ and ${\hat{H}}_{f}$ are in different topological phases.

To clarify the issue of topological order far from equilibrium, we appeal to a measurable physical quantity, namely the Hall conductance in Chern insulators. We first study a paradigmatic model, i.e. the Dirac model [20], and then extend our results to a general two-band Chern insulator. We find that the Chern number of the unitarily-evolving wave function is conserved and uniquely determined by ${\hat{H}}_{i}$. However, while the Hall conductance of the quenched state is a continuous function of the energy gap in ${\hat{H}}_{f}$, the derivative of this function displays a logarithmic divergence whenever the Chern number for the ground state of the post-quench Hamiltonian changes. We thus identify a new class of phase transitions with an exotic critical behavior, which is quite different from the orthodox one in which the Hall conductance is discontinuous but its derivative is zero everywhere in the phase diagram. The discrepancy between the Chern number based on unitary time evolution and the Hall conductance is attributed to the fact that the latter must be calculated from the diagonal ensemble, in which the coherence between different eigenstates of ${\hat{H}}_{f}$ in the wave function is lost in the long-time limit.

The contents of the paper are arranged as follows. In section 2, we discuss the Chern number of the unitarily-evolving wave function. We derive the formula for the Hall conductance of the diagonal ensemble in section 3. In section 4, we first discuss the nonanalytic behavior of the Hall conductance in the Dirac model and then extend our results to general two-band Chern insulators. At last, a concluding section summarizes our results.

The Hamiltonian of a two-band Chern insulator in two-dimensions is expressed as

$$\begin{eqnarray}&&\hat{H}=\displaystyle \sum _{\vec{k}}{\hat{c}}_{\vec{k}}^{\dagger }{{ \mathcal H }}_{\vec{k}}{\hat{c}}_{\vec{k}},\end{eqnarray} \tag{ 1 }$$

where ${\hat{c}}_{\vec{k}}={({\hat{c}}_{\vec{k}1},{\hat{c}}_{\vec{k}2})}^{T}$ is the fermionic operator and ${\sum }_{\vec{k}}$ sums over a single Brillouin zone. The single-particle Hamiltonian ${{ \mathcal H }}_{\vec{k}}$ can be decomposed into ${{ \mathcal H }}_{\vec{k}}={\vec{d}}_{\vec{k}}\cdot \vec{\sigma }$, where $\vec{\sigma }$ denotes the Pauli matrices, reading

$$\begin{eqnarray}{\sigma }_{1}=\left(\begin{array}{cc}0 & 1\\ 1 & 0\end{array}\right),\quad {\sigma }_{2}=\left(\begin{array}{cc}0 & -i\\ i & 0\end{array}\right),\quad {\sigma }_{3}=\left(\begin{array}{cc}1 & 0\\ 0 & -1\end{array}\right).\end{eqnarray} \tag{ 2 }$$

The Dirac model is a paradigm for two-band Chern insulators [20]. In the Dirac model, the coefficients of the Pauli matrices are ${\vec{d}}_{\vec{k}}=({k}_{x},{k}_{y},M-{{Bk}}^{2})$ with two parameters M and B, and ${\sum }_{\vec{k}}$ sums over the whole momentum plane. The ground state is well known to be classified by the Chern number $C=\tfrac{1}{2}({\bf{sgn}}(M)+{\bf{sgn}}(B))$, which is quantized and changes only at the phase boundary M = 0 or B = 0. The Hall conductance of the ground state is simply the Chern number in units of ${e}^{2}/h$.

At the time t = 0, we suddenly change the Hamiltonian from ${\hat{H}}_{i}=\hat{H}({M}_{i},{B}_{i})$ to ${\hat{H}}_{f}=\hat{H}({M}_{f},{B}_{f})$. Then the wave function evolves according to $| {\rm{\Psi }}(t)\rangle ={{\rm{e}}}^{-{\rm{i}}{\hat{H}}_{f}t}| {\rm{\Psi }}(0)\rangle ={\prod }_{\vec{k}}\otimes | {u}_{\vec{k}}(t)\rangle$, where $| {u}_{\vec{k}}(t)\rangle$ is the single-particle wave function obeying ${{ \mathcal H }}_{\vec{k}}^{f}| {u}_{\vec{k}}(t)\rangle ={\rm{i}}\tfrac{\partial }{\partial t}| {u}_{\vec{k}}(t)\rangle$. The momentum is a good quantum number both in ${\hat{H}}_{i}$ and ${\hat{H}}_{f}$. Therefore, it is natural to generalize the definition of the Chern number for the time-dependent wave function in the following way:

$$\begin{eqnarray}&&C(t)=\displaystyle \frac{{\rm{i}}}{2\pi }\int {\rm{d}}{\vec{k}}^{2}\left(\left\langle \left.\displaystyle \frac{\partial {u}_{\vec{k}}(t)}{\partial {k}_{x}}\right|\displaystyle \frac{\partial {u}_{\vec{k}}(t)}{\partial {k}_{y}}\right\rangle -{\rm{h.c.}}\right),\end{eqnarray} \tag{ 3 }$$

where the integral is over the whole momentum plane. This real-time Chern number characterizes the topological property of the wave function $| {\rm{\Psi }}(t)\rangle$, and can be re-expressed as $C(t)=\tfrac{{\rm{i}}}{2\pi }\int {\rm{d}}\vec{S}\cdot ({\nabla }_{\vec{k}}\times \vec{A}(t))$, where $\vec{S}$ denotes the k_x-k_y plane oriented in the k_z-direction and $\vec{A}(t)=\langle {u}_{\vec{k}}(t)| {\nabla }_{\vec{k}}| {u}_{\vec{k}}(t)\rangle$ is the Berry connection with its curl in the k_z-direction being $\left(\tfrac{\partial {A}_{y}}{\partial {k}_{x}}-\tfrac{\partial {A}_{x}}{\partial {k}_{y}}\right)$. C(t) is determined by the poles of $\vec{A}(t)$ and must remain quantized at all times since locally deforming $\vec{A}(t)$ cannot change it. In fact, the two poles of $\vec{A}(t)$ at k = 0 and $k=\infty$ have conserved residues under a unitary evolution (see appendix A for the proof), so that for arbitrary ${\hat{H}}_{i}$ and ${\hat{H}}_{f}$ we have $C(t)\equiv C(0)$. The Chern number of the wave function never changes although the system is driven out of equilibrium, which agrees with the previous results in different models [18, 19].

The observation that C(t) is independent of ${\hat{H}}_{f}$ does not imply the absence of nonequilibrium phase transitions because C(t) is not a measurable physical quantity. In this paper, a nonequilibrium phase transition is unambiguously indicated by the nonanalytic behavior of observables as the post-quench Hamiltonian ${\hat{H}}_{f}$ changes. We choose the Hall conductance as the indicator of nonequilibrium phase transitions. Notice that in the ground state the Hall conductance is directly related to the Chern number.

It is well known that the Hall conductance cannot be expressed as the expectation value of a local operator, but must be written as the long-time response to an external electric field in linear response theory. And an important observation is that in the long-time limit the far-from-equilibrium system can be described by the diagonal ensemble. This motivates us to study the Hall conductance in the diagonal ensemble, which is the experimentally relevant setting. In the long-time limit, the off-diagonal terms of the density matrix in the eigenbasis of ${\hat{H}}_{f}$ are averaged out [14]. The time-averaged expectation value of an operator $\hat{O}$ can be expressed as

$$\begin{eqnarray}&&\underset{T\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{T}{\int }_{0}^{T}{\rm{d}}t\langle {\rm{\Psi }}(t)| \hat{O}| {\rm{\Psi }}(t)\rangle =\displaystyle \sum _{E}| \langle E| {\rm{\Psi }}(0)\rangle {| }^{2}\langle E| \hat{O}| E\rangle ={\bf{Tr}}[\hat{O}\hat{\rho }],\end{eqnarray} \tag{ 4 }$$

where $| E\rangle$ are the eigenstates of ${\hat{H}}_{f}$ and $\hat{\rho }$ is diagonal in the basis $| E\rangle$ with the elements $| \langle E| {\rm{\Psi }}(0)\rangle {| }^{2}$. If the long-time limit of $\langle {\rm{\Psi }}(t)| \hat{O}| {\rm{\Psi }}(t)\rangle$ exists, it must be determined by $\hat{\rho }$, the so-called diagonal ensemble [14]. While this argument is based on nondegenerate eigenenergies, the applicability of the diagonal ensemble has also been shown in many integrable quantum many-body models [21, 22].

We build our formalism on the diagonal ensemble with the density matrix written as

$$\begin{eqnarray}&&\hat{\rho }=\displaystyle \prod _{\vec{k}}\otimes \left(\displaystyle \sum _{\alpha =\pm }{n}_{\vec{k}\alpha }| {u}_{\vec{k}\alpha }^{f}\rangle \langle {u}_{\vec{k}\alpha }^{f}| \right),\end{eqnarray} \tag{ 5 }$$

where $| {u}_{\vec{k}\alpha }^{f}\rangle$ is the eigenvector of ${{ \mathcal H }}_{\vec{k}}^{f}$ and $\alpha =\pm$ denotes the upper and lower bands with the positive and negative eigenvalues $\pm | {\vec{d}}_{\vec{k}}^{f}|$, respectively. ${n}_{\vec{k}\alpha }$ is the occupation number of the band α and can be expressed as the overlap ${n}_{\vec{k}\alpha }=| \langle {u}_{\vec{k}\alpha }^{f}| {u}_{\vec{k}-}^{i}\rangle {| }^{2}$, where $| {u}_{\vec{k}-}^{i}\rangle$ is the lower-band eigenvector of the initial Hamiltonian ${{ \mathcal H }}_{\vec{k}}^{i}$, which is in fact the initial wave function. The total occupation at each $\vec{k}$ is conserved to be ${n}_{\vec{k}+}+{n}_{\vec{k}-}\equiv 1$. Equation (5) is obtained by averaging out the off-diagonal elements in $(| {u}_{\vec{k}}(t)\rangle \langle {u}_{\vec{k}}(t)| )$.

Now we calculate the Hall conductance of the diagonal ensemble in linear response theory [23], i.e., we replace the equilibrium density matrix in linear response theory by the diagonal ensemble $\hat{\rho }$. This replacement does not cause any problem in the formalism because $\hat{\rho }$ is time-independent satisfying $[\hat{\rho },{\hat{H}}_{f}]=0$. We can then express the Hall conductance as the current–current correlation in the diagonal ensemble:

$$\begin{eqnarray}&&{\sigma }_{{\rm{H}}}=\underset{\omega \to 0}{\mathrm{lim}}\displaystyle \frac{1}{S\omega }{\int }_{0}^{\infty }{\rm{d}}t{{\rm{e}}}^{{\rm{i}}\omega t}{\bf{Tr}}(\hat{\rho }[{\hat{J}}_{y},{\hat{J}}_{x}(t)]),\end{eqnarray} \tag{ 6 }$$

where S denotes the area of the system and is conventionally set to unity. ${\hat{J}}_{\beta }=e{\sum }_{\vec{k}}{\hat{c}}_{\vec{k}}^{\dagger }\tfrac{\partial {{ \mathcal H }}_{\vec{k}}^{f}}{\partial {k}_{\beta }}{\hat{c}}_{\vec{k}}$ is the current operator along the β-direction with e denoting the charge of the particle. Following the process for obtaining the celebrated TKNN number [3], we re-express the dimensionless Hall conductance ${C}_{\mathrm{neq}}:= {\sigma }_{{\rm{H}}}/({e}^{2}/h)$ as (see appendix B for the derivation)

$$\begin{eqnarray}&&{C}_{\mathrm{neq}}=\displaystyle \frac{{\rm{i}}}{2\pi }\displaystyle \sum _{\alpha }\int {\rm{d}}{\vec{k}}^{2}{n}_{\vec{k}\alpha }\left(\left\langle \left.\displaystyle \frac{\partial {u}_{\vec{k}\alpha }^{f}}{\partial {k}_{x}}\right|\displaystyle \frac{\partial {u}_{\vec{k}\alpha }^{f}}{\partial {k}_{y}}\right\rangle -{\rm{h.c.}}\right),\end{eqnarray} \tag{ 7 }$$

which is the integral of the weighted mixture of Berry curvatures in different bands of the post-quench Hamiltonian. We note that a similar formula has been obtained by studying the unitary evolution of a closed system in Floquet states under an infinitesimal external electric field [24, 25]. In the case of ${\hat{H}}_{i}={\hat{H}}_{f}$ (no quench), the occupation is ${n}_{\vec{k}-}=1$ and ${n}_{\vec{k}+}=0$ everywhere in the Brillouin zone, and C_neq is just the Chern number of the initial state as we expect. But for ${\hat{H}}_{i}\ne {\hat{H}}_{f}$, ${n}_{\vec{k}\alpha }\in [0,1]$ becomes a continuous function of $\vec{k}$ so that C_neq is not quantized any more but can take an arbitrary value.

It is worth comparing the real time Chern number C(t) in equation (3) with the dimensionless Hall conductance C_neq in equation (7). The former reflects the topology of the wave function, being quantized but not measurable, while the latter is a true observable but not quantized. They are both integrals of the Berry curvature, but C(t) is derived from the wave function while C_neq follows from the diagonal ensemble where the coherence is lost.

In the Dirac model, it is straightforward to determine the Hall conductance as (see appendix C for the calculation)

$$\begin{eqnarray}&&{C}_{\mathrm{neq}}={\int }_{0}^{\infty }{\rm{d}}\tilde{k}\displaystyle \frac{(\tilde{k}+({B}_{i}\tilde{k}-{M}_{i})({B}_{f}\tilde{k}-{M}_{f}))({B}_{f}\tilde{k}+{M}_{f})}{4{d}_{\vec{k}}^{i}{({d}_{\vec{k}}^{f})}^{4}}\end{eqnarray} \tag{ 8 }$$

with ${d}_{\vec{k}}^{i/f}=\sqrt{\tilde{k}+{({B}_{i/f}\tilde{k}-{M}_{i/f})}^{2}}$. The Hall conductance C_neq is a function of $({M}_{i},{B}_{i},{M}_{f},{B}_{f})$, i.e., the parameters of ${\hat{H}}_{i}$ and ${\hat{H}}_{f}$. This function satisfies the properties:

$$\begin{eqnarray}&&{C}_{\mathrm{neq}}({M}_{i},{B}_{i},{M}_{f},{B}_{f})={C}_{\mathrm{neq}}({B}_{i},{M}_{i},{B}_{f},{M}_{f})=-{C}_{\mathrm{neq}}(-{M}_{i},-{B}_{i},-{M}_{f},-{B}_{f}).\end{eqnarray} \tag{ 9 }$$

Let us study this function as $({M}_{f},{B}_{f})$ changes, while $({M}_{i},{B}_{i})$ is invariant, i.e., the initial state is fixed. Due to equation (9), we only consider the cases ${M}_{i},{B}_{i}\gt 0$ and ${M}_{i}\gt 0,{B}_{i}\lt 0$. As shown in figure 1, C_neq is a continuous function of $({M}_{f},{B}_{f})$ in the whole parameter space (we strictly prove this in appendix C). Recall that the Hall conductance of the ground state has a jump whenever M or B change sign. By driving the system out of equilibrium, we smoothen the function of the Hall conductance. The shape of ${C}_{\mathrm{neq}}({M}_{f},{B}_{f})$ is insensitive to $({M}_{i},{B}_{i})$. As ${M}_{f},{B}_{f}\gg 0$ (${M}_{f},{B}_{f}\ll 0$), C_neq takes a positive (negative) value, while C_neq is close to zero as M_f and B_f have different signs. Even if the initial state is topologically trivial (see figure 1, the top panel), the Hall conductance is finite as ${\hat{H}}_{f}$ is in the nontrivial regime, but it cannot reach the quantized values $\pm {e}^{2}/h$. When the initial state is nontrivial (see figure 1, the bottom panel), the Hall conductance is suppressed as ${\hat{H}}_{f}$ deviates from ${\hat{H}}_{i}$, and can even change the sign as M_f and B_f both change their signs.

Zoom In Zoom Out Reset image size

While ${C}_{\mathrm{neq}}({M}_{f},{B}_{f})$ is continuous, the key finding is that whenever the post-quench Hamiltonian crosses the boundary at M_f = 0 (B_f = 0), the derivative of the Hall conductance $\tfrac{\partial {C}_{\mathrm{neq}}}{\partial {M}_{f}}$ ($\tfrac{\partial {C}_{\mathrm{neq}}}{\partial {B}_{f}}$) diverges to $+\infty$ in a logarithmic way (see appendix C for the proof):

$$\begin{eqnarray}&&\underset{{M}_{f}\to 0}{\mathrm{lim}}\displaystyle \frac{\partial {C}_{\mathrm{neq}}}{\partial {M}_{f}}\sim \displaystyle \frac{-1}{2| {M}_{i}| }\mathrm{ln}\;| {M}_{f}| ,\\ &&\underset{{B}_{f}\to 0}{\mathrm{lim}}\displaystyle \frac{\partial {C}_{\mathrm{neq}}}{\partial {B}_{f}}\sim \displaystyle \frac{-1}{2| {B}_{i}| }\mathrm{ln}\;| {B}_{f}| .\end{eqnarray} \tag{ 10 }$$

As ${M}_{f}\to 0$, $\tfrac{\partial {C}_{\mathrm{neq}}}{\partial {M}_{f}}$ as a function of $(\mathrm{ln}\;| {M}_{f}| )$ asymptotically approaches a straight line with the slope $-1/(2| {M}_{i}| )$, which is independent of ${B}_{f},{B}_{i}$ and the side from which M_f goes to zero. As ${B}_{f}\to 0$, $\tfrac{\partial {C}_{\mathrm{neq}}}{\partial {B}_{f}}$ has a similar divergence since C_neq is invariant under the exchange of M_f and B_f according to equation (9). We identify a nonequilibrium phase transition when the Chern number for the ground state of ${\hat{H}}_{f}$ changes. The critical behavior of this phase transition is exotic, compared to that of ground-state phase transitions in which the Hall conductance has a zero derivative everywhere but displays a discontinuity at the phase boundary.

This phase transition reveals different nonequilibrium phases which share the common symmetries of the Dirac model. Therefore, the broken symmetry picture does not account for this transition. Interestingly, the topological invariant of the unitarily-evolving wave function C(t) is independent of $({M}_{f},{B}_{f})$, and then fails to characterize different phases in this nonequilibrium phase transition. One can assign the Chern number for the ground-state wave function of ${\hat{H}}_{f}$ to each nonequilibrium phase to distinguish them.

Now let us consider a general two-band Chern insulator in two-dimensions with the Hamiltonian given by equation (1). The coefficient vector ${\vec{d}}_{\vec{k}}=({d}_{1\vec{k}},{d}_{2\vec{k}},{d}_{3\vec{k}})$ is different from model to model. But the nonanalytic behavior of Hall conductance is insensitive to the change of ${\vec{d}}_{\vec{k}}$. Instead, it depends only upon the lowest-order expansion of ${\vec{d}}_{\vec{k}}$ at the momentum $\vec{q}$ where the energy gap closes (${d}_{\vec{q}}=0$) at a phase boundary. In a generic model, two components of ${\vec{d}}_{\vec{q}}$ must be zero. Let us suppose them to be ${d}_{1\vec{q}}$ and ${d}_{2\vec{q}}$ without loss of generality. The energy gap is ${d}_{\vec{q}}=| {d}_{3\vec{q}}|$ when the system is close to the phase boundary. ${d}_{3\vec{q}}$ is a free parameter in the Hamiltonian (the gap parameter), which is denoted by m. Note that m = M in the Dirac model.

Suppose that the system is initially in a ground state with the gap parameter $m={m}_{i}$, before we suddenly change m in the Hamiltonian from m_i to m_f. We measure the Hall conductance in the long time limit. The Hall conductance C_neq is a function of m_f, while we fix m_i to be nonzero. The function ${C}_{\mathrm{neq}}({m}_{f})$ is continuous but nonanalytic at m_f = 0, where the gap of the post-quench Hamiltonian ${\hat{H}}_{f}$ closes. The derivative of ${C}_{\mathrm{neq}}({m}_{f})$ satisfies (see appendix D for the proof)

$$\begin{eqnarray}&&\underset{{m}_{f}\to 0}{\mathrm{lim}}\displaystyle \frac{{\rm{d}}{C}_{\mathrm{neq}}}{{\rm{d}}{m}_{f}}\sim \displaystyle \frac{\underset{{m}_{f}\to {0}^{-}}{\mathrm{lim}}C({m}_{f})-\underset{{m}_{f}\to {0}^{+}}{\mathrm{lim}}C({m}_{f})}{2| {m}_{i}| }\mathrm{ln}\;| {m}_{f}| ,\end{eqnarray} \tag{ 11 }$$

where $C({m}_{f})$ denotes the Chern number in the ground-state wave function of ${\hat{H}}_{f}$. The derivative of the Hall conductance with respect to the energy gap in ${\hat{H}}_{f}$ is logarithmically divergent as the gap closes. And the coefficient of this logarithmic function is the ratio of the change of Chern number for the ground state of ${\hat{H}}_{f}$ to the energy gap in the initial state. Equation (11) relates the nonequilibrium phase transition in quenched states to the topological phase transition in ground states. Equation (10) for the Dirac model is a special case of equation (11) as the change of Chern number is −1.

In summary, we find a new class of phase transitions in quenched states of two-band Chern insulators, which are characterized by the Hall conductance as a continuous function of the energy gap in the post-quench Hamiltonian ${\hat{H}}_{f}$ with a logarithmically divergent derivative. The asymptotic behavior of the Hall conductance is determined by the ratio of the change of the Chern number for the ground state of ${\hat{H}}_{f}$ to the energy gap in the initial state, which is universal in two-band Chern insulators. We obtain the Hall conductance by applying linear response theory to the diagonal ensemble of the system, which is the physically correct description of the long-time limit in a far-from-equilibrium quench setup. The topological invariant of the real-time wave function fails to predict this phase transition, which can only be correctly identified in the diagonal ensemble where decoherence effects are taken into account. Our finding indicates the possibility of exotic topological phase transitions in systems far from equilibrium.

Finally, we discuss the conditions for observing this phase transition in experiments. The nonequilibrium distribution of particles is responsible for the logarithmically divergent derivative of the Hall conductance. Ultracold atomic gases are known to be well isolated from the environment and suitable for studying the quench dynamics of many-body quantum systems [26]. The Haldane model [27] was recently realized with cold atoms in an optical lattice [28]. The Haldane model is a two-band Chern insulator, in which the quenched-state Hall conductance displays the nonanalytic behavior in equation (11). The direct measurement of the Hall conductance in cold atoms is difficult. But the Hall conductance can be obtained from the Faraday rotation angle [29]. Therefore, we expect our prediction to be verified by measuring the Faraday rotation angle of the cold atoms on an optical lattice that simulates the Haldane model.

We thank Professor Q Niu for inspiring discussions. We thank J Oberreuter and M Medvedyeva for their help in preparing the paper. Pei Wang is supported by NSFC under Grant No. 11304280, and by China Scholarship Council. SK was supported through SFB 1073 (project B03) of the Deutsche Forschungsgemeinschaft (DFG).

We express the real-time Chern number as

$$\begin{eqnarray}&&C(t)=\displaystyle \frac{{\rm{i}}}{2\pi }\int {\rm{d}}\vec{S}\cdot ({\nabla }_{\vec{k}}\times \vec{A}(t)),\end{eqnarray} \tag{ A1 }$$

where $\vec{A}(t)=\langle {u}_{\vec{k}}(t)| {\nabla }_{\vec{k}}| {u}_{\vec{k}}(t)\rangle$ is the Berry connection with $| {u}_{\vec{k}}{(t)\rangle =({\phi }_{1},{\phi }_{2})}^{T}$ denoting the single-particle wave function. In the Dirac model, it is straightforward to calculate the wave function and obtain

$$\begin{eqnarray}&&{\phi }_{1}(t)=\displaystyle \frac{1}{\sqrt{2{d}_{\vec{k}}^{i}({d}_{\vec{k}}^{i}-{d}_{3\vec{k}}^{i})}}\left[-\mathrm{cos}({d}_{\vec{k}}^{f}t)({d}_{\vec{k}}^{i}-{d}_{3\vec{k}}^{i})+{\rm{i}}\;\mathrm{sin}({d}_{\vec{k}}^{f}t)\displaystyle \frac{{d}_{3\vec{k}}^{f}({d}_{\vec{k}}^{i}-{d}_{3\vec{k}}^{i})-{k}^{2}}{{d}_{\vec{k}}^{f}}\right]\end{eqnarray} \tag{ A2 }$$

and

$$\begin{eqnarray}&&{\phi }_{2}(t)=\displaystyle \frac{{k}_{+}}{2{d}_{\vec{k}}^{f}\sqrt{2{d}_{\vec{k}}^{i}({d}_{\vec{k}}^{i}-{d}_{3\vec{k}}^{i})}}[2{d}_{\vec{k}}^{f}\mathrm{cos}({d}_{\vec{k}}^{f}t)+2{\rm{i}}\;\mathrm{sin}({d}_{\vec{k}}^{f}t)({d}_{\vec{k}}^{i}-{d}_{3\vec{k}}^{i}+{d}_{3\vec{k}}^{f})],\end{eqnarray} \tag{ A3 }$$

where ${k}_{+}={k}_{x}+{\rm{i}}{k}_{y}$, ${\vec{d}}_{\vec{k}}^{i/f}=({d}_{1\vec{k}}^{i/f},{d}_{2\vec{k}}^{i/f},{d}_{3\vec{k}}^{i/f})$ is the coefficient vector in the initial and post-quench Hamiltonians, respectively, and ${d}_{\vec{k}}^{i/f}$ is the length of ${\vec{d}}_{\vec{k}}^{i/f}$. We divide the Berry connection into $\vec{A}(t)={\vec{A}}_{1}(t)+{\vec{A}}_{2}(t)$ with ${\vec{A}}_{\alpha }={\phi }_{\alpha }^{*}{\nabla }_{\vec{k}}{\phi }_{\alpha }$. Noticing that ${\phi }_{1}$ is a function of $k=\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}$, we immediately know that $({\nabla }_{\vec{k}}\times {\vec{A}}_{1})$ must be zero, so that ${\vec{A}}_{1}$ does not contribute to C(t). We again divide ${\vec{A}}_{2}$ into the irrelevant term ${\vec{A}}_{2a}$ with a zero curl and the relevant term ${\vec{A}}_{2b}$ with its imaginary part written as

$$\begin{eqnarray}&&{\bf{Im}}[{\vec{A}}_{2b}]=\left[{\mathrm{cos}}^{2}({d}_{\vec{k}}^{f}t)+{\mathrm{sin}}^{2}({d}_{\vec{k}}^{f}t)\displaystyle \frac{{({d}_{\vec{k}}^{i}-{d}_{3\vec{k}}^{i}+{d}_{3\vec{k}}^{f})}^{2}}{{({d}_{\vec{k}}^{f})}^{2}}\right]\times \displaystyle \frac{-{k}_{y}\vec{x}+{k}_{x}\vec{y}}{2{d}_{\vec{k}}^{i}({d}_{\vec{k}}^{i}-{d}_{3\vec{k}}^{i})},\end{eqnarray} \tag{ A4 }$$

where $\vec{x}$ and $\vec{y}$ denote the unit vectors in the momentum plane.

Now we re-express the Chern number by the vector field ${\bf{Im}}[{\vec{A}}_{2b}]$ as

$$\begin{eqnarray}&&C(t)=\displaystyle \frac{-1}{2\pi }\int {\rm{d}}\vec{S}\cdot ({\nabla }_{\vec{k}}\times {\bf{Im}}[{\vec{A}}_{2b}(t)]).\end{eqnarray} \tag{ A5 }$$

${\bf{Im}}[{\vec{A}}_{2b}]$ is a vortex field with two poles at zero and infinity, respectively. Applying the Kelvin–Stokes theorem in an annulus with inner radius r and outer radius R, and then taking the limit $r\to 0$ and $R\to \infty$, we obtain

$$\begin{eqnarray}-2\pi C(t) & = & \underset{r\to 0,R\to \infty }{\mathrm{lim}}{\displaystyle \int }_{r\leqslant k\leqslant R}{\rm{d}}\vec{S}\cdot ({\nabla }_{\vec{k}}\times {\bf{Im}}[{\vec{A}}_{2b}])\\ & = & \left(\underset{R\to \infty }{\mathrm{lim}}{\oint }_{k=R}-\underset{r\to 0}{\mathrm{lim}}{\oint }_{k=r}\right){\bf{Im}}[{\vec{A}}_{2b}]\cdot {\rm{d}}\vec{k}\\ & = & \underset{R\to \infty }{\mathrm{lim}}(2\pi R| {\bf{Im}}[{\vec{A}}_{2b}]{| }_{k=R})\\ & & -\underset{r\to 0}{\mathrm{lim}}(2\pi r| {\bf{Im}}[{\vec{A}}_{2b}]{| }_{k=r}),\end{eqnarray} \tag{ A6 }$$

where $| {\bf{Im}}[{\vec{A}}_{2b}]{| }_{k=R}$ denotes the length of the vector ${\bf{Im}}[{\vec{A}}_{2b}]$ at the circle of radius $(k=R)$. The first limit evaluates $(\pi (1-{\bf{sgn}}({B}_{i})))$, while the second limit evaluates $(\pi (1+{\bf{sgn}}({M}_{i})))$, being both time-independent. In other words, the residues of ${\bf{Im}}[{\vec{A}}_{2b}]$ at zero and infinity are both time-invariant, which leads to a conserved Chern number:

$$\begin{eqnarray}&&C(t)\equiv \displaystyle \frac{1}{2}({\bf{sgn}}({M}_{i})+{\bf{sgn}}({B}_{i})).\end{eqnarray} \tag{ A7 }$$

In this section, we first show how to express the Hall conductance of quenched states as the integral of the Berry curvature. Our derivation is a straightforward extension of the work by Thouless et al [3]. We then express the Hall conductance by using the coefficient vectors in two-band Chern insulators.

In linear response theory, the Hall conductance is written as

$$\begin{eqnarray}&&{\sigma }_{{\rm{H}}}=\underset{\omega \to 0}{\mathrm{lim}}\displaystyle \frac{1}{S\omega }{\int }_{0}^{\infty }{\rm{d}}t{{\rm{e}}}^{{\rm{i}}\omega t-\eta | t| }{\bf{Tr}}(\hat{\rho }[{\hat{J}}_{y},{\hat{J}}_{x}(t)]),\end{eqnarray} \tag{ B1 }$$

where η is an infinitesimal number corresponding to the adiabatic switch-on of an external electric field, and ω is the frequency of the electric field with the limit $\omega \to 0$ corresponding to the dc-conductance. The diagonal ensemble is known to be $\hat{\rho }={\prod }_{\vec{k}}\otimes \left({\sum }_{\alpha }{n}_{\vec{k}\alpha }| {u}_{\vec{k}\alpha }^{f}\rangle \langle {u}_{\vec{k}\alpha }^{f}| \right)$, which is a product state. Due to the conversation law ${\sum }_{\alpha }{n}_{\vec{k}\alpha }=1$, the state of the system is limited in a subspace of the Fock space in which the empty or doubly-occupied states at each momentum are excluded. We can then re-express ${\sigma }_{{\rm{H}}}$ in the first-quantization language as

$$\begin{eqnarray}&&{\sigma }_{{\rm{H}}}=\underset{\omega \to 0}{\mathrm{lim}}\displaystyle \frac{1}{S\omega }{\int }_{0}^{\infty }{\rm{d}}t{{\rm{e}}}^{{\rm{i}}\omega t-\eta | t| }\times \displaystyle \sum _{\vec{k},\alpha }{n}_{\vec{k}\alpha }\langle {u}_{\vec{k}\alpha }^{f}| [{\hat{J}}_{\vec{k}}^{y},{\hat{J}}_{\vec{k}}^{x}(t)]| {u}_{\vec{k}\alpha }^{f}\rangle ,\end{eqnarray} \tag{ B2 }$$

where the momentum-resolved current operator is ${\hat{J}}_{\vec{k}}^{\beta }=e\tfrac{\partial {{ \mathcal H }}_{\vec{k}}^{f}}{\partial {k}_{\beta }}$. Since we are interested in the dc Hall conductance which is a real number, we take the real part of ${\sigma }_{{\rm{H}}}$ and obtain

$$\begin{eqnarray}{\bf{Re}}{\sigma }_{{\rm{H}}} & = & \displaystyle \frac{{\sigma }_{{\rm{H}}}+{\sigma }_{{\rm{H}}}^{*}}{2}\\ & = & \displaystyle \frac{-{\rm{i}}{e}^{2}}{S}\displaystyle \sum _{\vec{k},\alpha ,\beta }\displaystyle \frac{{n}_{\vec{k}\alpha }}{{({\epsilon }_{\vec{k}\alpha }-{\epsilon }_{\vec{k}\beta })}^{2}}\\ & & \times \left[\langle {u}_{\vec{k}\alpha }^{f}| \displaystyle \frac{\partial {{ \mathcal H }}_{\vec{k}}^{f}}{\partial {k}_{y}}| {u}_{\vec{k}\beta }^{f}\rangle \langle {u}_{\vec{k}\beta }^{f}| \displaystyle \frac{\partial {{ \mathcal H }}_{\vec{k}}^{f}}{\partial {k}_{x}}| {u}_{\vec{k}\alpha }^{f}\rangle -{\rm{h.c.}}\right],\end{eqnarray} \tag{ B3 }$$

where ${\epsilon }_{\vec{k}\alpha }$ denotes the eigenvalue of ${{ \mathcal H }}_{\vec{k}}^{f}$. We make use of the relation ${{ \mathcal H }}_{\vec{k}}^{f}={\sum }_{\alpha }{\epsilon }_{\vec{k}\alpha }| {u}_{\vec{k}\alpha }^{f}\rangle \langle {u}_{\vec{k}\alpha }^{f}|$ and finally obtain

$$\begin{eqnarray}{C}_{\mathrm{neq}} & = & \displaystyle \frac{{\bf{Re}}{\sigma }_{{\rm{H}}}}{{e}^{2}/h}\\ & = & \displaystyle \frac{{\rm{i}}}{2\pi }\displaystyle \sum _{\alpha }\displaystyle \int {\rm{d}}{\vec{k}}^{2}{n}_{\vec{k}\alpha }\left(\left\langle \left.\displaystyle \frac{\partial {u}_{\vec{k}\alpha }^{f}}{\partial {k}_{x}}\right|\displaystyle \frac{\partial {u}_{\vec{k}\alpha }^{f}}{\partial {k}_{y}}\right\rangle -{\rm{h.c.}}\right).\end{eqnarray} \tag{ B4 }$$

In a two-band Chern insulator with the Hamiltonian ${{ \mathcal H }}_{\vec{k}}={\vec{d}}_{\vec{k}}\cdot \vec{\sigma }$, the Berry curvatures in different bands are opposite to each other. By using this and the conservation law ${n}_{\vec{k}+}+{n}_{\vec{k}-}\equiv 1$, we reexpress equation (B4) as

$$\begin{eqnarray}&&{C}_{\mathrm{neq}}=\int {\rm{d}}{\vec{k}}^{2}\mathrm{cos}\theta \cdot { \mathcal C },\end{eqnarray} \tag{ B5 }$$

where ${ \mathcal C }$ denotes the Berry curvature in the lower-band of the post-quench Hamiltonian ${\hat{H}}_{f}$ and can be expressed as

$$\begin{eqnarray}&&{ \mathcal C }=\displaystyle \frac{\left(\tfrac{\partial {\vec{d}}_{\vec{k}}^{f}}{\partial {k}_{x}}\times \tfrac{\partial {\vec{d}}_{\vec{k}}^{f}}{\partial {k}_{y}}\right)\cdot {\vec{d}}_{\vec{k}}^{f}}{4\pi {({d}_{\vec{k}}^{f})}^{3}},\end{eqnarray} \tag{ B6 }$$

and $\mathrm{cos}\theta$ is the occupation factor defined as

$$\begin{eqnarray}&&\mathrm{cos}\theta := (2{n}_{\vec{k}-}-1)=({\vec{d}}_{\vec{k}}^{f}\cdot {\vec{d}}_{\vec{k}}^{i})/({d}_{\vec{k}}^{f}{d}_{\vec{k}}^{i})\end{eqnarray} \tag{ B7 }$$

with θ denoting the angle between ${\vec{d}}_{\vec{k}}^{i}$ and ${\vec{d}}_{\vec{k}}^{f}$. ${\vec{d}}_{\vec{k}}^{i}$ and ${\vec{d}}_{\vec{k}}^{f}$ are the coefficients of the Pauli matrices in the initial and post-quench Hamiltonians, respectively, and ${d}_{\vec{k}}^{i}$ and ${d}_{\vec{k}}^{f}$ are the length of ${\vec{d}}_{\vec{k}}^{i}$ and ${\vec{d}}_{\vec{k}}^{f}$, respectively.

In this section, we show how to prove the continuity of ${C}_{\mathrm{neq}}({M}_{f},{B}_{f})$ and the logarithmic divergence of its derivative at the phase boundary. We only prove the case at M_f = 0 when B_f is fixed to be nonzero, since ${C}_{\mathrm{neq}}({M}_{f},{B}_{f})$ is invariant under the exchange of M_f and B_f.

In the Dirac model, both ${ \mathcal C }$ and $\mathrm{cos}\theta$ are rotationally invariant in the k_x–k_y plane. We can then carry out the azimuthal integration in equation (B5). By making a substitution $\tilde{k}={k}^{2}$, we express the Hall conductance as

$$\begin{eqnarray}&&{C}_{\mathrm{neq}}={\int }_{0}^{\infty }{\rm{d}}\tilde{k}\mathrm{cos}\theta \cdot { \mathcal C },\end{eqnarray} \tag{ C1 }$$

where the Berry curvature is

$$\begin{eqnarray}&&{ \mathcal C }=\displaystyle \frac{1}{4}\displaystyle \frac{{B}_{f}\tilde{k}+{M}_{f}}{{({d}_{\vec{k}}^{f})}^{3}}\end{eqnarray} \tag{ C2 }$$

and the occupation factor is

$$\begin{eqnarray}&&\mathrm{cos}\theta =\displaystyle \frac{\tilde{k}+({B}_{i}\tilde{k}-{M}_{i})({B}_{f}\tilde{k}-{M}_{f})}{{d}_{\vec{k}}^{i}{d}_{\vec{k}}^{f}}\end{eqnarray} \tag{ C3 }$$

with ${d}_{\vec{k}}^{i/f}=\sqrt{\tilde{k}+{({B}_{i/f}\tilde{k}-{M}_{i/f})}^{2}}$. At ${M}_{f}\ne 0$, we can express the derivative of C_neq as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {C}_{\mathrm{neq}}}{\partial {M}_{f}}={\int }_{0}^{\infty }{\rm{d}}\tilde{k}\displaystyle \frac{\partial (\mathrm{cos}\theta \cdot { \mathcal C })}{\partial {M}_{f}}.\end{eqnarray} \tag{ C4 }$$

A straightforward observation is that both $\mathrm{cos}\theta (\tilde{k})$ and ${ \mathcal C }(\tilde{k})$ are smooth functions for $\tilde{k}\in (0,\infty )$. However, they do not uniformly converge to ${(\mathrm{cos}\theta )}_{{M}_{f}=0}$ or ${{ \mathcal C }}_{{M}_{f}=0}$ as ${M}_{f}\to 0$. The unique singularity is $\tilde{k}=0$, at which we have ${\mathrm{lim}}_{\tilde{k}\to 0}{\mathrm{lim}}_{{M}_{f}\to 0}\mathrm{cos}\theta =0$ but ${\mathrm{lim}}_{{M}_{f}\to 0}{\mathrm{lim}}_{\tilde{k}\to 0}\mathrm{cos}\theta ={\bf{sgn}}({M}_{i}){\bf{sgn}}({M}_{f})$. And ${ \mathcal C }(\tilde{k}=0)={\bf{sgn}}({M}_{f})/(4{M}_{f}^{2})$ is divergent as ${M}_{f}\to 0$. We divide the integral into two parts: ${\int }_{0}^{\infty }{\rm{d}}\tilde{k}={\int }_{0}^{\eta }{\rm{d}}\tilde{k}+{\int }_{\eta }^{\infty }{\rm{d}}\tilde{k}$ with $\eta \gt 0$ a number that can be arbitrarily small. The second integral is a smooth function of M_f, which can be proved by studying the asymptotic behavior of $(\mathrm{cos}\theta \cdot { \mathcal C })$ in the limit $\tilde{k}\to \infty$, or more precisely, by making a substitution $\tilde{k}\to 1/\tilde{k}$ in the integral. In fact, $\tilde{k}=\infty$ is a true singularity at the boundary B_f = 0, where $\tilde{k}=0$ is a regular point, since $\mathrm{cos}\theta$ and ${ \mathcal C }$ are invariant under the substitution $\tilde{k}\leftrightarrow 1/\tilde{k}$ and ${M}_{i/f}\leftrightarrow {B}_{i/f}$. If there is any nonanalytic behavior in the function ${C}_{\mathrm{neq}}({M}_{f})$, it must come from the first integral denoted by ${C}_{\mathrm{neq}}^{\eta }$ next. Interestingly, we can choose an arbitrarily small η so that ${d}_{\vec{k}}^{i}$ in $\mathrm{cos}\theta$ converges to a constant ${d}_{\vec{k}}^{i}=| {M}_{i}|$. We then obtain

$$\begin{eqnarray}&&{C}_{\mathrm{neq}}^{\eta }={\int }_{0}^{\eta }{\rm{d}}\tilde{k}\displaystyle \frac{(\tilde{k}+({B}_{i}\tilde{k}-{M}_{i})({B}_{f}\tilde{k}-{M}_{f}))({B}_{f}\tilde{k}+{M}_{f})}{4| {M}_{i}| {(\tilde{k}+{({B}_{f}\tilde{k}-{M}_{f})}^{2})}^{2}}.\end{eqnarray} \tag{ C5 }$$

The calculation of this integral is straightforward since the integrand is rational.

We express the result as ${C}_{\mathrm{neq}}^{\eta }=F(\eta )-F(0)$ with F denoting the original function. The expression of F is lengthy, but it is an elementary function. $F(\eta )$ is a smooth function of M_f, while F(0) is expressed as

$$\begin{eqnarray}F(0) & = & \displaystyle \frac{1}{8| {M}_{i}| {B}_{f}^{2}}\left[\displaystyle \frac{2{B}_{f}({B}_{i}+{B}_{f}){M}_{f}-{B}_{i}}{\sqrt{1-4{B}_{f}{M}_{f}}}\right.\\ & & \times \mathrm{ln}\displaystyle \frac{(2{B}_{f}^{2}{M}_{f}^{2}-4{B}_{f}{M}_{f}+1)\sqrt{1-4{B}_{f}{M}_{f}}-8{B}_{f}^{2}{M}_{f}^{2}+6{B}_{f}{M}_{f}-1}{{M}_{f}^{2}}\\ & & +{B}_{i}\mathrm{ln}{M}_{f}^{2}+2{B}_{i}-2{B}_{f}].\end{eqnarray} \tag{ C6 }$$

We are interested in F(0) as a function of M_f in the neighborhood of the phase boundary M_f = 0. We notice that $\sqrt{1-4{B}_{f}{M}_{f}}$ can be expanded at M_f = 0 into

$$\begin{eqnarray}&&\sqrt{1-4{B}_{f}{M}_{f}}=1-2{B}_{f}{M}_{f}-2{B}_{f}^{2}{M}_{f}^{2}-4{B}_{f}^{3}{M}_{f}^{3}-10{B}_{f}^{4}{M}_{f}^{4}+{ \mathcal O }({M}_{f}^{5}).\end{eqnarray} \tag{ C7 }$$

We substitute this expression into equation (C6) and obtain

$$\begin{eqnarray}F(0) & = & \displaystyle \frac{2{B}_{i}-2{B}_{f}-{B}_{i}\mathrm{ln}(2{B}_{f}^{4})}{8| {M}_{i}| {B}_{f}^{2}}+\displaystyle \frac{{M}_{f}}{4| {M}_{i}| }\mathrm{ln}{M}_{f}^{2}\\ & & +{ \mathcal O }({M}_{f})-\displaystyle \frac{{B}_{i}}{8| {M}_{i}| {B}_{f}^{2}}\mathrm{ln}(1+{ \mathcal O }({M}_{f}))\\ & & +{ \mathcal O }({M}_{f})\mathrm{ln}(1+{ \mathcal O }({M}_{f}))+{ \mathcal O }({M}_{f}^{2})\mathrm{ln}\;| {M}_{f}| .\end{eqnarray} \tag{ C8 }$$

The first term is independent of M_f. The second term is a continuous function of M_f, but its derivative with respect to M_f is divergent as ${M}_{f}\to 0$. All the other terms are continuous functions of M_f, and their derivatives with respect to M_f are finite at M_f = 0. The asymptotic behavior of $\partial {C}_{\mathrm{neq}}/\partial {M}_{f}$ is uniquely determined by the second term. The function ${C}_{\mathrm{neq}}({M}_{f})$ then asymptotically approaches $(-{M}_{f}\mathrm{ln}\;| {M}_{f}| /(2| {M}_{i}| )+\mathrm{const}.)$ in the limit ${M}_{f}\to 0$. This immediately leads to our results that ${C}_{\mathrm{neq}}({M}_{f})$ is continuous³ and $\partial {C}_{\mathrm{neq}}/\partial {M}_{f}$ is logarithmically divergent as

$$\begin{eqnarray}&&\underset{{M}_{f}\to 0}{\mathrm{lim}}\displaystyle \frac{\partial {C}_{\mathrm{neq}}}{\partial {M}_{f}}\sim \displaystyle \frac{-1}{2| {M}_{i}| }\mathrm{ln}\;| {M}_{f}| .\end{eqnarray} \tag{ C9 }$$

Furthermore, we calculate the Hall conductance by numerically integrate equation (B5). We plot $\tfrac{\partial {C}_{\mathrm{neq}}}{\partial {M}_{f}}$ as a function of $\mathrm{ln}\;| {M}_{f}|$ in figure 2. In the limit ${M}_{f}\to 0$, the curves asymptotically approach straight lines with the slope $-1/2| {M}_{i}|$, which is independent of B_i, B_f and the side from which M_f goes to zero. The numerical result coincides well with our analysis.

Zoom In Zoom Out Reset image size

Let us consider a general two-band Chern insulator in two-dimensions with the Hamiltonian expressed as

$$\begin{eqnarray}&&\hat{H}=\displaystyle \sum _{\vec{k}}{\hat{c}}_{\vec{k}}^{\dagger }{{ \mathcal H }}_{\vec{k}}{\hat{c}}_{\vec{k}},\end{eqnarray} \tag{ D1 }$$

where the single-particle Hamiltonian can be decomposed into ${{ \mathcal H }}_{\vec{k}}={\vec{d}}_{\vec{k}}\cdot \vec{\sigma }$ with $\vec{\sigma }$ denoting the Pauli matrices. Examples include the Dirac model, the Haldane model [27] or the Kitaev honeycomb model in the fermionic basis [30, 31]. The coefficient vector ${\vec{d}}_{\vec{k}}=({d}_{1\vec{k}},{d}_{2\vec{k}},{d}_{3\vec{k}})$ is different from model to model. But the nonanalytic behavior of Hall conductance is insensitive to the change of ${\vec{d}}_{\vec{k}}$, but depends only upon the lowest-order expansion of ${\vec{d}}_{\vec{k}}$ at the singularities of the Berry curvature.

Let us first show how the Chern number of the ground-state wave function is related to the expansion of ${\vec{d}}_{\vec{k}}$. The Chern number is expressed by the Berry connection as

$$\begin{eqnarray}&&C=\displaystyle \frac{-1}{2\pi }\int {\rm{d}}\vec{S}\cdot ({\nabla }_{\vec{k}}\times {\bf{Im}}\vec{A})\end{eqnarray} \tag{ D2 }$$

with

$$\begin{eqnarray}&&{\bf{Im}}\vec{A}=\displaystyle \frac{{d}_{1\vec{k}}{\nabla }_{\vec{k}}{d}_{2\vec{k}}-{d}_{2\vec{k}}{\nabla }_{\vec{k}}{d}_{1\vec{k}}}{2{d}_{\vec{k}}({d}_{\vec{k}}-{d}_{3\vec{k}})}.\end{eqnarray} \tag{ D3 }$$

The Chern number must be zero when $\vec{A}$ has no singularity in the Brillouin zone. According to Kelvin–Stokes theorem, the Chern number can be expressed as the line integral of $\vec{A}$ over the boundaries of the infinitesimal neighborhoods of singularities. Suppose that $\vec{A}$ has a set of singularities ${\vec{q}}_{1},{\vec{q}}_{2},\ldots ,{\vec{q}}_{N}$ in a single Brillouin zone. The Chern number can be expressed as

$$\begin{eqnarray}&&C=\displaystyle \sum _{j=1}^{N}{C}^{({\vec{q}}_{j})}\end{eqnarray} \tag{ D4 }$$

with

$$\begin{eqnarray}&&{C}^{({\vec{q}}_{j})}=\displaystyle \frac{1}{2\pi }\underset{\eta \to 0}{\mathrm{lim}}{\oint }_{\partial {B}_{\eta }({\vec{q}}_{j})}{\bf{Im}}\vec{A}\cdot {\rm{d}}\vec{k},\end{eqnarray} \tag{ D5 }$$

where $\partial {B}_{\eta }({\vec{q}}_{j})$ denotes the boundary of a circle of radius $\sqrt{\eta }$ centered at ${\vec{q}}_{j}$, and the integral is along the anticlockwise direction. Here we do not consider the singularity at infinity, since the Brillouin zone is finite in a generic model.

In general, a singularity of $\vec{A}$ is a point $\vec{q}$ at which $({d}_{\vec{q}}-{d}_{3\vec{q}}=0)$ and then ${d}_{1\vec{q}}={d}_{2\vec{q}}=0$. In a generic model, ${d}_{\vec{q}}=| {d}_{3\vec{q}}|$ is the energy gap when the system is close to the phase boundary. ${d}_{3\vec{q}}$ is a free parameter in the Hamiltonian, which is denoted by m next. Note that m = M in the Dirac model. m is zero if and only if the energy gap closes accompanied by a change of the Chern number.

The Berry connection can be reexpressed as

$$\begin{eqnarray}&&{\bf{Im}}\vec{A}=\left(\displaystyle \frac{{d}_{\vec{k}}+{d}_{3\vec{k}}}{2{d}_{\vec{k}}}\right)\left(\displaystyle \frac{{d}_{1\vec{k}}{\nabla }_{\vec{k}}{d}_{2\vec{k}}-{d}_{2\vec{k}}{\nabla }_{\vec{k}}{d}_{1\vec{k}}}{{({d}_{1\vec{k}})}^{2}+{({d}_{2\vec{k}})}^{2}}\right).\end{eqnarray} \tag{ D6 }$$

Since ${d}_{\vec{q}}$ and ${d}_{3\vec{q}}$ are finite at $m\ne 0$, we can replace $\left(\tfrac{{d}_{\vec{k}}+{d}_{3\vec{k}}}{2{d}_{\vec{k}}}\right)$ by its value at $\vec{k}=\vec{q}$, which is $(1+{\bf{sgn}}(m))/2$ with ${\bf{sgn}}(m)$ denoting the sign of m. This replacement will not change the integral in equation (D5) in the limit $\eta \to 0$. The value of ${d}_{3\vec{k}}$ at $\vec{k}\ne \vec{q}$ has nothing to do with the Chern number.

From equation (D5), we know that the Chern number depends only upon ${\vec{d}}_{\vec{k}}$ around the singularities of $\vec{A}$. We then expand ${d}_{1\vec{k}}$ and ${d}_{2\vec{k}}$ into power series of ${\rm{\Delta }}\vec{k}=\vec{k}-\vec{q}$. Without loss of generality, we have

$$\begin{eqnarray}{d}_{1\vec{k}} & = & {a}_{1x}{\rm{\Delta }}{k}_{x}+{a}_{1y}{\rm{\Delta }}{k}_{y}+{ \mathcal O }({\rm{\Delta }}{k}^{2}),\\ {d}_{2\vec{k}} & = & {a}_{2x}{\rm{\Delta }}{k}_{x}+{a}_{2y}{\rm{\Delta }}{k}_{y}+{ \mathcal O }({\rm{\Delta }}{k}^{2}),\\ {d}_{3\vec{k}} & = & m+{ \mathcal O }({\rm{\Delta }}k).\end{eqnarray} \tag{ D7 }$$

It is straight forward to prove that the higher-order terms in this expansion do not contribute to the integral in equation (D5) in the limit $\eta \to 0$, which evaluates

$$\begin{eqnarray}&&{C}^{(\vec{q})}=\displaystyle \frac{1}{2}(1+{\bf{sgn}}(m)){\bf{sgn}}({a}_{1x}{a}_{2y}-{a}_{2x}{a}_{1y}).\end{eqnarray} \tag{ D8 }$$

It is worth mentioning that the three components of ${\vec{d}}_{\vec{k}}$ are on an equal footing. Depending on the basis that is chosen, the components of ${\vec{d}}_{\vec{k}}$ could be exchanged in some models.

Notice that, in equation (D7), the coefficients a_1x, a_1y, a_2x and a_2y are $\vec{q}$-dependent. While m at different ${\vec{q}}_{j}$ may represent different parameters in the Hamiltonian, i.e. the gap parameters at different phase boundaries. An example is the Haldane model [27]. In a single Brillouin zone, $\vec{A}$ has two singularities. And the energy gap closes at one of them as the system is at some phase boundary, but closes at the other singularity as the system is at the different phase boundary. On the other hand, if the system has some symmetries so that at a specific phase boundary the gap closes simultaneously at several ${\vec{q}}_{j}$, m at these ${\vec{q}}_{j}$ must be the same parameter.

Now let us discuss the Hall conductance of quenched states when the parameters in the initial and post-quench Hamiltonians are both nearby a specific phase boundary where the gap parameter is denoted by m. Suppose that the system is initially in a ground state with the gap parameter $m={m}_{i}$, before we suddenly change m in the Hamiltonian from m_i to m_f. We then measure the Hall conductance in the long time limit. The Hall conductance C_neq is a function of m_f, while we fix m_i to be nonzero.

Noting ${\vec{d}}_{\vec{k}}^{i}={\vec{d}}_{\vec{k}}({m}_{i})$ and ${\vec{d}}_{\vec{k}}^{f}={\vec{d}}_{\vec{k}}({m}_{f})$, we express the Hall conductance as

$$\begin{eqnarray}&&{C}_{\mathrm{neq}}=\displaystyle \frac{1}{4\pi }\int {\rm{d}}{\vec{k}}^{2}\left[\displaystyle \frac{\left(\tfrac{\partial {\vec{d}}_{\vec{k}}({m}_{f})}{\partial {k}_{x}}\times \tfrac{\partial {\vec{d}}_{\vec{k}}({m}_{f})}{\partial {k}_{y}}\right)\cdot {\vec{d}}_{\vec{k}}({m}_{f})}{{({d}_{\vec{k}}({m}_{f}))}^{4}}\times \displaystyle \frac{{\vec{d}}_{\vec{k}}({m}_{i})\cdot {\vec{d}}_{\vec{k}}({m}_{f})}{{d}_{\vec{k}}({m}_{i})}\right],\end{eqnarray} \tag{ D9 }$$

where the integral is over a single Brillouin zone. In a generic model, the components of ${\vec{d}}_{\vec{k}}$ are all analytic functions of $\vec{k}$. According to equation (D9), ${C}_{\mathrm{neq}}({m}_{f})$ is nonanalytic only if ${d}_{\vec{k}}({m}_{f})$ in the denominator of the integral vanishes at some $\vec{k}$, i.e., the singularities of the Berry curvature. This is the case at m_f = 0 when the gap of the post-quench Hamiltonian closes at some singularities of the Berry connection $\vec{A}$. Without loss of generality, we suppose that these singularities are ${\vec{q}}_{1},{\vec{q}}_{2},\ldots ,{\vec{q}}_{{N}^{\prime }}$ with ${N}^{\prime }\leqslant N$. The nonanalyticity of ${C}_{\mathrm{neq}}({m}_{f})$ comes from the integral over the neighborhoods of ${\vec{q}}_{1},{\vec{q}}_{2},\ldots ,{\vec{q}}_{{N}^{\prime }}$. We then divide C_neq into the analytic part and the nonanalytic part as we did in the Dirac model. The latter is written as

$$\begin{eqnarray}&&{C}_{\mathrm{neq}}^{\eta }=\displaystyle \sum _{j=1}^{{N}^{\prime }}{C}_{\mathrm{neq}}^{({\vec{q}}_{j})}\end{eqnarray} \tag{ D10 }$$

with

$$\begin{eqnarray}&&{C}_{\mathrm{neq}}^{({\vec{q}}_{j})}=\displaystyle \frac{1}{4\pi }{\int }_{{B}_{\eta }({\vec{q}}_{j})}{\rm{d}}{\vec{k}}^{2}\left[\displaystyle \frac{\left(\tfrac{\partial {\vec{d}}_{\vec{k}}({m}_{f})}{\partial {k}_{x}}\times \tfrac{\partial {\vec{d}}_{\vec{k}}({m}_{f})}{\partial {k}_{y}}\right)\cdot {\vec{d}}_{\vec{k}}({m}_{f})}{{({d}_{\vec{k}}({m}_{f}))}^{4}}\times \displaystyle \frac{{\vec{d}}_{\vec{k}}({m}_{i})\cdot {\vec{d}}_{\vec{k}}({m}_{f})}{{d}_{\vec{k}}({m}_{i})}\right],\end{eqnarray} \tag{ D11 }$$

where ${B}_{\eta }({\vec{q}}_{j})$ is a circle of radius $\sqrt{\eta }$ centered at ${\vec{q}}_{j}$ with η a positive number that can be arbitrarily small.

In the neighborhood of the singularity $\vec{q}$, we can expand the components of ${\vec{d}}_{\vec{k}}$ into power series. Let us first consider the lowest-order term given by equation (D7). We substitute equation (D7) into (D11). We replace ${d}_{\vec{k}}({m}_{i})$ by its value at $\vec{k}=\vec{q}$, that is ${d}_{\vec{q}}({m}_{i})=| {m}_{i}|$. This replacement will not change the nonanalytic behavior of ${C}_{\mathrm{neq}}^{(\vec{q})}$ since m_i is nonzero. While the denominator of the integrand becomes

$$\begin{eqnarray}&&{({d}_{\vec{k}}({m}_{f}))}^{4}={\left({m}_{f}^{2}+\displaystyle \sum _{j=1}^{2}{({a}_{{jx}}{\rm{\Delta }}{k}_{x}+{a}_{{jy}}{\rm{\Delta }}{k}_{y})}^{2}\right)}^{2}.\end{eqnarray} \tag{ D12 }$$

We change the coordinate system so that the function ${({d}_{\vec{k}}({m}_{f}))}^{4}$ has rotational symmetry around $\vec{q}$. In the new coordinate system we have

$$\begin{eqnarray}&&\displaystyle \sum _{j=1}^{2}{({a}_{{jx}}{\rm{\Delta }}{k}_{x}+{a}_{{jy}}{\rm{\Delta }}{k}_{y})}^{2}={\rm{\Delta }}{{k}^{\prime }}^{2}.\end{eqnarray} \tag{ D13 }$$

This transformation is always possible. Otherwise, the coefficients before ${\rm{\Delta }}{k}_{x}^{2}$ and ${\rm{\Delta }}{k}_{y}^{2}$ have different signs, which contradicts the proposition that $\vec{q}$ is an isolated singularity. In the new coordinate system, we carry out the azimuthal integration and obtain

$$\begin{eqnarray}&&{C}_{\mathrm{neq}}^{(\vec{q})}=\displaystyle \frac{{m}_{f}{\bf{sgn}}({a}_{1x}{a}_{2y}-{a}_{2x}{a}_{1y})}{4| {m}_{i}| }\times {\int }_{0}^{\eta }{\rm{d}}({\rm{\Delta }}{{k}^{\prime }}^{2})\displaystyle \frac{{m}_{i}{m}_{f}+{\rm{\Delta }}{{k}^{\prime }}^{2}}{{({m}_{f}^{2}+{\rm{\Delta }}{{k}^{\prime }}^{2})}^{2}}.\end{eqnarray} \tag{ D14 }$$

In the numerator of the integrand, only the 2nd-order term ${\rm{\Delta }}{{k}^{\prime }}^{2}$ has a contribution to the nonanalyticity of ${C}_{\mathrm{neq}}^{(\vec{q})}({m}_{f})$. It is trivial to find the original function of this integral, whose value is an analytic function of m_f at ${\rm{\Delta }}{{k}^{\prime }}^{2}=\eta$ but a nonanalytic one at ${\rm{\Delta }}{{k}^{\prime }}^{2}=0$. This coincides with our expectation that the nonanalytic behavior of ${C}_{\mathrm{neq}}^{(\vec{q})}({m}_{f})$ should be independent of the choice of η. The nonanalytic part of ${C}_{\mathrm{neq}}^{(\vec{q})}({m}_{f})$ is

$$\begin{eqnarray}&&{C}_{\mathrm{neq}}^{(\vec{q})}\sim \displaystyle \frac{-{\bf{sgn}}({a}_{1x}{a}_{2y}-{a}_{2x}{a}_{1y})}{2| {m}_{i}| }{m}_{f}\mathrm{ln}\;| {m}_{f}| .\end{eqnarray} \tag{ D15 }$$

First, ${C}_{\mathrm{neq}}^{(\vec{q})}$ is a continuous function of m_f, and then the Hall conductance ${C}_{\mathrm{neq}}({m}_{f})$ must be continuous. Second, the derivative of ${C}_{\mathrm{neq}}^{(\vec{q})}$ with respect to m_f is logarithmically divergent in the limit ${m}_{f}\to 0$, i.e.,

$$\begin{eqnarray}&&\underset{{m}_{f}\to 0}{\mathrm{lim}}\displaystyle \frac{{\rm{d}}{C}_{\mathrm{neq}}^{(\vec{q})}}{{\rm{d}}{m}_{f}}\sim \displaystyle \frac{-{\bf{sgn}}({a}_{1x}{a}_{2y}-{a}_{2x}{a}_{1y})}{2| {m}_{i}| }\mathrm{ln}\;| {m}_{f}| .\end{eqnarray} \tag{ D16 }$$

Comparing equation (D8) with (D16), we immediately find that the $\vec{q}$-dependent coefficient in $\tfrac{{\rm{d}}{C}_{\mathrm{neq}}^{(\vec{q})}}{{\rm{d}}{m}_{f}}$ is equal to the change of ${C}^{(\vec{q})}({m}_{f})$ at the phase boundary m_f = 0. ${C}_{\mathrm{neq}}^{\eta }$ is the sum of ${C}_{\mathrm{neq}}^{({\vec{q}}_{j})}$ at the singularities ${\vec{q}}_{1},{\vec{q}}_{2},\ldots ,{\vec{q}}_{{N}^{\prime }}$, while the Chern number C is the sum of ${C}^{({\vec{q}}_{j})}$ at all the singularities of $\vec{A}$. But ${C}^{({\vec{q}}_{j})}$ at $j\gt {N}^{\prime }$ does not change at m_f = 0, since the corresponding gap parameter is different from m_f. We finally obtain

$$\begin{eqnarray}&&\underset{{m}_{f}\to 0}{\mathrm{lim}}\displaystyle \frac{{\rm{d}}{C}_{\mathrm{neq}}}{{\rm{d}}{m}_{f}}\sim \displaystyle \frac{\underset{{m}_{f}\to {0}^{-}}{\mathrm{lim}}C({m}_{f})-\underset{{m}_{f}\to {0}^{+}}{\mathrm{lim}}C({m}_{f})}{2| {m}_{i}| }\mathrm{ln}\;| {m}_{f}| ,\end{eqnarray} \tag{ D17 }$$

which is the central result of this paper.

Equation (D17) is obtained by considering only the lowest-order term in the expansion of ${\vec{d}}_{\vec{k}}$. Next we prove that the higher-order terms do not change the continuity of C_neq or the asymptotic behavior of ${\rm{d}}{C}_{\mathrm{neq}}/{\rm{d}}{m}_{f}$ in the limit ${m}_{f}\to 0$. This is true if the higher-order terms do not change the continuity of ${C}_{\mathrm{neq}}^{(\vec{q})}$ or the asymptotic behavior of ${\rm{d}}{C}_{\mathrm{neq}}^{(\vec{q})}/{\rm{d}}{m}_{f}$ at an arbitrary singularity.

A linear term like $({a}_{3x}{\rm{\Delta }}{k}_{x}+{a}_{3y}{\rm{\Delta }}{k}_{y})$ is not allowed in the expansion of ${d}_{3\vec{k}}$ in equation (D7). Otherwise, ${d}_{\vec{q}}$ is not the energy gap, or the minimum point of ${d}_{\vec{k}}$ is not at $\vec{k}=\vec{q}$, but changes with m, which contradicts our proposition. In a generic model like the Dirac model, the Haldane model or the Kitaev honeycomb model, the minimum point of ${d}_{\vec{k}}$ is determined by the symmetry of the model and then keeps invariant as the system is in the vicinity of the phase boundary.

Let us add the 2nd-order term into ${d}_{3\vec{k}}$, i.e., $({b}_{3x}{\rm{\Delta }}{k}_{x}^{2}+{b}_{3y}{\rm{\Delta }}{k}_{y}^{2}+{b}_{3m}{\rm{\Delta }}{k}_{x}{\rm{\Delta }}{k}_{y})$ without loss of generality. The denominator in the integrand of ${C}_{\mathrm{neq}}^{(\vec{q})}$ becomes

$$\begin{eqnarray}{({d}_{\vec{k}}({m}_{f}))}^{4} & = & \left[{m}_{f}^{2}+\displaystyle \sum _{j=1}^{2}{({a}_{{jx}}{\rm{\Delta }}{k}_{x}+{a}_{{jy}}{\rm{\Delta }}{k}_{y})}^{2}\right.\\ & & +{2{m}_{f}({b}_{3x}{\rm{\Delta }}{k}_{x}^{2}+{b}_{3y}{\rm{\Delta }}{k}_{y}^{2}+{b}_{3m}{\rm{\Delta }}{k}_{x}{\rm{\Delta }}{k}_{y})+{ \mathcal O }({\rm{\Delta }}{k}^{4})]}^{2}.\end{eqnarray} \tag{ D18 }$$

${C}_{\mathrm{neq}}^{(\vec{q})}$ is an integral over the infinitesimal neighborhood of $\vec{q}$, where the 4th-order term ${ \mathcal O }({\rm{\Delta }}{k}^{4})$ is much smaller than the 2nd-order term and can be neglected. At the same time, the additional 2nd-order term that is proportional to m_f has no contribution to the asymptotic behavior of ${C}_{\mathrm{neq}}^{(\vec{q})}$ and ${\rm{d}}{C}_{\mathrm{neq}}^{(\vec{q})}/{\rm{d}}{m}_{f}$ in the limit ${m}_{f}\to 0$. Therefore, the effective denominator is the same as equation (D12). The numerator of the integrand becomes

$$\begin{eqnarray}&&\left[\left(\displaystyle \frac{\partial {\vec{d}}_{\vec{k}}({m}_{f})}{\partial {k}_{x}}\times \displaystyle \frac{\partial {\vec{d}}_{\vec{k}}({m}_{f})}{\partial {k}_{y}}\right)\cdot {\vec{d}}_{\vec{k}}({m}_{f})\right]({\vec{d}}_{\vec{k}}({m}_{i})\cdot {\vec{d}}_{\vec{k}}({m}_{f}))\\ &&\quad =({a}_{1x}{a}_{2y}-{a}_{2x}{a}_{1y})\left[{m}_{i}{m}_{f}^{2}+{m}_{f}\displaystyle \sum _{j=1}^{2}{({a}_{{jx}}{\rm{\Delta }}{k}_{x}+{a}_{{jy}}{\rm{\Delta }}{k}_{y})}^{2}\right.\\ &&\quad +{m}_{f}^{2}({b}_{3x}{\rm{\Delta }}{k}_{x}^{2}+{b}_{3y}{\rm{\Delta }}{k}_{y}^{2}+{b}_{3m}{\rm{\Delta }}{k}_{x}{\rm{\Delta }}{k}_{y})+{ \mathcal O }({\rm{\Delta }}{k}^{4})\phantom{\rule{0ex}{4ex}}].\end{eqnarray} \tag{ D19 }$$

The 4th-order term ${ \mathcal O }({\rm{\Delta }}{k}^{4})$ can be neglected in the limit $\eta \to 0$. This can be easily verified by adding ${\rm{\Delta }}{{k}^{\prime }}^{4}$ in the numerator of the integrand in equation (D14) and checking the output. The additional 2nd-order term that is proportional to m_f² leads to a correction of ${C}_{\mathrm{neq}}^{(\vec{q})}\sim {m}_{f}^{2}\mathrm{ln}| {m}_{f}|$, which does not change the asymptotic behavior of ${C}_{\mathrm{neq}}^{(\vec{q})}$ and ${\rm{d}}{C}_{\mathrm{neq}}^{(\vec{q})}/{\rm{d}}{m}_{f}$ in the limit ${m}_{f}\to 0$. In the power series of ${d}_{3\vec{k}}$, any term in order higher than two leads to a correction to numerator or denominator of the integrand which is at least in the 3rd order of ${\rm{\Delta }}k$ and can then be neglected in the limit $\eta \to 0$. Therefore, the higher-order terms in ${d}_{3\vec{k}}$ do not affect the asymptotic behavior of ${\rm{d}}{C}_{\mathrm{neq}}^{(\vec{q})}/{\rm{d}}{m}_{f}$ or the continuity of ${C}_{\mathrm{neq}}^{(\vec{q})}$.

Similarly, we can prove that the higher-order terms in ${d}_{1\vec{k}}$ or ${d}_{2\vec{k}}$ have no contribution. In fact, the terms in order higher than one lead to a correction of ${ \mathcal O }({\rm{\Delta }}{k}^{3})$ in the denominator. The terms in order higher than three also lead to a correction of ${ \mathcal O }({\rm{\Delta }}{k}^{3})$ in the numerator, which can be neglected. The 2nd- and 3rd-order terms in ${d}_{1\vec{k}}$ or ${d}_{2\vec{k}}$ generate a linear term and a 2nd-order term that is proportional to m_f² in the numerator. The latter does not contribute to the asymptotic behavior of ${\rm{d}}{C}_{\mathrm{neq}}^{(\vec{q})}/{\rm{d}}{m}_{f}$ due to the similar reason mentioned above. While the linear term in the numerator is an odd function of ${\rm{\Delta }}{k}_{x}$ or ${\rm{\Delta }}{k}_{y}$, and then has no contribution to the integral since both the denominator and the integration boundary have rotational symmetry with respect to the singularity.