Emergent topology and symmetry-breaking order in correlated quench dynamics

TL;DR

This paper develops a dynamical theory to characterize topology and symmetry-breaking order in correlated quantum systems through quenching the Haldane-Hubbard model, revealing universal behaviors on band-inversion surfaces that encode these properties.

Contribution

It introduces a novel approach using quench dynamics and flow equation methods to extract topological and magnetic order information in correlated systems beyond mean-field approximations.

Findings

Universal behaviors on band-inversion surfaces (BISs) reveal topology and magnetic order.

Dynamical topological patterns are robust against dephasing and heating.

Universal scaling behavior of quench dynamics indicates symmetry-breaking orders.

Abstract

Quenching a quantum system involves three basic ingredients: the initial phase, the post-quench target phase, and the non-equilibrium dynamics which carries the information of the former two. Here we propose a dynamical theory to characterize both the topology and symmetry-breaking order in correlated quantum system, through quenching the Haldane-Hubbard model from an initial magnetic phase to topologically nontrivial regime. The equation of motion for the complex pseudospin dynamics is obtained with the flow equation method, with the pseudospin evolution shown to obey a microscopic Landau-Lifshitz-Gilbert-Bloch equation. We find that the correlated quench dynamics exhibit robust universal behaviors on the so-called band-inversion surfaces (BISs), from which the nontrivial topology and magnetic order can be extracted. In particular, the topology of the post-quench regime can be characterized by an emergent dynamical topological pattern of quench dynamics on BISs, which is robust against dephasing and heating induced by interactions; the pre-quench symmetry-breaking orders is read out from a universal scaling behavior of the quench dynamics emerging on the BIS, which is valid beyond the mean-field regime. This work opens a way to characterize both the topology and symmetry-breaking orders by correlated quench dynamics.

Figures

Equations

$$H$$

$$H_0$$

$$+e^{-\text{mathbbm}i\varphi }{b}_{\vec{i}\sigma }^{†}{b}_{\vec{j}\sigma }+h.c.)+M\sum _{\vec{i},\sigma }(a_{\vec{i}\sigma }^{†}{a}_{\vec{i}\sigma }-b_{\vec{i}\sigma }^{†}{b}_{\vec{i}\sigma }).$$

$$H(l)$$

$$\sum _{\begin{array}{c}{\mathbf{p}}^{\prime }\mathbf{p}{\mathbf{q}}^{\prime }\mathbf{q}\\ s_1{s}_2{s}_3{s}_4\end{array}}{U}_{{\mathbf{p}}^{\prime }\mathbf{p}{\mathbf{q}}^{\prime }\mathbf{q}}^{s_1{s}_2{s}_3{s}_4}(l):c_{{\mathbf{p}}^{\prime },s_1\uparrow }^{†}{c}_{\mathbf{p},s_2\uparrow }{c}_{{\mathbf{q}}^{\prime },s_3\downarrow }^{†}{c}_{\mathbf{q},s_4\downarrow }:,$$

$$\frac{d{\mathbf{S}}_{\sigma }(t)}{dt}={\mathbf{S}}_{\sigma }(t)\times 2\mathbf{h}-{\eta }_1^{\sigma }{\mathbf{S}}_{\sigma }(t)\times \frac{d{\mathbf{S}}_{\sigma }(t)}{dt}-{\eta }_2^{\sigma }\frac{{\mathbf{S}}_{\sigma }(t)}{T_{\mathrm{g}}},$$

$${\mathbf{S}}_{\sigma }(t)={\mathbf{S}}_{\sigma }^{(0)}+{\mathbf{S}}_{\sigma }^{(c)}(t)+{\mathbf{S}}_{\sigma }^{(h)}(t)+{\mathbf{S}}_{\sigma }^{(l)}(t),$$

$$f(m_{\sigma })=-\frac{\mathrm{sgn}(Z_0^{\sigma })}{g(Z_0^{\sigma })}+O(U^4),$$

$$\mathcal{H}(\mathbf{k})=\mathbf{h}(\mathbf{k})\cdot \mathbf{\tau }=h_x(\mathbf{k}){\tau }_x+h_y(\mathbf{k}){\tau }_y+h_z(\mathbf{k}){\tau }_z,$$

$${\mathcal{E}}_{\pm }(\mathbf{k})=\pm \sqrt{h_x^2(\mathbf{k})+h_y^2(\mathbf{k})+h_z^2(\mathbf{k})}\equiv \pm E_0(\mathbf{k}).$$

$$a_{\mathbf{k}\sigma }={\chi }_{+}(\mathbf{k})c_{\mathbf{k},+\sigma }+{\chi }_{-}(\mathbf{k})c_{\mathbf{k},-\sigma },b_{\mathbf{k}\sigma }={\xi }_{+}(\mathbf{k})c_{\mathbf{k},+\sigma }+{\xi }_{-}(\mathbf{k})c_{\mathbf{k},-\sigma }.$$

$$a_{{\mathbf{p}}^{\prime }\uparrow }^{†}{a}_{\mathbf{p}\uparrow }{a}_{{\mathbf{q}}^{\prime }\downarrow }^{†}{a}_{\mathbf{q}\downarrow }=$$

$$\left[{\chi }_{+}^{\ast }({\mathbf{q}}^{\prime })c_{{\mathbf{q}}^{\prime },+\downarrow }^{†}+{\chi }_{-}^{\ast }({\mathbf{q}}^{\prime })c_{{\mathbf{q}}^{\prime },-\downarrow }^{†}\right]\left[{\chi }_{+}(\mathbf{q})c_{\mathbf{q},+\downarrow }+{\chi }_{-}(\mathbf{q})c_{\mathbf{q},-\downarrow }\right]$$

$$=$$

$$b_{{\mathbf{p}}^{\prime }\uparrow }^{†}{b}_{\mathbf{p}\uparrow }{b}_{{\mathbf{q}}^{\prime }\downarrow }^{†}{b}_{\mathbf{q}\downarrow }=\sum _{s_1{s}_2{s}_3{s}_4}{\xi }_{s_1}^{\ast }({\mathbf{p}}^{\prime }){\xi }_{s_2}(\mathbf{p}){\xi }_{s_3}^{\ast }({\mathbf{q}}^{\prime }){\xi }_{s_4}(\mathbf{q})c_{{\mathbf{p}}^{\prime },s_1\uparrow }^{†}{c}_{\mathbf{p},s_2\uparrow }{c}_{{\mathbf{q}}^{\prime },s_3\downarrow }^{†}{c}_{\mathbf{q},s_4\downarrow }$$

$$H_I=U\sum _{{\mathbf{p}}^{\prime }\mathbf{p}{\mathbf{q}}^{\prime }\mathbf{q}}{\delta }_{{\mathbf{p}}^{\prime }+{\mathbf{q}}^{\prime }}^{\mathbf{p}+\mathbf{q}}(a_{{\mathbf{p}}^{\prime }\uparrow }^{†}{a}_{\mathbf{p}\uparrow }{a}_{{\mathbf{q}}^{\prime }\downarrow }^{†}{a}_{\mathbf{q}\downarrow }+b_{{\mathbf{p}}^{\prime }\uparrow }^{†}{b}_{\mathbf{p}\uparrow }{b}_{{\mathbf{q}}^{\prime }\downarrow }^{†}{b}_{\mathbf{q}\downarrow })=\sum _{\begin{array}{c}{\mathbf{p}}^{\prime }\mathbf{p}{\mathbf{q}}^{\prime }\mathbf{q}\\ s_1{s}_2{s}_3{s}_4\end{array}}{U}_{{\mathbf{p}}^{\prime }\mathbf{p}{\mathbf{q}}^{\prime }\mathbf{q}}^{s_1{s}_2{s}_3{s}_4}{c}_{{\mathbf{p}}^{\prime },s_1\uparrow }^{†}{c}_{\mathbf{p},s_2\uparrow }{c}_{{\mathbf{q}}^{\prime },s_3\downarrow }^{†}{c}_{\mathbf{q},s_4\downarrow },$$

$${\Lambda }_{{\mathbf{p}}^{\prime }\mathbf{p}{\mathbf{q}}^{\prime }\mathbf{q}}^{s_1{s}_2{s}_3{s}_4}\stackrel{\mathrm{def}}{=}{\chi }_{s_1}^{\ast }({\mathbf{p}}^{\prime }){\chi }_{s_2}(\mathbf{p}){\chi }_{s_3}^{\ast }({\mathbf{q}}^{\prime }){\chi }_{s_4}(\mathbf{q})+{\xi }_{s_1}^{\ast }({\mathbf{p}}^{\prime }){\xi }_{s_2}(\mathbf{p}){\xi }_{s_3}^{\ast }({\mathbf{q}}^{\prime }){\xi }_{s_4}(\mathbf{q}).$$

$$H_I=U\sum _i\left(⟨a_{i\uparrow }^{†}{a}_{i\uparrow }⟩a_{i\downarrow }^{†}{a}_{i\downarrow }+a_{i\uparrow }^{†}{a}_{i\uparrow }⟨a_{i\downarrow }^{†}{a}_{i\downarrow }⟩+⟨b_{i\uparrow }^{†}{b}_{i\uparrow }⟩b_{i\downarrow }^{†}{b}_{i\downarrow }+b_{i\uparrow }^{†}{b}_{i\uparrow }⟨b_{i\downarrow }^{†}{b}_{i\downarrow }⟩\right)$$

$$\displaystyle{\cal H}_{\rm MF}=\left(\begin{array}[]{cc}{\cal H}({\bf k})+m_{\uparrow}\tau_{z}&\\ &{\cal H}({\bf k})+m_{\downarrow}\tau_{z}\end{array}\right),$$

$$c_{\mathbf{k},+\sigma }=f_{++}^{\sigma }(\mathbf{k}){\bar{c}}_{\mathbf{k},+\sigma }+f_{+-}^{\sigma }(\mathbf{k}){\bar{c}}_{\mathbf{k},-\sigma },c_{\mathbf{k},-\sigma }=f_{-+}^{\sigma }(\mathbf{k}){\bar{c}}_{\mathbf{k},+\sigma }+f_{--}^{\sigma }(\mathbf{k}){\bar{c}}_{\mathbf{k},-\sigma },$$

$$|f_{++}^{\sigma }{|}^2+|f_{+-}^{\sigma }{|}^2=|f_{-+}^{\sigma }{|}^2+|f_{--}^{\sigma }{|}^2=1,$$

$$|f_{++}^{\sigma }{|}^2=|f_{--}^{\sigma }{|}^2,|f_{+-}^{\sigma }{|}^2=|f_{-+}^{\sigma }{|}^2,$$

$$f_{++}^{\sigma \ast }{f}_{-+}^{\sigma }+f_{+-}^{\sigma \ast }{f}_{--}^{\sigma }=0.$$

$$⟨c_{\mathbf{k},s_1\sigma }^{†}{c}_{\mathbf{k},s_2\sigma }⟩=f_{s_1-}^{\sigma \ast }(\mathbf{k})f_{s_2-}^{\sigma }(\mathbf{k}),⟨c_{\mathbf{k},s_1\sigma }{c}_{\mathbf{k},s_2\sigma }^{†}⟩=f_{s_1+}^{\sigma }(\mathbf{k})f_{s_2+}^{\sigma \ast }(\mathbf{k}).$$

$$H(l)$$

$$[:c_{\mathbf{k},s\sigma }^{†}{c}_{\mathbf{k},s\sigma }:,:c_{{\mathbf{p}}^{\prime },s_1\uparrow }^{†}{c}_{\mathbf{p},s_2\uparrow }{c}_{{\mathbf{q}}^{\prime },s_3\downarrow }^{†}{c}_{\mathbf{q},s_4\downarrow }:]=\left(-{\delta }_s^{s_2}{\delta }_{\sigma }^{\uparrow }{\delta }_{\mathbf{k}}^{\mathbf{p}}-{\delta }_s^{s_4}{\delta }_{\sigma }^{\downarrow }{\delta }_{\mathbf{k}}^{\mathbf{q}}+{\delta }_s^{s_1}{\delta }_{\sigma }^{\uparrow }{\delta }_{\mathbf{k}}^{{\mathbf{p}}^{\prime }}+{\delta }_s^{s_3}{\delta }_{\sigma }^{\downarrow }{\delta }_{\mathbf{k}}^{{\mathbf{q}}^{\prime }}\right):c_{{\mathbf{p}}^{\prime },s_1\uparrow }^{†}{c}_{\mathbf{p},s_2\uparrow }{c}_{{\mathbf{q}}^{\prime },s_3\downarrow }^{†}{c}_{\mathbf{q},s_4\downarrow }:$$

$$\eta (l)$$

$$[\eta (l),H_0(l)]=-\sum _{\begin{array}{c}{\mathbf{p}}^{\prime }\mathbf{p}{\mathbf{q}}^{\prime }\mathbf{q}\\ s_1{s}_2{s}_3{s}_4\end{array}}{U}_{{\mathbf{p}}^{\prime }\mathbf{p}{\mathbf{q}}^{\prime }\mathbf{q}}^{s_1{s}_2{s}_3{s}_4}(l)({\Delta }_{{\mathbf{p}}^{\prime }\mathbf{p}{\mathbf{q}}^{\prime }\mathbf{q}}^{s_1{s}_2{s}_3{s}_4}{)}^2:c_{{\mathbf{p}}^{\prime },s_1\uparrow }^{†}{c}_{\mathbf{p},s_2\uparrow }{c}_{{\mathbf{q}}^{\prime },s_3\downarrow }^{†}{c}_{\mathbf{q},s_4\downarrow }:,$$

$$U_{{\mathbf{p}}^{\prime }\mathbf{p}{\mathbf{q}}^{\prime }\mathbf{q}}^{s_1{s}_2{s}_3{s}_4}(l)=U{\delta }_{{\mathbf{p}}^{\prime }+{\mathbf{q}}^{\prime }}^{\mathbf{p}+\mathbf{q}}{\Lambda }_{{\mathbf{p}}^{\prime }\mathbf{p}{\mathbf{q}}^{\prime }\mathbf{q}}^{s_1{s}_2{s}_3{s}_4}\exp [-l({\Delta }_{{\mathbf{p}}^{\prime }\mathbf{p}{\mathbf{q}}^{\prime }\mathbf{q}}^{s_1{s}_2{s}_3{s}_4}{)}^2],$$

$$[:c_{{\mathbf{p}}^{\prime },s_1\uparrow }^{†}{c}_{\mathbf{p},s_2\uparrow }{c}_{{\mathbf{q}}^{\prime },s_3\downarrow }^{†}{c}_{\mathbf{q},s_4\downarrow }:,c_{\mathbf{k},s\sigma }^{†}]={\delta }_s^{s_2}{\delta }_{\sigma }^{\uparrow }{\delta }_{\mathbf{k}}^{\mathbf{p}}:c_{{\mathbf{p}}^{\prime },s_1\uparrow }^{†}{c}_{{\mathbf{q}}^{\prime },s_3\downarrow }^{†}{c}_{\mathbf{q},s_4\downarrow }:+{\delta }_s^{s_4}{\delta }_{\sigma }^{\downarrow }{\delta }_{\mathbf{k}}^{\mathbf{q}}:c_{{\mathbf{p}}^{\prime },s_1\uparrow }^{†}{c}_{\mathbf{p},s_2\uparrow }{c}_{{\mathbf{q}}^{\prime },s_3\downarrow }^{†}:,$$

$${\mathcal{A}}_{\mathbf{k}\uparrow }^{†}(l)$$

$${\mathcal{A}}_{\mathbf{k}\downarrow }^{†}(l)$$

$$[:c_{{\mathbf{p}}^{\prime },s_1\uparrow }^{†}{c}_{\mathbf{p},s_2\uparrow }{c}_{{\mathbf{q}}^{\prime },s_3\downarrow }^{†}{c}_{\mathbf{q},s_4\downarrow }:,:c_{1^{\prime },\mu \uparrow }^{†}{c}_{2^{\prime },\nu \downarrow }^{†}{c}_{2,\gamma \downarrow }:]$$

$$=$$

$$-{\delta }_{s_4}^{\nu }{\delta }_{\mathbf{q}}^{2^{\prime }}:c_{{\mathbf{p}}^{\prime },s_1\uparrow }^{†}{c}_{\mathbf{p},s_2\uparrow }{c}_{{\mathbf{q}}^{\prime },s_3\downarrow }^{†}{c}_{1^{\prime },\mu \uparrow }^{†}{c}_{2,\gamma \downarrow }:+(f_{s_3-}^{\downarrow \ast }{f}_{\gamma -}^{\downarrow }{\delta }_{{\mathbf{q}}^{\prime }}^2{f}_{s_2+}^{\uparrow }{f}_{\mu +}^{\uparrow \ast }{\delta }_{\mathbf{p}}^{1^{\prime }}-f_{\mu -}^{\uparrow \ast }{f}_{s_2-}^{\uparrow }{\delta }_{\mathbf{p}}^{1^{\prime }}{f}_{\gamma +}^{\downarrow }{f}_{s_3+}^{\downarrow \ast }{\delta }_{{\mathbf{q}}^{\prime }}^2)$$

Full text

Quenching a quantum system across a phase transition, the induced far-from-equilibrium dynamics carries the information of both the initial and final phases. Quantum quench has been extensively applied to study non-equilibrium physics from the real-time dynamics Polkovnikov2011 ; Eisert2015 ; Giamarchi_book ; Heyl2018 . In condensed matter physics, the melting or creation of long-range order can be investigated in the dynamics of symmetry-breaking states, e.g. the survival of magnetic order following an interaction quench in the Hubbard models Tsuji2013 ; Sandri2013 ; Balzer2015 . For topological systems, characterization of topology by quench dynamics has also attracted particular interest very recently Tarnowski2017 ; Song2017 ; Sun2018 ; Wang2017 ; Zhanglin2018 ; Zhanglong2018 ; Yu2018 ; McGinley2018a ; McGinley2018b .

So far the dynamical characterization theories are applicable to noninteracting topological systems Tarnowski2017 ; Sun2018 ; Song2017 ; Wang2017 ; Zhanglin2018 ; Zhanglong2018 ; Yu2018 ; McGinley2018a ; McGinley2018b . For an interacting system, the more challenging but interesting issues could arise. First, the single-particle quantum numbers in correlated systems are not conserved. It is unclear how to define the dynamical topology for the characterization. Second, the interaction can bring about complex effects Polkovnikov2011 , such as dephasing and heating. Their influence on topology remains an open question. Third, symmetry-breaking orders can emerge in correlated systems. The compositive characterization of both topology and symmetry-breaking orders in a correlated quantum quench is, at present, an outstanding issue.

We propose for the first time a dynamical theory to characterize both topology and symmetry-breaking orders, through quenching the spin-$1/2$ Haldane model Haldane1988 ; Jotzu2014 with Hubbard interaction from an initial magnetic phase (Fig. 1a), which exists in strongly interacting regime Hickey2015 ; Zheng2015 ; Liu2016 ; Arun2016 ; Hickey2016 ; Vanhala2016 ; JXWu2016 ; Rubio2018 , to topological regime with relatively weak interaction. We show that the particle-particle interaction has nontrivial correlation effects on the pseudospin dynamics which, after being projected onto the momentum space, follow a novel microscopic Landau-Lifshitz-Gilbert-Bloch equation. With this the dephasing and heating effects are explicitly predicted. Importantly, we find that the correlated quench dynamics exhibit emergent robust topological pattern and universal scaling on the one-dimensional (1D) momentum subspaces called band-inversion surfaces (BISs) Zhanglin2018 ; Zhanglong2018 . These exotic features manifest a deep dynamical bulk-surface correspondence relating both topology and symmetry-breaking orders to correlated quench dynamics on BISs.

*The model.—*The full Hamiltonian of the 2D Haldane-Hubbard model with onsite interaction $U$ reads

[TABLE]

Here $a_{\vec{i}\sigma}$ $(b_{\vec{i}\sigma})$ and $a_{\vec{i}\sigma}^{\dagger}$ ($b_{\vec{i}\sigma}^{\dagger}$) are annihilation and creation operators, respectively, for fermions of spin $\sigma=\uparrow,\downarrow$ on $A$ ($B$) sites. The nearest- ($t_{1}$) and next-nearest-neighbor ($t_{2}$) hopping is considered, with the latter having a phase $\pm\phi$. $M$ is an energy imbalance between $A$ and $B$ sites.

The noninteracting Hamiltonian $H_{0}=\sum_{\bf k,\sigma}{\bf h}({\bf k})\cdot{\bm{\tau}}^{\sigma}$, where ${\bf h}({\bf k})=(h_{x},h_{y},h_{z})$ mimics an effective Zeeman field in Bloch $\mathbb{k}$ space Suppl , with the pseudospin operators $\tau_{z}^{\sigma}=a_{{\bf k}\sigma}^{\dagger}a_{{\bf k}\sigma}-b_{{\bf k}\sigma}^{\dagger}b_{{\bf k}\sigma}$, $\tau_{x}^{\sigma}=a_{{\bf k}\sigma}^{\dagger}b_{{\bf k}\sigma}+{\rm h.c.}$, and $\tau_{y}^{\sigma}=-\mathbbm{i}[\tau_{z}^{\sigma},\tau_{x}^{\sigma}]$. It has been widely studied Hickey2015 ; Zheng2015 ; Liu2016 ; Arun2016 ; Hickey2016 ; Vanhala2016 ; JXWu2016 ; Rubio2018 that in the ground state $|\Psi_{G}\rangle$, an antiferromagnetic (AF) order $m_{\rm AF}=(m_{\downarrow}-m_{\uparrow})/2$ arises for strong repulsive interaction, and the energy imbalance $M$ further leads to a charge order $m_{\rm C}=(m_{\uparrow}+m_{\downarrow})/2$ corrected by Hubbard interaction, characterizing the population difference in the two sublattices. Here the initial orders $m_{\sigma}\equiv\langle a^{\dagger}_{\vec{i}\sigma}a_{\vec{i}\sigma}-b^{\dagger}_{\vec{i}\sigma}b_{\vec{i}\sigma}\rangle U_{\rm in}/4$, with $U_{\rm in}$ the initial strong interaction, and the expectation $\langle\cdot\rangle$ computed in ground state $|\Psi_{G}\rangle$. For the quench study, we can write down initial ground state in the mean-field form $|\Psi_{G}\rangle\rightarrow|\Psi_{\rm MF}\rangle$ which solely depends on the order parameters Suppl , or as the Gutzwiller many-body wave function which is beyond mean field picture S (7, 8, 9), and investigate its evolution under the post-quench Hamiltonian with relatively weak interactions. As studied below and detailed in Supplementary Material Suppl , the central results in this work are valid beyond mean-field regime.

We solve the quench dynamics by the flow equation method Wilson1993 ; Wilson1994 ; Wegner1994 . The process is below. First, through a unitary transformation that changes continuously with a flow parameter $l$, we (nearly) diagonalize the post-quench Hamiltonian at $l\to\infty$ Kehrein_book . Accordingly, the transformation of an operator ${\cal O}(l)$ (including the Hamiltonian) follows the flow equation $d{\cal O}(l)/dl=[\eta(l),{\cal O}(l)]$, where the canonical generator $\eta(l)=[H_{0}(l),H_{I}(l)]=-\eta(l)^{\dagger}$ is anti-Hermitian, with $H_{I}$ the interacting term of the full Hamiltonian. Second, the time-evolved operator ${\cal O}(l\to\infty,t)$ is obtained straightforwardly in the diagonal bases. Finally, we perform the backward transformation so that the operator flows back as ${\cal O}(l\to\infty,t)\rightarrow{\cal O}(0,t)$ Moeckel2008 ; Moeckel2009 . The time evolution is then given in the original bases.

We apply this method to the present system (details are given in supplementary material Suppl ). We consider the ansatz below for post-quench regime

[TABLE]

where ${\cal E}_{\pm}({\bf k})$ are the band energies of $H_{0}$, the normal ordering is with respect to the initial state $|\Psi_{G}\rangle$, and $c_{{\bf k},\pm\sigma}^{\dagger}$ ($c_{{\bf k},\pm\sigma}$) are the creation (annihilation) operators of spin $\sigma=\uparrow\downarrow$ for the upper and lower band states of $H_{0}$ Suppl . The interaction strength $U_{\bf p^{\prime}pq^{\prime}q}^{s_{1}s_{2}s_{3}s_{4}}(l)$ is defined for momentum-conserved scattering channels. With the interaction weaker than the band width, only the leading order contributions from scatterings up to $U^{2}$ will be considered. Taking the previously defined canonical generator $\eta(l)$, the interaction $U(l)$ decays exponentially with $l$ and flows to zero at $l\rightarrow\infty$. We then work out the flow of creation and annihilation operators with respect to $A$ and $B$ sites, with ${\cal A}_{{\bf k}\uparrow}^{\dagger}(l=0)=a_{{\bf k}\uparrow}^{\dagger}$ and ${\cal B}_{{\bf k}\uparrow}^{\dagger}(l=0)=b_{{\bf k}\uparrow}^{\dagger}$, from the same generator $\eta(l)$. Finally we obtain time evolution of pseudospin polarization at momentum $\mathbb{k}$, calculated by $\langle\tau^{\sigma}_{z}({\bf k},t)\rangle=\langle\Psi_{G}|{\cal A}^{\dagger}_{{\bf k}\sigma}(l=0,t){\cal A}_{{\bf k}\sigma}(l=0,t)-{\cal B}^{\dagger}_{{\bf k}\sigma}(l=0,t){\cal B}_{{\bf k}\sigma}(l=0,t)|\Psi_{G}\rangle$, similar for $\langle\tau^{\sigma}_{x,y}({\bf k},t)\rangle$. Note that the single-particle $\mathbb{k}$ is no longer conserved. We project the results onto the single-particle momentum space to study the pseudospin dynamics (Fig. 1b,c).

Equation of motion for pseudospin dynamics.— We show that the essential physics of pseudospin dynamics can be captured by the equation of motion in the projected $\mathbb{k}$-space with ${\bf S}_{\sigma}({\bf k},t)\equiv\frac{1}{2}\langle{\bm{\tau}}^{\sigma}({\bf k},t)\rangle$. Taking into account the leading-order contributions we obtain Suppl

[TABLE]

where the first $2\mathbb{h}$-term corresponds to the single-particle precession, and the second and third terms are induced by correlation effects. The $\eta_{1}^{\sigma}$-term represents the interaction induced damping of precession, and $\eta_{2}^{\sigma}$-term leads to dephasing and heating, with $T_{\rm g}\equiv 1/(2E_{0})$ and $E_{0}({\bf k})=[h^{2}_{x}({\bf k})+h^{2}_{y}({\bf k})+h^{2}_{z}({\bf k})]^{1/2}$. This equation renders a novel mixed microscopic form of Landau-Lifshitz-Gilbert Gilbert2004 and Bloch equations Bloch1946 for magnetization, and is not altered in characterizing the initial phase in mean-field theory or as Gutzwiller state, albeit the beyond-mean-field effects can correct the coefficients of the equation Suppl . The solution reads generically

[TABLE]

where ${\bf S}_{\sigma}^{(0)}({\bf k})=\delta n_{\sigma}({\bf k}){\bf h}({\bf k})/E_{0}({\bf k})$ is the incoherent time-independent part, with $\delta n_{\sigma}({\bf k})=n^{\sigma}_{+-}({\bf k})-n^{\sigma}_{--}({\bf k})$ being the density difference of the initial state populated in the upper ($n^{\sigma}_{+-}$) and lower ($n^{\sigma}_{--}$) eigen-bands, ${\bf S}_{\sigma}^{(c)}({\bf k},t)\sim\cos(t/T_{\rm g})$ is the single-particle coherent oscillation, ${\bf S}_{\sigma}^{(h)}({\bf k},t)\approx-\lambda^{\sigma}_{1}({\bf k},t){\bf S}_{\sigma}^{(c)}({\bf k},t)$ ($\lambda^{\sigma}_{1}\propto U^{2}/E_{0}^{2}$) represents the interaction-induced high-frequency fluctuation, which reduces the single-particle procession, and ${\bf S}_{\sigma}^{(l)}({\bf k},t)\approx-2\lambda^{\sigma}_{2}({\bf k},t){\bf h}({\bf k})/E_{0}({\bf k})$ ($\lambda^{\sigma}_{2}\propto U^{2}/E_{0}^{2}$) denotes the low-frequency interaction effect, which equilibrates the density distribution on upper and lower bands (heating). The coefficients $\lambda_{1,2}^{\sigma}$ are related to the factors $\eta_{1,2}^{\sigma}$ (see later). Note that the entire many-body system evolves unitary. The dephasing and heating arise in the projected quench dynamics at fixed $\mathbb{k}$, since all the particles with $\mathbb{k}^{\prime}\neq\mathbb{k}$ act as a bath scattering the $\mathbb{k}$ state.

The $\eta_{1,2}^{\sigma}$ terms are momentum dependent. For comparison, we first define the BIS for single-particle Hamiltonian $H_{0}$ Zhanglin2018 ; Zhanglong2018 ; Sun2018 , being the momentum subspace where band inversion occurs and time-averaged spin polarizations $\overline{\mathbb{S}_{\sigma}(\mathbb{k},t)}|_{U=0}=0$, equivalent to ${\bf S}_{\sigma}^{(0)}({\bf k})=0$. On the single-particle BIS, we have $\eta_{1}^{\sigma}\simeq-4(d\lambda^{\sigma}_{2}/dt)T_{\rm g}$ and $\eta_{2}^{\sigma}\simeq 4(d\lambda^{\sigma}_{1}/dt)T_{\rm g}n^{\sigma}_{+-}n^{\sigma}_{--}$, where $d\lambda^{\sigma}_{1,2}/dt$ are approximately constant in early time Suppl . Fig. 2c shows that near the BIS (dashed line) $\eta_{1}^{\sigma}$ is small (due to the cancelling of the two-band contributions) and the heating due to $\eta_{2}^{\sigma}$-term dominates the correlation effect. In comparison, the damping enhances at $\mathbb{k}$ away from BIS. The heating shortens the pseudospin vector while the damping drags the vector towards the $\mathbb{h}$ axis (see Fig. 2a-b).

*Topology emerging on BIS.–*We show now the correlated pseudospin dynamics on BISs exhibiting emergent topological pattern, which corresponds to the post-quench topology. From Eq. (3) one finds that the damping $\eta_{1}$-term modifies the procession. Thus the BISs in the presence of interactions, with $\overline{\mathbb{S}_{\sigma}(\mathbb{k},t)}|_{U\neq 0}=0$, is deformed from the single-particle BISs where $\mathbb{h}$ is perpendicular to ${\mathbb{S}}_{\sigma}$. In contrast, one can prove that the positions $\mathbb{k}=\mathbb{k}_{\rho}$ of topological charges, with $\mathbb{h}_{\rm so}(\mathbb{k}_{\rho})\equiv(h_{y},h_{x})=0$ in the noninteracting regime, is unchanged from the equation of motion (3). As shown below, the emergent topology of quench dynamics on BISs reflects the total topological charges enclosed by the BISs.

To characterize the topology emerging on the BISs, we introduce a dynamical field ${\bf g}^{\sigma}({\bf k})$, with the components $\mathbb{g}^{\sigma}(\mathbb{k})=\pm\frac{1}{{\cal N}_{k}}\partial_{k_{\perp}}\overline{\mathbb{S}_{\sigma}(\mathbb{k},t)}$. It takes $+$ (or $-$) for $\sigma=\uparrow$ (or $\downarrow$), the momentum $k_{\perp}$ is perpendicular to the BIS, and ${\cal N}_{k}$ is the normalization factor. Due to the damping and heating effects, the ${\bf g}^{\sigma}({\bf k})$ vector is generally not in the $x$-$y$ plane. To characterize the topology, we project the dynamical field onto the $x$-$y$ plane such that ${\bf g}_{\parallel}^{\sigma}({\bf k})=\hat{e}_{\parallel}\cdot{\bf g}^{\sigma}({\bf k})=(g^{\sigma}_{y},g^{\sigma}_{x})$, and can prove that ${\bf g}^{\sigma}_{\parallel}({\bf k})\simeq{\bf h}_{\rm so}({\bf k})$ on the interacting BISs Suppl . Therefore, the winding of ${\bf g}^{\sigma}_{\parallel}({\bf k})$ on BISs quantifies the total topological charges (at zeros of $\mathbb{h}_{\rm so}$-vector) enclosed by the BISs, corresponding to the topology of the post-quench regime and valid for the present interacting system. This new characterization is different from the free-fermion regime, where the topology emerges in the bare dynamical field $\mathbb{g}^{\sigma}(\mathbb{k})$ Zhanglin2018 ; Zhanglong2018 , not directly applicable to the regime with interactions. A typical example for the topology of the dynamical field is illustrated in Fig. 3d,h.

*Magnetic order from quench dynamics on BIS.–*The AF and charge orders are closely related to the spin and density distributions in $A$ and $B$ sites, hence related to the pseudospin dynamics, in which the BISs also play the pivotal role. For the initial phase characterized by the mean-field theory that $|\Psi_{G}\rangle\rightarrow|\Psi_{\rm MF}\rangle$, the BIS defined by $\overline{{\bf S}_{\sigma}(\mathbb{k},t)}=0$ is alternatively interpreted as the momenta satisfying $E_{0}^{2}({\bf k})+m_{\sigma}h_{z}({\bf k})=-(d\lambda^{\sigma}_{2}/dt)TE_{0}({\bf k})E_{0}^{\sigma}({\bf k})$ with $E_{0}^{\sigma}\equiv\sqrt{E_{0}^{2}+2m_{\sigma}h_{z}+m_{\sigma}^{2}}$. Here $T$ is the interval for time averaging and the right-hand side represents shift of BISs by interaction. This formula shows that BISs depends on both the pre-quench phase ($m_{\sigma}$) and the post-quench Hamiltonian. Further, the half amplitude, defined as $Z_{0}^{\sigma}({\bf k})\equiv\langle\tau_{z}^{\sigma}({\bf k},t=0)\rangle$, on BISs reads $Z^{\sigma}_{0}=(d\lambda^{\sigma}_{2}/dt)Th_{z}/E_{0}-m_{\sigma}(E_{0}^{2}-h_{z}^{2})/(E_{0}^{2}E_{0}^{\sigma})$, which also relates the magnetization to quench dynamics. With these results and up to the leading order correction from interaction, we show the scaling Suppl

[TABLE]

where $f(m_{\sigma})=m_{\sigma}T_{0}$ and $g(Z^{\sigma})=\sqrt{1-Z^{\sigma\,2}_{0}}/\pi$, with $T_{0}({\bf k})=\pi/E_{0}$. The result in Eq. (5) gives a universal scaling at any $\mathbb{k}$ on BISs, insensitive to interactions.

We provide numerical results in Fig. 4a for spin-up component. By identifying the BIS (the dashed purple curve), we record the short-time dynamics at momenta of three kinds: inside (${\bf k}^{<}_{1}$), outside (${\bf k}^{>}_{1}$), and right on the BIS (${\bf k}^{=}_{1,2,3}$). We measure both $Z_{0}(\mathbb{k})$ and $T_{0}(\mathbb{k})$ versus order parameters $m_{\uparrow}/t_{1}$ (Fig. 4b). The results are plotted as points $(|m_{\uparrow}T_{0}|,\sqrt{1-Z_{0}^{2}})$ in Fig. 4c, showing that the data measured on the BIS all satisfy the scaling (5). In experiment, one can obtain $m_{\sigma}$ by measuring only the first one or two oscillations. The AF order is then obtained by $m_{\rm AF}=(m_{\downarrow}-m_{\uparrow})/2$, and the charge order is $m_{\rm C}=(m_{\uparrow}+m_{\downarrow})/2$. Finally, we emphasize that the scaling (5) is satisfied beyond the mean-field theory. For the initial phase described with the correlated Gutzwiller wave function $|\Psi_{G}\rangle$, the same scaling holds, with $Z_{0}^{\sigma}$ and the order parameters $m_{\sigma}$ renormalized by correlations in the more precise Gutzwiller ground state. For simplicity we put the details in the Supplementary Material Suppl .

*Conclusion.–*We have proposed a dynamical theory to characterize both the topology and symmetry-breaking orders by quantum quench in correlated topological system. By quenching the Haldane-Hubbard model from initial symmetry-breaking phase into topological regime, the induced quantum dynamics on band inversion surfaces (BISs) exhibits emergent topology and universal scaling, which uniquely correspond to and thus characterize the post-quench equilibrium topological state and pre-quench symmetry-breaking orders, respectively. Our results are shown to be valid beyond mean-field theory Suppl , hence reveal a deep dynamical bulk-surface correspondence for topology and symmetry-breaking orders. These results are expected to be generically seen in Chern-Hubbard insulators and 1D topological-Hubbard systems. Note that the pseudospin dynamics can be measured by the tomography of Bloch states Hauke2014 ; Flaschner2016 . This work opens an avenue to explore profound correlation physics with novel topology by quench dynamics.

This work was supported by the National Key R&D Program of China (2016YFA0301604, 2017YFA0304203), National Natural Science Foundation of China (11574008, 11761161003, 11825401, and 11874038), and the Strategic Priority Research Program of Chinese Academy of Science (Grant No. XDB28000000). Y.H. also acknowledges support from the National Thousand-Young-Talents Program, and Changjiang Scholars and Innovative Research Team (Grant No. IRT13076).

SUPPLEMENTAL MATERIAL

In the Supplementary Materials, we first provide the details in deriving the flow equation, dynamical topology, and universal scaling on BISs when the initial ground state is characterized with mean-field theory. Then, we show in detail that these results can be further obtained beyond the mean-field theory.

I I. Hamiltonians and the mean-field ground state

The Bloch Hamiltonian of the noninteracting Haldane model regardless of the spin can be written as

[TABLE]

with $h_{x}({\bf k})=-t_{1}\sum_{j}\cos({\bf k}\cdot{\bf e}_{j})$, $h_{y}({\bf k})=-t_{1}\sum_{j}\sin({\bf k}\cdot{\bf e}_{j})$ and $h_{z}({\bf k})=M-2t_{2}\sin\phi\sum_{j}\sin({\bf k}\cdot{\bf v}_{j})$. Here we have removed the trivial identity matrix term, with ${\bf e}_{1}=(0,a_{0})$, ${\bf e}_{2}=(-\frac{\sqrt{3}a_{0}}{2},-\frac{a_{0}}{2})$, ${\bf e}_{3}=(\frac{\sqrt{3}a_{0}}{2},-\frac{a_{0}}{2})$ and ${\bf v}_{1}=(\sqrt{3}a_{0},0)$, ${\bf v}_{2}=(-\frac{\sqrt{3}a_{0}}{2},\frac{3a_{0}}{2})$, ${\bf v}_{3}=-{\bf v}_{2}-{\bf v}_{1}$ ($a_{0}$ is the lattice constant). Moreover, we set the energy difference between the two sublattices $M/t_{1}=-0.5$ if it is considered. Thus, with $\phi=\pi/2$ and $t_{2}=0.3t_{1}$, the noninteracting system lies in the topological phase with Chern number $C=1$ S (1). The two energy bands read

[TABLE]

We write $H_{0}=\sum_{{\bf k},s=\pm\sigma}{\cal E}_{s}({\bf k})c_{{\bf k},s\sigma}^{\dagger}c_{{\bf k},s\sigma}$, with

[TABLE]

One can easily obtain $\chi_{+}=\sqrt{\frac{1}{2E_{0}(E_{0}-h_{z})}}(h_{x}-\mathbbm{i}h_{y})$, $\chi_{-}=-\sqrt{\frac{1}{2E_{0}(E_{0}+h_{z})}}(h_{x}-\mathbbm{i}h_{y})$, $\xi_{+}=\sqrt{\frac{E_{0}-h_{z}}{2E_{0}}}$ and $\xi_{-}=\sqrt{\frac{E_{0}+h_{z}}{2E_{0}}}$. We further have

[TABLE]

and, similarly,

[TABLE]

Thus, the on-site interaction

[TABLE]

where $U_{\bf p^{\prime}pq^{\prime}q}^{s_{1}s_{2}s_{3}s_{4}}\equiv U\delta^{\bf p+q}_{\bf p^{\prime}+q^{\prime}}\Lambda_{\bf p^{\prime}pq^{\prime}q}^{s_{1}s_{2}s_{3}s_{4}}$ and

[TABLE]

For large $U$, we consider the symmetry-breaking order in $z$ direction, and write the Hubbard interaction in the mean-field form

[TABLE]

Here $\langle\cdot\rangle$ is taken with respect to the mean-field ground state of the total Hamiltonian, which can be solved self-consistently. We define the antiferromagnetic (AF) order $m_{\rm AF}\equiv\langle b^{\dagger}_{i\uparrow}b_{i\uparrow}-b^{\dagger}_{i\downarrow}b_{i\downarrow}\rangle U/2=-\langle a^{\dagger}_{i\uparrow}a_{i\uparrow}-a^{\dagger}_{i\downarrow}a_{i\downarrow}\rangle U/2$ and the charge order $m_{\rm C}=\langle a^{\dagger}_{i\uparrow}a_{i\uparrow}+a^{\dagger}_{i\downarrow}a_{i\downarrow}-b^{\dagger}_{i\uparrow}b_{i\uparrow}-b^{\dagger}_{i\downarrow}b_{i\downarrow}\rangle U/4$. After the Fourier transform, we have the Bloch Hamiltonian

[TABLE]

where the magnetic order $m_{\uparrow}=m_{\rm C}-m_{\rm AF}$ and $m_{\downarrow}=m_{\rm C}+m_{\rm AF}$. Regarding the orders $m_{\rm C}$ and $m_{\rm AF}$ as the input parameters, the Hamiltonian can be diagonalized as ${\cal H}_{\rm MF}({\bf k})=\sum_{s=\pm,\sigma=\uparrow\downarrow}\bar{{\cal E}}_{s\sigma}({\bf k})\bar{c}_{{\bf k},s\sigma}^{\,\dagger}\bar{c}_{{\bf k},s\sigma}$, where $\bar{{\cal E}}_{\pm\sigma}=\pm E_{0}^{\sigma}$, with $E_{0}^{\sigma}\equiv\sqrt{E_{0}^{2}+2m_{\sigma}h_{z}+m_{\sigma}^{2}}$. One can find the relation between the noninteracting and mean-field solutions ($\sigma=\uparrow\downarrow$):

[TABLE]

with

[TABLE]

The AF state at half-filling can be denoted as $|\Psi_{\rm MF}\rangle=\prod_{{\bf k}}\bar{c}_{{\bf k},-\uparrow}^{\,\dagger}\bar{c}_{{\bf k},-\downarrow}^{\,\dagger}|0\rangle$, and $|\Psi_{\rm MF}\rangle\simeq\prod_{{\bf k}}a_{{\bf k}\uparrow}^{\,\dagger}b_{{\bf k}\downarrow}^{\,\dagger}|0\rangle$ when $m_{\rm AF}\to\infty$.

II II. Flow equations

We study the interacting Haldane model by flow equation method. The expansion parameter is the (small) interaction $U$ and normal ordering is with respect to the AF state $|\Psi_{\rm MF}\rangle$, with

[TABLE]

We start with the ansatz

[TABLE]

where the interaction $U_{\bf p^{\prime}pq^{\prime}q}^{s_{1}s_{2}s_{3}s_{4}}(l)$ is responsible for the flow of the Hamiltonian and the flow of band energies and higher order terms are neglected. Since

[TABLE]

we have the generator

[TABLE]

where $\Delta_{\bf p^{\prime}pq^{\prime}q}^{s_{1}s_{2}s_{3}s_{4}}\overset{\rm def}{=}{\cal E}_{s_{1}}({\bf p}^{\prime})-{\cal E}_{s_{2}}({\bf p})+{\cal E}_{s_{3}}({\bf q}^{\prime})-{\cal E}_{s_{4}}({\bf q})$ is the energy difference before and after scattering. Since

[TABLE]

the flow of the interaction is given by

[TABLE]

which decays to zero when the flow parameter $l\to\infty$.

Next we work out the flow equation transformation for the creation operators. Since $a^{\dagger}_{{\bf k}\sigma}=\chi^{*}_{+}c^{\dagger}_{{\bf k},+\sigma}+\chi^{*}_{-}c^{\dagger}_{{\bf k},-\sigma}$, $b^{\dagger}_{{\bf k}\sigma}=\xi^{*}_{+}c^{\dagger}_{{\bf k},+\sigma}+\xi^{*}_{-}c_{{\bf k},-\sigma}$ and the relations

[TABLE]

we assume

[TABLE]

Here $h_{{\bf k},+}(l=0)=g_{{\bf k},+}(l=0)=\chi^{*}_{+}({\bf k})$, $h_{{\bf k},-}(l=0)=g_{{\bf k},-}(l=0)=\chi^{*}_{-}({\bf k})$, and $M_{{\bf k},{\bf p^{\prime}q^{\prime}q}}^{\mu\nu\gamma}(l=0)=W_{{\bf k},{\bf p^{\prime}pq^{\prime}}}^{\mu\nu\gamma}(l=0)=0$. The operators ${\cal B}_{{\bf k}\sigma}^{\dagger}(l)$ take the same form as in Eq. (II) but with $h_{{\bf k},+}(l=0)=g_{{\bf k},+}(l=0)=\xi^{*}_{+}({\bf k})$ and $h_{{\bf k},-}(l=0)=g_{{\bf k},-}(l=0)=\xi^{*}_{-}({\bf k})$. With

[TABLE]

and

[TABLE]

we obtain the leading-order flow equations for the creation operators

[TABLE]

and

[TABLE]

where $F_{\bf p^{\prime}qq^{\prime}}^{s_{2}s_{3}s_{4},\mu\gamma\nu}=f_{s_{2}+}^{\uparrow}({\bf p}^{\prime})f_{\mu+}^{\uparrow\,*}({\bf p}^{\prime})f_{s_{3}-}^{\downarrow\,*}({\bf q})f_{\gamma-}^{\downarrow}({\bf q})f_{s_{4}+}^{\downarrow}({\bf q}^{\prime})f_{\nu+}^{\downarrow\,*}({\bf q}^{\prime})+f_{s_{2}-}^{\uparrow}({\bf p}^{\prime})f_{\mu-}^{\uparrow\,*}({\bf p}^{\prime})f_{s_{3}+}^{\downarrow\,*}({\bf q})f_{\gamma+}^{\downarrow}({\bf q})f_{s_{4}-}^{\downarrow}({\bf q}^{\prime})f_{\nu-}^{\downarrow\,*}({\bf q}^{\prime})$ and $G_{\bf pp^{\prime}q^{\prime}}^{s_{1}s_{2}s_{4},\nu\mu\gamma}=f_{s_{1}-}^{\uparrow\,*}({\bf p})f_{\nu-}^{\uparrow}({\bf p})f_{s_{2}+}^{\uparrow}({\bf p}^{\prime})f_{\mu+}^{\uparrow\,*}({\bf p}^{\prime})f_{s_{4}+}^{\downarrow}({\bf q}^{\prime})f_{\gamma+}^{\downarrow\,*}({\bf q}^{\prime})+f_{s_{1}+}^{\uparrow\,*}({\bf p})f_{\nu+}^{\uparrow}({\bf p})f_{s_{2}-}^{\uparrow}({\bf p}^{\prime})f_{\mu-}^{\uparrow\,*}({\bf p}^{\prime})f_{s_{4}-}^{\downarrow}({\bf q}^{\prime})f_{\gamma-}^{\downarrow\,*}({\bf q}^{\prime})$.

We adopt the forward-backward transformation S (2, 3) to calculate the time-evolved operators. The forward (or backward) transformations are derived by integrating the flow equations (II) and (II) from $l=0$ to $\infty$ (or from $l=\infty$ to 0) with different initial conditions. We keep the terms up to second order in $U$ and obtain the approximate analytic solutions. Take the number operator ${\cal N}^{A}_{{\bf k}\uparrow}(l)={\cal A}_{{\bf k}\uparrow}^{\dagger}(l){\cal A}_{{\bf k}\uparrow}(l)$ as an example. Time evolution yields

[TABLE]

Since

[TABLE]

we obtain the distribution of spin-up particles at $A$ sites

[TABLE]

The computation of $h_{{\bf k},\pm}(l=0,t)$ and $M_{{\bf k},{\bf p^{\prime}q^{\prime}q}}^{\mu\nu\gamma}(l=0,t)$ is achieved by composing the forward transformation (FT), the time evolution (TE) and the backward transformation (BT), such as

[TABLE]

Up to now, the analytic solutions are very complicated despite the neglect of higher order terms. It is mainly due to the various possible scattering channels in the flow equations (II) and (II). To simplify the analysis, we take into account only the major contribution, i.e.

[TABLE]

where ${\cal F}_{\bf p^{\prime}qq^{\prime}}^{\mu\gamma\nu}\overset{\rm def}{=}n_{\mu+}^{\uparrow}({\bf p}^{\prime})n_{\gamma-}^{\downarrow}({\bf q})n_{\nu+}^{\downarrow}({\bf q}^{\prime})+n_{\mu-}^{\uparrow}({\bf p}^{\prime})n_{\gamma+}^{\downarrow}({\bf q})n_{\nu-}^{\downarrow}({\bf q}^{\prime})$ and ${\cal G}_{\bf pp^{\prime}q^{\prime}}^{\nu\mu\gamma}\overset{\rm def}{=}n_{\nu-}^{\uparrow}({\bf p})n_{\mu+}^{\uparrow}({\bf p}^{\prime})n_{\gamma+}^{\downarrow}({\bf q}^{\prime})+n_{\nu+}^{\uparrow}({\bf p})n_{\mu-}^{\uparrow}({\bf p}^{\prime})n_{\gamma-}^{\downarrow}({\bf q}^{\prime})$ with $n_{s_{1}s_{2}}^{\sigma}({\bf k})\equiv|f_{s_{1}s_{2}}^{\sigma}({\bf k})|^{2}$.

III III. Pseudospin dynamics

For convenience’s sake, we denote ${\cal E}_{+}({\bf k})-{\cal E}_{-}({\bf k})=2E_{0}({\bf k})\equiv 1/T_{\rm g}({\bf k})$, $X_{\sigma}({\bf k},t)\equiv\langle\tau_{x}^{\sigma}({\bf k},t)\rangle$, $Y_{\sigma}({\bf k},t)\equiv\langle\tau_{y}^{\sigma}({\bf k},t)\rangle$ and $Z_{\sigma}({\bf k},t)\equiv\langle\tau_{z}^{\sigma}({\bf k},t)\rangle$ in the following. By the forward-backward transformation, we have the results

[TABLE]

where the incoherent part

[TABLE]

the coherent time-dependent oscillation

[TABLE]

the interaction-induced high-frequency fluctuation

[TABLE]

with

[TABLE]

and the interaction-induced low-frequency fluctuation

[TABLE]

with

[TABLE]

Here we have denoted ($s=\pm$)

[TABLE]

The high-frequency fluctuations come from the scattering processes from the upper to the lower band or the other way round via the background. Note that the expressions in Eqs. (III) and (III) resemble the structure of transition probability in time-dependent perturbation theory (see, e.g. Ref. S (4)), and the sinusoidal time dependence $\sin^{2}(\Delta_{\bf pp^{\prime}qq^{\prime}}^{s_{1}s_{2}s_{3}s_{4}}t/2)/(\Delta_{\bf pp^{\prime}qq^{\prime}}^{s_{1}s_{2}s_{3}s_{4}})^{2}$ determines the contribution of each scattering process in the time evolution. One can find that the dependence $\sin^{2}(\omega t/2)/\omega^{2}$, as a function of $\omega$ for fixed $t$, has a major peak with the height $\propto t^{2}$ and the width $\propto 1/t$ S (4). Hence, after a summation, the parameters $\lambda^{\sigma}_{1,2}(t)$ are approximately linear in $t$, i.e., $\lambda^{\sigma}_{1,2}(t)\propto t$ (see Fig. S1).

We define the pseudospin vector

[TABLE]

and assume the equation of motion takes the form

[TABLE]

where the first term in the right-hand side corresponds to the precessional motion, $\eta_{1}^{\sigma}$-term represents the damping effect, and $\eta_{2}^{\sigma}$-term describes the heating. According to Eqs. (III-III), we obtain

[TABLE]

Note that $\eta^{\sigma}_{1,2}({\bf k})$ are approximated as time-independent in short time due to the linear time dependence of $\lambda^{\sigma}_{1,2}$ (Fig. S1). In Fig. S2, we show the calculated results of $\eta_{1,2}^{\sigma}$ for $M=0$, $m_{\rm AF}\to\infty$, which have a symmetrical (for $\eta_{2}^{\sigma}$) or antisymmetrical (for $\eta_{1}^{\sigma}$) distribution.

Finally, we discuss the reliability of our method. It should be pointed out that secular terms may arise from our zeroth order approximation for the time evolution, e.g. in Eq. (II), we take $H(l\to\infty)\approx H_{0}$. When $\Delta_{\bf p^{\prime}pq^{\prime}q}^{s_{1}s_{2}s_{3}s_{4}}=0$, the canonical generator (S17) vanishes. Hence the energy-diagonal contributions of $H_{I}$ cannot be erased by flow equations. The perturbation solutions would fail on long-time scales. This failure can be also indicated by the sinusoidal time dependence in Eqs. (III) and (III). For a large $t$, the function $\sin^{2}(\omega t/2)/\omega^{2}$ has a very narrow peak, and approximate energy conservation is required, which means the energy-diagonal contributions can not be neglected for a long time evolution. Fortunately, we only need to focus on short-time pseudospin dynamics, from which the topology as well as magnetic order can be measured. Furthermore, from the early stage dynamics, we can qualitatively analyze which interaction effect dominates even for a relatively long time.

IV IV. Detecting the topology

In Ref. S (5), we have developed a dynamical classification theory, which is applicable to noninteracting topological systems and to the situation that the quench starts from a deep trivial regime. Here we first generalize the theory to the shallow quench case, which corresponds to initializing finite magnetization in the interaction quench, and then discuss its feasibility in the current interacting system.

IV.0.1 A. Projection approach

In this subsection, we consider the noninteracting system and generalize the dynamical classification theory in Ref. S (5) to the situation that the pre-quench state is not completely polarized. For a post-quench Hamitonian ${\cal H}({\bf k})={\bf h}({\bf k})\cdot{\bm{\tau}}$, the spin texture reads ($i=x,y,z$)

[TABLE]

where $\rho_{0}$ is the density matrix of the initial state. Although it is defined for an infinite time period, the time average can be taken over several oscillations for each ${\bf k}$, and the results are unchanged. The band inversion surfaces (BISs) are defined as

[TABLE]

This implies that on the (noninteracting) BIS, the spin vector ${\bf S}({\bf k})\equiv\frac{1}{2}\left(\langle\tau_{x}\rangle,\langle\tau_{y}\rangle,\langle\tau_{z}\rangle\right)$ is perpendicular to the field ${\bf h}({\bf k})$, i.e., ${\bf S}({\bf k})\cdot{\bf h}({\bf k})=0$. We denote by $k_{\perp}$ the direction perpendicular to the contour ${\rm Tr}[\rho_{0}{\cal H}]$. For the contours infinitely close to the BIS, we have ${\rm Tr}[\rho_{0}{\cal H}]\simeq\pm k_{\perp}$ and the variation of $h_{i}$ is of order $\mathcal{O}(k_{\perp})$. Therefore, the directional derivative on the BIS reads

[TABLE]

Without loss of generality, we consider quenching ${h}_{z}$. When the initial state $\rho_{0}$ is fully polarized ($|h_{z}|\to\infty$ for $t<0$), the BIS conincides with the surfaces with $h_{z}({\bf k})=0$, and $\partial_{k_{\perp}}\overline{\langle\tau_{z}\rangle}$ vanishes. Hence $\partial_{k_{\perp}}\overline{\langle{\bm{\tau}}\rangle}$ is a vector in the $x$-$y$ plane, and the bulk topology is well defined by the winding of the spin-orbit (SO) field ${\bf h}_{\rm so}\equiv(h_{y},h_{x})$ along BISs, which is characterized by the dynamical field $\partial_{k_{\perp}}\overline{\langle{\bm{\tau}}\rangle}=(\partial_{k_{\perp}}\overline{\langle\tau_{y}\rangle},\partial_{k_{\perp}}\overline{\langle\tau_{x}\rangle})$ S (5). From the viewpoint of topological charges, which are located at ${\bf h}_{\rm so}({\bf k})=0$, the winding of $\partial_{k_{\perp}}\overline{\langle{\bm{\tau}}\rangle}$ counts the total charges enclosed by BISs S (6). The BIS where ${\bf S}({\bf k})\cdot{\bf h}({\bf k})=0$ divides the charges into two categories: ${\bf S}({\bf k})\cdot{\bf h}({\bf k})>0$ and ${\bf S}({\bf k})\cdot{\bf h}({\bf k})<0$. The winding of $\partial_{k_{\perp}}\overline{\langle{\bm{\tau}}\rangle}$ in fact characterizes the charges of the same category.

Now we consider the case that the initial state is not completely polarized, i.e., at $t=0$, $\langle\tau_{z}({\bf k})\rangle=1$ (or $-1$) does not hold for all ${\bf k}$ but $\langle\tau_{z}({\bf k})\rangle>0$ (or $<0$) does. In this case, the topological charges enclosed by BISs are unchanged. The reason is as follows: First, by definition their locations are irrelevant to the initial state. Second, the category, i.e., the condition ${\bf S}({\bf k})\cdot{\bf h}({\bf k})=\frac{1}{2}h_{z}({\bf k})\langle\tau_{z}({\bf k})\rangle>0$ or $<0$, remains the same.The BIS encloses the same charges as in the completely polarized case. Note that the topological charges are characterized by the winding of ${\bf h}_{\rm so}$. Although the vector $\partial_{k_{\perp}}\overline{\langle{\bm{\tau}}\rangle}$ is not in the $x$-$y$ plane in general, we can define the topological invariant by the winding of a projected dynamical field $(\partial_{k_{\perp}}\overline{\langle\tau_{y}\rangle},\partial_{k_{\perp}}\overline{\langle\tau_{x}\rangle})$. The dynamical field defined in the completely polarized case can be also regarded as a projection of $\partial_{k_{\perp}}\overline{\langle{\bm{\tau}}\rangle}$ but with $\partial_{k_{\perp}}\overline{\langle\tau_{z}\rangle}=0$.

IV.0.2 B. Dynamical characterization in interacting systems

In this subsection, we will show that the dynamical characterization theory discussed above is also applicable to the interacting Haldane model. According to the results shown in Eqs. (III-S36), the time-averaged pseudospin textures in the presence of interaction are ($i=x,y,z$)

[TABLE]

where $T$ is the period over which the time average is taken and $\lambda^{\sigma}_{2}\propto U^{2}$ defined in Eq. (III) represents the interaction shift. Thus, the (interacting) BIS is determined by $\overline{\langle\tau^{\sigma}_{i}({\bf k})\rangle}=0$, which leads to

[TABLE]

Note that in the interacting system, $k_{\perp}$ is defined to be perpendicular to the contour of $\delta n^{\sigma}_{I}({\bf k})$. For the contours infinitely close to $\delta n^{\sigma}_{I}({\bf k}_{0})=0$, we have $\delta n^{\sigma}_{I}({\bf k}_{0}\pm\hat{e}_{\perp}k_{\perp})\simeq\pm c_{I}k_{\perp}/E_{0}$, with $c_{I}$ being a coefficient dependent on $m_{\sigma}$, $U$ and $T$. Therefore, we have

[TABLE]

which means the emergent gradient field $\partial_{k_{\perp}}\overline{\langle{\bm{\tau}}^{\sigma}\rangle}$ on the (interacting) BIS still characterizes the vector field ${\bf h}({\bf k})$ despite of the interaction effect. Due to the AF order, quenches for the two spins $\sigma=\uparrow\downarrow$ are along opposite directions. Thus, according to Ref. S (5), we define the projected dynamical fields on the BIS ${\bf g}_{\parallel}^{\sigma}({\bf k})=(g^{\sigma}_{y},g^{\sigma}_{x})$ with components given by

[TABLE]

Here the sign $+$ (or $-$) is for $\sigma=\uparrow$ (or $\downarrow$) and ${\cal N}_{k}$ is the normalization factor. The topological invariant is then defined by the winding of the projected dynamical field with $\sigma=\uparrow$ or $\downarrow$:

[TABLE]

Here a special case is shown in Fig. S3 with $M=0$ and $m_{\rm AF}\to\infty$. In this case, one can see that the damping factor $\eta^{\sigma}_{1}=-4T_{\rm g}d\lambda^{\sigma}_{2}/dt$ vanishes right on the noninteracting BIS where $n_{--}({\bf k})=n_{+-}({\bf k})$ (Fig. S2), which is due to the exact cancelling of the two contributions in Eq (III). Thus, the BIS does not move in the presence of interaction. Moreover, the distributions of $\eta^{\sigma}_{1}$ are antisymmetrical. As shown in Fig. S3, the time-averaged textures $\overline{\langle\tau_{i}^{\sigma}({\bf k})\rangle}$ take the same distributions as in the noninteracing case, except for a small reduction of polarization values. The time averages are taken over 10 times of oscillation period for each ${\bf k}$. An example of a general case with finite $m_{\rm AF}$ is discussed in Fig. 3 of the main text.

V V. Measuring the magnetic order

We aim to obtain the magnetic order $m_{\sigma}$ by measuring the pseudospin dynamics. In the presence of interaction, the BIS is given by Eq. (S46). Here we assume $|d\lambda^{\sigma}_{2}/dt|T\ll 1$. Since

[TABLE]

where $E_{0}^{\sigma}\equiv\sqrt{E_{0}^{2}+2m_{\sigma}h_{z}+m_{\sigma}^{2}}$, the BIS can be alternatively interpreted as the momenta satisfying

[TABLE]

Note that when $h_{x}=h_{y}=0$, the above equation becomes $(1+Td\lambda^{\sigma}_{2}/dt)(h_{z}^{2}+m_{\sigma}h_{z})=0$, which fails to hold for $|m_{\sigma}|>|h_{z}|$ except that $h_{z}=0$. That is to say when we consider the quench from a trivial phase ($|m_{\sigma}|>|h_{z}|$), the BIS would not move across a charge where ${\bf h}_{\rm so}=0$ except that topological phase transition occurs. Furthermore, half of the amplitude in the early time $Z_{0}^{\sigma}({\bf k})\overset{\rm def}{=}\langle\tau_{z}^{\sigma}({\bf k},t=0)\rangle$ reads

[TABLE]

on the BIS. Equations (S51) and (S69) provide two relations for the derivation of the magnetic order $m_{\sigma}$. We regard the interaction effect as a perturbation and approximate $m_{\sigma}$ and $h_{z}$ to the first order of $\varepsilon\equiv Td\lambda^{\sigma}_{2}/dt$, i.e. $m_{\sigma}=m_{\sigma}^{(0)}+\varepsilon m_{\sigma}^{(1)}$ and $h_{z}=h_{z}^{(0)}+\varepsilon h_{z}^{(1)}$. We then have

[TABLE]

From Eqs. (S53a) and (S53c), we obtain

[TABLE]

Substituting the results into Eqs. (S53b) and (S53d) leads to

[TABLE]

Finally, to the second order of $U$, we have

[TABLE]

which is the universal scaling behavior immune to the interaction. The AF and charge orders are finally given by $m_{\rm AF}=(m_{\downarrow}-m_{\uparrow})/2$ and $m_{\rm C}=(m_{\uparrow}+m_{\downarrow})/2$.

VI VI. Results beyond mean field theory

In this section, we examine the beyond-mean-field (BMF) effect of the initial state. To this end, we take a Gutzwiller wave function S (7), instead of a mean-field ground state, to describe the AF order phase, and examine the corrections to our results.

VI.0.1 A. Gutzwiller wave function

We take the Gutzwiller ansatz (see, e.g., Refs S (8, 9) and references therein)

[TABLE]

where $\alpha$ is the variational parameter with $0\leq\alpha\leq 1$. Here we have denoted $A_{i}\equiv a^{\dagger}_{i\uparrow}a^{\dagger}_{i\downarrow}a_{i\downarrow}a_{i\uparrow}$ and $B_{i}\equiv b^{\dagger}_{i\uparrow}b^{\dagger}_{i\downarrow}b_{i\downarrow}b_{i\uparrow}$, and $|\Psi_{\rm MF}\rangle$ is the mean-field ground state. The $\alpha$-terms suppress the double occupation of the mean-field ground state. To explore the correlation effect, we keep the leading-order terms and obtain that

[TABLE]

After transforming into the momentum space, we have the unnormalized wave function

[TABLE]

with $\Lambda_{\bf p^{\prime}pq^{\prime}q}^{s_{1}s_{2}s_{3}s_{4}}$ given by Eq. (S7). One can see that the constructed wave function is in fact a superposition of the mean-field state and a series of excited states via single scattering. Thus the correlation effects are explicitly incorporated in the Gutzwiller mang-body ground state. The initial state $|\Psi_{\rm MF}\rangle$ is updated by diagonalizing the mean-field Hamiltonian with the renormalized order parameters $m_{\uparrow}=m_{\rm C}-m_{\rm AF}$ and $m_{\downarrow}=m_{\rm C}+m_{\rm AF}$, where

[TABLE]

Here $\langle\cdot\rangle=\langle\Psi_{\rm G}|\cdot|\Psi_{\rm G}\rangle$ is now calculated based on $|\Psi_{\rm G}\rangle$. We note that to solve the Gutzwiller ground state of the initial Hamiltonian, one needs to determine the variational parameter $\alpha$ and the wave function $|\Psi_{\rm G}\rangle$ with renormalized $m_{\uparrow,\downarrow}$ iteratively until the energy $\langle H\rangle_{\rm G}$ is minimized. For the purpose of the present study, we only need to formally write down the ground state in the Gutzwiller form, and show that the quench dynamics can extract the information of the ground state.

VI.0.2 B. Correlations and flow equations

Compared with Eq. (S14), we have nonzero correlations with respect to the Gutzwiller wave function ($\mu,\nu=\pm$, ${\bf k}_{1}\neq{\bf k}_{2}$)

[TABLE]

where the BMF corrections

[TABLE]

Here we denote $\bar{\uparrow}=\downarrow$ and $\bar{\downarrow}=\uparrow$. Based on the above results, one can easily check that the generator $\eta(l)$ [see Eq. (S17)] and the ansatzes for creation operators [Eq. (II)] remain unchanged. Moreover, the leading-order flow equations take the same forms as Eqs. (II) and (II), except that the parameters $F_{\bf p^{\prime}qq^{\prime}}^{s_{2}s_{3}s_{4},\mu\gamma\nu}$ and $G_{\bf pp^{\prime}q^{\prime}}^{s_{1}s_{2}s_{4},\nu\mu\gamma}$ should be, respectively, replaced by $\widetilde{F}_{\bf p^{\prime}qq^{\prime}}^{s_{2}s_{3}s_{4},\mu\gamma\nu}$ and $\widetilde{G}_{\bf pp^{\prime}q^{\prime}}^{s_{1}s_{2}s_{4},\nu\mu\gamma}$, where

[TABLE]

Similar to Eq. (II), we can further simplify the calculations by only taking into account the major contribution, i.e.

[TABLE]

with $\widetilde{\cal F}_{\bf p^{\prime}qq^{\prime}}^{\mu\gamma\nu}\simeq{\cal F}_{\bf p^{\prime}qq^{\prime}}^{\mu\gamma\nu}+\alpha{\cal C}_{\mu\mu}^{\uparrow}({\bf p^{\prime}})[n_{\gamma-}^{\downarrow}({\bf q})-n_{\nu-}^{\downarrow}({\bf q}^{\prime})]+\alpha{\cal C}_{\gamma\gamma}^{\downarrow}({\bf q})[n_{\mu-}^{\uparrow}({\bf p}^{\prime})-n_{\nu+}^{\downarrow}({\bf q}^{\prime})]-\alpha{\cal C}_{\nu\nu}^{\downarrow}({\bf q^{\prime}})[n_{\mu-}^{\uparrow}({\bf p}^{\prime})-n_{\gamma-}^{\downarrow}({\bf q})]+{\cal O}(\alpha^{2})$ and $\widetilde{\cal G}_{\bf pp^{\prime}q^{\prime}}^{\nu\mu\gamma}\simeq{\cal G}_{\bf pp^{\prime}q^{\prime}}^{\nu\mu\gamma}+\alpha{\cal C}_{\nu\nu}^{\uparrow}({\bf p})[n_{\mu-}^{\uparrow}({\bf p}^{\prime})-n_{\gamma+}^{\downarrow}({\bf q}^{\prime})]+\alpha{\cal C}_{\mu\mu}^{\uparrow}({\bf p}^{\prime})[n_{\nu-}^{\uparrow}({\bf p})-n_{\gamma-}^{\downarrow}({\bf q}^{\prime})]+\alpha{\cal C}_{\gamma\gamma}^{\downarrow}({\bf q}^{\prime})[n_{\nu-}^{\uparrow}({\bf p})-n_{\mu-}^{\uparrow}({\bf p}^{\prime})]+{\cal O}(\alpha^{2})$. With these results we conclude that the forms of flow equations remain the same, while the BMF corrections only modify the parameters that depend on the initial momentum distribution. Therefore, we find that the equation of motion keeps the same form as Eq. (S40), but the parameters are corrected by the interactions. Also, the initial spin length at each single particle momentum can be less than one, since in the initial Gutzwiller ground state the different momentum states are correlated.

VI.0.3 C. Universal scaling on the BIS

From the results above, one conclude that the quench dynamics is governed by the equation of motion of the same form as Eq. (S40) but with modified parameters. With this we can expect that the emergent topology of the quench dynamics is not affected by the corrections, as long as the initial ground state is topologically trivial. We then focus on the BMF correction to the universal scaling. The BIS is determined by

[TABLE]

where $T$ denotes the evolution time, the mean-field density

[TABLE]

with $\widetilde{E}_{0}^{\sigma}\equiv\sqrt{E_{0}^{2}+2m_{\sigma}h_{z}+m_{\sigma}^{2}}$, and $\widetilde{\lambda}^{\sigma}_{2}({\bf k})$ are defined as in Eqs. (III) and (III) with $n^{\sigma}_{s\pm}({\bf k})$ being replaced by $n^{\sigma}_{s\pm}({\bf k})\pm\alpha{\cal C}_{ss}^{\sigma}({\bf k})$. Since

[TABLE]

where we define $\zeta\overset{\rm def}{=}\sum_{\begin{subarray}{c}{\bf q,p}\\ s_{1}s_{2}s_{3}s_{4}\end{subarray}}\Re\left[\Lambda_{\bf ppqq}^{s_{1}s_{2}s_{3}s_{4}}f_{s_{1}-}^{\uparrow\,*}({\bf p})f_{s_{2}-}^{\uparrow}({\bf p})f_{s_{3}-}^{\downarrow\,*}({\bf q})f_{s_{4}-}^{\downarrow}({\bf q})\right]$, we then know that on the BIS,

[TABLE]

and half of the amplitude reads

[TABLE]

If we define $\widetilde{Z}^{\sigma}_{0}\equiv Z^{\sigma}_{0}/(1-2\alpha\zeta)$, we still have the universal scaling for quench dynamics on the BIS

[TABLE]

where $f(m_{\sigma})=m_{\sigma}T_{0}$ and $g(\widetilde{Z}^{\sigma}_{0})=\sqrt{1-\widetilde{Z}_{0}^{\sigma\,2}}/\pi$. Note that the unknown constant $\alpha\zeta$ can be determined by measuring $Z^{\sigma}_{0}$ and $T_{0}$ of the dynamics at any two points on BIS, which satisfy the same scaling in the above Eq. (S70).