Quantum dynamics of impenetrable SU(N) fermions in one-dimensional lattices

TL;DR

This paper investigates the quantum quench dynamics of strongly interacting SU(N) fermions in one-dimensional lattices, revealing emergent correlations and distributions during expansion, and establishing a framework for analyzing their equilibration.

Contribution

It introduces an exact approach to study SU(N) fermion dynamics, uncovering phenomena like Gaussian correlations and quasimomentum to rapidity distribution transformation, and discusses equilibration via a generalized Gibbs ensemble.

Findings

Gaussian one-body correlations emerge during domain wall melting.

Quasimomentum distribution transforms into rapidity distribution over time.

Observables after equilibration are described by a generalized Gibbs ensemble.

Abstract

We study quantum quench dynamics in the Fermi-Hubbard model, and its SU($N$) generalizations, in one-dimensional lattices in the limit of infinite onsite repulsion between all flavors. We consider families of initial states with generalized Neel order, namely, initial state in which there is a periodic $N$-spin pattern with consecutive fermions carrying distinct spin flavors. We introduce an exact approach to describe the quantum evolution of those systems, and study two unique transient phenomena that occur during expansion dynamics in finite lattices. The first one is the dynamical emergence of Gaussian one-body correlations during the melting of sharp (generalized) Neel domain walls. Those correlations resemble the ones in the ground state of the SU($N$) model constrained to the same spin configurations. This is explained using an emergent eigenstate solution to the quantum dynamics. The second phenomenon is the transformation of the quasimomentum distribution of the expanding strongly interacting SU($N$) gas into the rapidity distribution after long times. Finally, we study equilibration in SU($N$) gasses and show that observables after equilibration are described by a generalized Gibbs ensemble. Our approach can be used to benchmark analytical and numerical calculations of dynamics of strongly correlated SU($N$) fermions at large $U$.

Equations

$${\hat{H}}_N=-J\sum _{l=1}^{L-1}\sum _{\sigma =1}^N\left[{\hat{f}}_l^{(\sigma )†}{\hat{f}}_{l+1}^{(\sigma )}+{\hat{f}}_{l+1}^{(\sigma )†}{\hat{f}}_l^{(\sigma )}\right],$$

$$\underline{\sigma }=\{\{{\sigma }_j\};j=1,...,N_p;{\sigma }_j\ne {\sigma }_{j+1}\forall j<N_p\}.$$

$$\underline{\sigma }=\{\{{\sigma }_j\};j=1,...,N_p;{\sigma }_j=\left[(j-1)\text{mbox}modN\right]+1\},$$

$$C_l^{\sigma }(x;t)$$

$$C_l(x;t)=\sum _{\sigma =1}^N{C}_l^{\sigma }(x;t),$$

$$C_l^{\sigma }(x;t)=⟨{\Psi }_{\mathrm{SF}}(t)|{\hat{c}}_{l+x}^{†}{\hat{c}}_l{\hat{\mathcal{P}}}_{l,x}^{(\sigma )}|{\Psi }_{\mathrm{SF}}(t)⟩.$$

$${\hat{H}}_{\mathrm{SF}}=-\sum _{l=1}^{L-1}({\hat{c}}_l^{†}{\hat{c}}_{l+1}+{\hat{c}}_{l+1}^{†}{\hat{c}}_l),$$

$${\hat{\mathcal{P}}}_{l,x}^{(\sigma )}={\hat{\mathcal{M}}}_{l,x}{\hat{\mathcal{R}}}_l^{(\sigma )}.$$

$${\hat{\mathcal{M}}}_{l,x}=\prod _{j=l+1}^{l+x-1}\left(1-{\hat{c}}_j^{†}{\hat{c}}_j\right),$$

$${\hat{\mathcal{R}}}_l^{(\sigma )}(N)=\frac{1}{N}\sum _{k=0}^{N-1}{e}^{-\frac{2\pi i}{N}\sigma k}\exp \left[\frac{2\pi i}{N}k\sum _{j=1}^l{\hat{c}}_j^{†}{\hat{c}}_j\right].$$

$$C_l(x;t)=⟨{\Psi }_{\mathrm{SF}}(t)|{\hat{c}}_{l+x}^{†}{\hat{c}}_l{\hat{\mathcal{M}}}_{l,x}|{\Psi }_{\mathrm{SF}}(t)⟩.$$

$$|{\Psi }_{\mathrm{SF}}(t)⟩=e^{-i{\hat{H}}_{\mathrm{SF}}t}|{\Psi }_{\mathrm{SF}}^I⟩=\prod _{j=1}^{N_p}\sum _{m=1}^L{G}_{mj}(t){\hat{c}}_m^{†}|\varnothing ⟩,$$

$${\hat{\mathcal{P}}}_{l,x}^{(\sigma )}|{\Psi }_{\mathrm{SF}}(t)⟩=\frac{1}{N}\sum _{k=0}^{N-1}{e}^{-\frac{2\pi i}{N}\sigma k}\prod _{j=1}^{N_p}\sum _{m=1}^L{G}_{mj}^k(t){\hat{c}}_m^{†}|\varnothing ⟩,$$

$$G^{k}_{mj}(t)=\Bigg{\{}\begin{tabular}[]{ccl}$e^{\frac{2\pi i}{N}k}\,G_{mj}(t)$,&&$m\leq l$\\ $0$,&&$l<m<l+x$\\ $G_{mj}(t)$,&&${\rm otherwise}$\end{tabular}\,$$

$$C_l^{\sigma }(x;t)=\frac{1}{N}\sum _{k=0}^{N-1}{e}^{-\frac{2\pi i}{N}\sigma k}$$

$$-\det [{\mathbf{G}}^{\prime }(t;l{)}^{†}{\mathbf{G}}^{\prime k}(t;l+x)]).$$

$$C_l(x;t)=$$

$$-\det [{\mathbf{G}}^{\prime }(t;l{)}^{†}{\mathbf{G}}^{\prime k=0}(t;l+x)],$$

$$|{\Psi }^I⟩=\prod _{j=1}^{N_p}{\hat{f}}_j^{({\sigma }_j)†}|\varnothing ⟩,$$

$$n(l;t{)}_{\infty }=\frac{1}{\pi }\arccos \left(\frac{l-l_0}{2t}\right),$$

$${\hat{\mathcal{H}}}_{\mathrm{SF}}(t)=-t\sum _{l=1}^{L-1}(i{\hat{c}}_l^{†}{\hat{c}}_{l+1}+H.c.)+\sum _{l=1}^{L}l{\hat{n}}_l,$$

$$\mu (s)={\mu }_0-s,$$

$${\hat{\mathcal{H}}}_{\mathrm{SF}}^{\mathrm{LDA}}(s)=-\sum _l(i{\hat{c}}_l^{†}{\hat{c}}_{l+1}+H.c.)-\mu (s)\hat{N}.$$

$$n(s)=\frac{1}{\pi }\arccos \left(\frac{s}{2}\right),$$

$$|C_l(x;t)|\stackrel{?}{=}|C_s(x)|.$$

$$|C_s(x)|=n(s)e^{-x^2/x_0^2},$$

$$m(k;t)=\frac{1}{L}\sum _{l,x}{e}^{-ikx}{C}_l(x;t),$$

$${\hat{H}}_N^{\prime }={\hat{H}}_N+\frac{1}{R^2}\sum _{l=1}^L\sum _{\sigma =1}^N(l-l_0{)}^2{\hat{f}}_l^{†(\sigma )}{\hat{f}}_l^{(\sigma )},$$

$${\hat{\mathcal{H}}}_{\mathrm{SF}}^{\prime }(t)=$$

$$-\frac{t^2}{R^2}\sum _{l=1}^{L-2}({\hat{c}}_{l+2}^{†}{\hat{c}}_l+H.c.)+\frac{1}{R^2}\sum _{l=1}^L(l-l_0{)}^2{\hat{c}}_l^{†}{\hat{c}}_l,$$

$${\hat{\mathcal{H}}}_{\mathrm{SF}}^{\prime \mathrm{LDA}}(t,l)=$$

$$-\frac{t^2}{R^2}\sum _j({\hat{c}}_{j+2}^{†}{\hat{c}}_j+H.c.)+\mu (l)\hat{N}.$$

Full text

I Introduction

Nonequilibrium dynamics of isolated many-body quantum systems is currently one of the most active research fields in many-body quantum physics polkovnikov_sengupta_review_11 ; eisert_friesdorf_review_15 ; dalessio_kafri_16 . A special class of systems, called integrable systems, represent an important cornerstone in this field Calabrese_2016 ; essler_fagotti_2016 ; calabrese_cardy_2016 ; cazalilla_chung_2016 ; bernard_doyon_2016 ; caux_2016 ; vidmar16 ; ilievski_medenjak_2016 ; langen_gasenzer_2016 ; vasseur_moore_2016 ; deluca_mussardo_2016 . They possess an extensive number of (quasi-local) conserved quantities, and after equilibration expectation values of few-body observables can be described by generalized Gibbs ensembles (GGEs) rigol_dunjko_07 ; vidmar16 . Several experimental groups have been exploring many-body quantum dynamics with ultracold atoms close to integrability kinoshita_wenger_06 ; gring_kuhnert_12 ; ronzheimer13 ; meinert13 ; fukuhara13a ; fukuhara13b ; hild_fukuhara_14 ; vidmar15 ; langen_erne_15 ; scherg_18 ; tang_kao_18 .

In this work we study the paradigmatic one-dimensional (1D) Fermi-Hubbard model essler-book and its generalized 1D SU($N$) versions cazalilla_2009 ; Gorshkov_2010 ; Taie_2012 ; Pagano_2014 ; Scazza_2014 ; Zhang_2014 ; cazalilla_rey_14 ; nataf_mila_14 ; volosniev_14 ; dufour_nataf_15 ; decamp_armagnat_16 ; decamp_junemann_16 ; laird_shi_17 ; jen_yip_18 ; zinner_18 . The Fermi-Hubbard model is a minimal model used to describe correlated electrons in solids, and provides an accurate description of ongoing experiments with ultracold fermions in optical lattices esslinger_review_10 . After early theoretical studies of quantum quenches in the Fermi-Hubbard model in dimensions higher than one moeckel_kehrein_08 ; moeckel_kehrein_09 ; eckstein_kollar_09 ; eckstein_kollar_10 ; schiro_fabrizio_10 ; schiro_fabrizio_11 , recent works have focused on quantum quenches in the 1D Fermi-Hubbard model hm08 ; hm09 ; kajala11 ; langer12 ; bolech12 ; hamerla_uhrig_13 ; hamerla_uhrig_13_b ; iyer_mondaini_14 ; riegger_orso_15 ; mei16 ; yin_radzihovsky_16 ; schlunzen_joost_17 ; bertini_tartaglia_2017 ; bleicker_uhrig_18 , which is a quintessential integrable model essler-book .

We introduce an exact approach to study quantum dynamics of 1D SU($N$) fermions, of which the Fermi-Hubbard model is the $N=2$ case, in the limit of infinite onsite repulsion. A key requirement of our approach is that the states of which we study the quantum evolution must have consecutive impenetrable SU($N$) fermions with different spin flavors. In those configurations, since the spin order is preserved at all times, one can think of each fermion as a distinguishable quantum particle. This allows us to use spin-charge separation to treat the charge degrees of freedom as spinless fermions, and the spin degrees of freedom by means of nonlocal strings of operators acting on spinless fermion wave functions. Using this decomposition, we develop an efficient (polynomial-time) way to compute one-body correlations. The ideas behind this approach were recently used to describe equilibrium properties of impenetrable 1D SU($N$) fermions zhang18 , and can be viewed as a fermionic generalization of the Jordan-Wigner transformation, which maps spin-$\frac{1}{2}$ models onto spinless fermions jordan_wigner_28 ; lieb61 . We focus on dynamics of states that, in addition of having consecutive fermions with distinct spin flavors, have periodic $N$-spin patterns that we call generalized Neel order.

As a first application of our approach, we study the expansion (“melting”) of generalized Neel domain walls such as the one shown in Fig. 1 for the SU(2) case. This type of setup has been widely studied in the literature, in particular for lattice models without internal degrees of freedom antal99 ; karevski02 ; ogata02 ; hunyadi04 ; rigol04 ; rigol05a ; gobert05 ; antal08 ; eisler09 ; lancaster10 ; mossel10 ; santos11 ; eisler13 ; sabetta_misguich_13 ; alba14 ; hauschild15 ; vasseur15 ; lancaster16 ; bertini16 ; eisler_maislinger_16 ; hauschild_heidrichmeisner_16 ; vidmar_iyer_17 ; deluca_collura_17 ; ljubotina_znidaric_17 ; kormos_17 ; ribeiro17 ; herbrych_feiguin_17 ; stephan_17 ; eisler_bauernfeind_17 ; fagotti_17 ; lorenzo_17 ; bastianello_deluca_18 ; collura_deluca_18 ; bertini_piroli_18 ; bertini_18 ; denardis_bernard_18 ; alba_18 ; alba_19 ; mastyan_bertini_19 ; ljubotina_sotiriadis_19 ; mazza_perfetto_18 . We show that the resulting current-carrying state exhibits a dynamical emergence of one-body correlations with a Gaussian decay. These correlations resemble the ones in the ground state of the same model, constrained to the same spin configurations zhang18 . We explain this observation using an emergent eigenstate solution to quantum dynamics vidmar_iyer_17 , combined with the local density approximation.

As a second application, we study the expansion of harmonically trapped SU($N$) fermions after suddenly turning off the trap. Those quenches provide a fertile playground to unveil remarkable properties of correlated many-body systems sutherland_98 ; rigol04 ; rigol05 ; rodriguez_manmana_06 ; kajala11 ; bolech12 ; vidmar13 ; boschi14 ; vidmar15 ; ren_wu_15 ; mei16 ; vidmar_iyer_17 ; xu_rigol_17 ; vidmar_xu_17 ; siegl_mistakidis_18 ; zhang_meurice_19 . They can be viewed as nontrivial time-of-flight expansions that occur in the presence of interactions. Here we demonstrate that, during the expansion, the quasimomentum distribution of impenetrable SU($N$) fermions evolves into the rapidity distribution, which is the initial quasimomentum distribution of the underlying spinless fermions to which the impenetrable SU($N$) fermions are mapped. Transformations of quasimomentum distributions of expanding fermions into rapidity distributions have been previously discussed in the context of the Fermi-Hubbard model bolech12 ; mei16 .

As a final application of our approach, we study the equilibration dynamics of SU($N$) fermions in a box trap. In that case, the quasimomentum distribution function after equilibration is described by means of a GGE in which the conserved quantities are the occupations of the single-particle eigenstates of the underlying spinless fermions to which the impenetrable SU($N$) fermions are mapped. This shows that the GGE can be used to describe observables after equilibration in integrable multiflavor fermionic systems with infinite onsite repulsion.

In addition to studying phenomena of relevance to current experiments with ultracold quantum gases, the approach introduced here can also be used to benchmark analytical and numerical calculations of quantum dynamics of strongly correlated SU($N$) fermions at large $U$, and specifically of the integrable SU(2) Fermi-Hubbard model. Quantum dynamics of integrable models non-mappable onto noninteracting ones have proved challenging. Analytic progress has been made in the context of the spin-$\frac{1}{2}$ XXZ chain, for which exact steady-state properties have been obtained for some classes of initial states either using special properties of those states (e.g, within the quench action approach caux_essler_13 ; wouters_denardis_14 ; brockmann_wouters_14 ; pozsgay_mestyan_14 ; mestyan_pozsgay_15 ; alba_calabrese_16 ; caux_2016 and the quantum transfer matrix formalism piroli_pozsgay_2017_a ; piroli_pozsgay_2017_b ), or using the (quasi-local) conserved quantities of the model fagotti_14 ; ilievski15 ; ilievski_quinn_2016 ; piroli_vernier_16 ; ilievski_quinn_17 ; piroli_vernier_17 ; pozsgay_vernier_17 . More recently, interest has grown in understanding systems with internal degrees of freedom, so-called nested systems, such as the Fermi-Hubbard model. In those systems, exact steady-state properties can also be obtained for a few classes of initial states. Examples include the two-component Lieb-Liniger gas robinson_caux_16 ; robinson_caux_16_b , the Fermi-Hubbard model bolech12 ; mei16 ; bertini_tartaglia_2017 ; enej_17 , and the Lai-Sutherland model mestyan_bertini_2017 ; piroli_vernier_19_a ; piroli_vernier_19_b . Nevertheless, their dynamics remains widely unexplored, and a systematic study of the GGE in such systems remains challenging due to difficulties in obtaining the complete set of conserved quantities.

The rest of the paper is organized as follows. In Sec. II, we introduce the nonequilibrium setup and the exact approach used to describe quantum dynamics after quenches. We then discuss two applications to transient phenomena: we study the melting of generalized Neel domain walls in Sec. III, and the sudden expansion of initially trapped systems in Sec. IV. In Sec. V, we discuss equilibration to the GGE. A summary of our results is presented in Sec. VI.

II Setup and quantum dynamics

Here, we discuss general considerations about the systems studied and introduce the exact approach developed to describe their dynamics.

II.1 Distinguishable quantum particles

We consider a generalized 1D Fermi-Hubbard model for SU($N$) fermions, $N$ is the number of flavors, with infinite onsite repulsion. For open chains with $L$ sites (on which we focus here), the Hamiltonian can be written as

[TABLE]

where $\sigma$ is the spin flavor ($\sigma\in\{1,\dots,N\}$ in our notation). Infinite onsite repulsion is enforced by the constraints $\hat{f}^{(\sigma)\dagger}_{l}\hat{f}^{\dagger(\sigma^{\prime})}_{l}=\hat{f}^{(\sigma)}_{l}\hat{f}^{(\sigma^{\prime})}_{l}=0$, where $\hat{f}^{(\sigma)\dagger}_{l}$ ($\hat{f}^{(\sigma)}_{l}$) is the creation (annihilation) operator of a fermion with flavor $\sigma$ at site $l$. The traditional Fermi-Hubbard Hamiltonian corresponds to $N=2$. We set the hopping amplitude $J$ and the lattice spacing to unity.

An important property of Hamiltonian $\hat{H}_{N}$, which is a consequence of infinite onsite repulsion, is that it preserves the order of spin flavor configurations $\underline{\sigma}=\{\sigma_{1},\dots,\sigma_{N_{p}}\}$, where $\sigma_{j}\in\{1,\dots,N\}$, and $N_{p}$ is the number of particles. In a sector with $\underline{\sigma}$, any wave function in the occupation number basis is a linear combination of base kets $|\varphi_{\underline{x},\underline{\sigma}}\rangle=\prod_{j=1}^{N_{p}}\hat{f}_{x_{j}}^{(\sigma_{j})\dagger}|\emptyset\rangle$, where every base ket denotes a different set of occupied sites $\underline{x}=\{x_{1},\dots,x_{N_{p}}\}$, with $x_{j}\in\{1,\dots,L\}$ and $x_{1}<x_{2}<\cdots<x_{N_{p}}$.

Here we are interested in quantum dynamics of initial states with a given spin configuration $\underline{\sigma}$. The most general states for which our approach is applicable are states with $\underline{\sigma}$ such that every pair of consecutive fermions carries distinct spin flavors

[TABLE]

We call the impenetrable SU($N$) model governed by $\hat{H}_{N}$ (1), constrained to a single spin configuration $\underline{\sigma}$ that obeys Eq. (2), a model of distinguishable quantum particles (DQPs).

In the SU(2) case, the DQP model is the traditional Fermi-Hubbard model within the sectors with states that have alternating spin flavors. An example of such an spin ordering is shown in Fig. 1(a). In the SU($N$) case with $N>2$, different realizations of spin configurations satisfying Eq. (2) are possible. For simplicity, we focus on periodic $N$-spin patterns with consecutive fermions carrying distinct spin flavors, as described by

[TABLE]

which we call generalized Neel order for SU($N$) particles. Sectors with desired spin configurations can be accessed in experiments with ultracold atoms in optical lattices via spin-resolved manipulation techniques weitenberg_11 ; hild_fukuhara_14 .

II.2 Spin-charge decomposition

Our main goal is to calculate nonequilibrium properties of spin-resolved one-body correlations $C^{\sigma}_{l}(x;t)$ between sites $l$ and $l+x$

[TABLE]

where $|\Psi(t)\rangle$ is the time-evolving state of the DQP model. We are also interested in computing the total one-body correlations

[TABLE]

which contain the contributions from all spin flavors.

Central to our approach is exploiting the separation of spin and charge degrees of freedom at infinite onsite repulsion. In particular, we use a compact representation of charge degrees of freedom as spinless fermions and spin degrees of freedom as nonlocal string of operators (projectors) expressed in terms of spinless fermions. Formal steps demonstrating the existence of such representations were reported in Ref. kumar_09 . Here we use an explicit representation of the spin and charge decomposition for states that exhibit a generalized Neel order [see Eq. (3)]. This representation was introduced in Ref. zhang18 and allows one to efficiently (in polynomial time) calculate spin-resolved one-body correlations numerically.

Our approach is not based on the Bethe ansatz solution of the infinitely repulsive Fermi-Hubbard model lieb_wu_68 ; ogata_shiba_90 ; parola_sorella_90 ; parola_sorella_92 ; penc_96 ; penc_hallberg_97 ; izergin_pronko_98 ; essler-book . It has analogies with the Jordan-Wigner transformation used to map spin-$\frac{1}{2}$ systems onto spinless fermions. Since the pioneering work by Lieb, Schultz, and Mattis lieb61 , the Jordan-Wigner transformation has been used as a standard tool to study quantum magnetism. It can be efficiently implemented numerically using properties of Slater determinants rigol05a , and has been used to study quantum quenches in 1D hard-core boson systems rigol04 ; rigol05 ; rigol05a ; rigol_dunjko_07 ; rigol_muramatsu_06 ; vidmar15 ; vidmar_iyer_17 ; xu_rigol_17 ; vidmar_xu_17 .

Using a compact representation of the spin-charge decomposition of states with a generalized Neel order zhang18 , we rewrite Eq. (4) as

[TABLE]

Here, $|\Psi_{\rm SF}(t)\rangle$ is the wave function that describes the charges (spinless fermions) at time $t$, while the time-independent projection operator ${\cal\hat{P}}^{(\sigma)}_{l,x}$ properly accounts for the spins. The Hamiltonian that governs the dynamics of the charge degrees of freedom $|\Psi_{\rm SF}(t)\rangle$ is the spinless fermion Hamiltonian

[TABLE]

where $\hat{c}^{\dagger}_{l}$ ($\hat{c}_{l}$) is the spinless fermion creation (annihilation) operator at site $l$.

The spin projection operator ${\cal\hat{P}}^{(\sigma)}_{l,x}$ is in general a nonlocal multi-body operator zhang18 . It is constructed as the product of two operators

[TABLE]

The role of the operator ${\cal\hat{M}}_{l,x}$ is to prevent exchange of fermions. It is defined as

[TABLE]

and guarantees that all lattice sites from ${l+1}$ through $l+x-1$ are empty. The role of ${\cal\hat{R}}^{(\sigma)}_{l}$ is to target spin flavor $\sigma$ at site $l$. For spin configurations that have generalized Neel order [Eq. (3)], ${\cal\hat{R}}^{(\sigma)}_{l}(N)$ is defined as

[TABLE]

This operator counts the number of fermions from sites $1$ through $l$, and ensures that it is consistent with a fermion with spin flavor $\sigma$ occupying site $l$.

For the total one-body correlations $C_{l}(x;t)$, the spin projection operators can be further simplified. In analogy to equilibrium calculations zhang18 , we get

[TABLE]

$C_{l}(x;t)$ is hence independent of the number of spin flavors $N$. Equation (11) is valid for any spin orderings $\underline{\sigma}$ that fulfill Eq. (2).

II.3 Implementation for quantum quenches

To evaluate Eq. (6), we follow the approach used for hard-core bosons in Refs. rigol04 ; rigol05a . We first note that, given an initial state whose charge degrees of freedom can be written as a Slater determinant $|\Psi_{\rm SF}^{I}\rangle=\prod_{j=1}^{N_{p}}\sum_{m=1}^{L}A_{mj}\hat{c}_{m}^{\dagger}|\emptyset\rangle$, where $A_{mj}$ are matrix elements of an $L\times N_{p}$ matrix ${\bf A}$, the charge degrees of freedom of the time-evolving state can be written as

[TABLE]

where $G_{mj}(t)$ are the matrix elements of an $L\times N_{p}$ matrix ${\bf G}(t)={\bf U}e^{-i{\bf E}_{\rm SF}t}{\bf U}^{\dagger}{\bf A}$. Here, ${\bf E}_{\rm SF}$ is a diagonal matrix that satisfies ${\bf H}_{\rm SF}{\bf U}={\bf U}{\bf E}_{\rm SF}$, and ${\bf H}_{\rm SF}$ is the single-particle matrix representation of $\hat{H}_{\rm SF}$.

A crucial next step is to evaluate ${\cal\hat{P}}^{(\sigma)}_{l,x}|\Psi_{\rm SF}(t)\rangle$. As in Ref. zhang18 , the spin projection operator defined in Eqs. (8)–(10) acting on $|\Psi_{\rm SF}(t)\rangle$ yields a linear combination of Slater determinants,

[TABLE]

where

[TABLE]

are matrix elements of a $L\times N_{p}$ matrix ${\bf G}^{k}(t)$.

The next to last step to evaluate Eq. (6) is to act with the spinless fermion annihilation and creation operators on ${\cal\hat{P}}^{(\sigma)}_{l,x}|\Psi_{\rm SF}(t)\rangle$. We first rewrite Eq. (6) as $C^{\sigma}_{l}(x;t)=\delta_{x,0}\langle\Psi_{\rm SF}(t)|{\cal\hat{P}}^{(\sigma)}_{l,x}|\Psi_{\rm SF}(t)\rangle-\langle\Psi_{\rm SF}(t)|\hat{c}_{l}\hat{c}^{\dagger}_{l+x}{\cal\hat{P}}^{(\sigma)}_{l,x}|\Psi_{\rm SF}(t)\rangle$. In the second term, $\hat{c}^{\dagger}_{j}|\Psi_{\rm SF}(t)\rangle$ and $\hat{c}^{\dagger}_{j}{\cal\hat{P}}^{(\sigma)}_{l,x}|\Psi_{\rm SF}(t)\rangle$ result in adding an extra column to the matrices representing the Slater determinants involved, yielding ${\bf G}(t)\to{\bf G}^{\prime}(t;j)$ and ${\bf G}^{k}(t)\to{\bf G}^{\prime k}(t;j)$, respectively, where $G_{l,N_{p}+1}^{\prime}=G_{l,N_{p}+1}^{\prime k}=\delta_{l,j}$ ($G_{li}^{\prime}=G_{li}$ and $G_{li}^{\prime k}=G_{li}^{k}$, for $i\leq N_{p}$). Finally, computing the inner product of the resulting Slater determinants, one obtains

[TABLE]

Note that, for time-evolving states after quantum quenches, Eq. (15) is equivalent to Eq. (19) in Ref. zhang18 .

We are particularly interested in the total one-body correlation function $C_{l}(x;t)$ defined in Eq. (5). It was shown in Ref. zhang18 that $C_{l}(x)$ in equilibrium exhibits weak finite-size effects and allows one to learn about universal properties of DQPs that are independent of the specific spin configuration $\underline{\sigma}$. $C_{l}(x;t)$ can be computed as

[TABLE]

where the matrix elements of ${\bf G}^{k=0}(t)$ are given by $G_{mj}^{k}(t)$ in Eq. (14) for $k=0$.

Before moving on to applications, we comment on a subtlety regarding correlations that are summed over all spin flavors. The total site occupations $n(l;t)=C_{l}(0;t)$ are identical to those of the noninteracting spinless fermions onto which we map the impenetrable SU($N$) fermions. On the other hand, for $x>1$, the total one-body correlations $C_{l}(x;t)$ are in general different from those of the spinless fermions. In fact, in Ref. zhang18 it was shown that $C_{l}(x)$ exhibits a Gaussian decay in the ground state of the DQP model, while it exhibits a power-law decay in the ground state of spinless fermions.

III Dynamical emergence of Gaussian one-body correlations

As a first application of the approach introduced in the previous section, we study the expansion of (generalized) Neel domain walls. The initial state $|\Psi^{I}\rangle$ for the Fermi-Hubbard model ($N=2$) is shown in Fig. 1(a). It is a product state with particles occupying the $N_{p}$ leftmost lattice sites,

[TABLE]

and the spin flavors exhibit the generalized Néel order defined in Eq. (3). A computational implementation of Eqs. (15) and (II.3) allows us to straightforwardly solve systems with of the order of $L\sim 10^{3}$ sites. We set $L=2N_{p}$ in our numerical calculations and only consider times at which particles moving to the right (or holes moving to the left) have not yet reached the lattice boundaries, namely, times $t<N_{p}/v_{\rm max}$, where $v_{\rm max}=2$ is the maximal group velocity in the lattice.

As $N_{p}\rightarrow\infty$, our setup is equivalent to the expansion of a semi-infinite domain wall in an infinite lattice. Such an expansion gives rise to a current-carrying steady state. In lattice Hamiltonians without internal degrees of freedom, properties of the current-carrying states emerging from an initial domain wall, as defined in Eq. (17), have been studied for spinless fermions antal99 , hard-core bosons rigol04 , and for the anisotropic Heisenberg model gobert05 .

Figures 2(a) and 2(b) show $n^{\sigma}(l;t=0)$ in the initial generalized Neel domain wall for SU(2) and SU(4) systems, respectively, with 24 fermions. Figures 2(c)–2(f) show the corresponding $n^{\sigma}(l;t)$ at time $t=8$. The latter make apparent an interesting feature of the domain-wall melting in our setup. At short times, the spin-resolved site occupations $n^{\sigma}(l;t)\equiv\langle\Psi(t)|\hat{f}^{(\sigma)\dagger}_{l}\hat{f}^{(\sigma)}_{l}|\Psi(t)\rangle$ resolve individual quantum particles quite well. The “resolution” increases with increasing $N$, as made apparent by comparing $N=2$ and $4$. Figures 2(c)–2(f) also show the total site occupations $n(l;t)=\sum^{N}_{\sigma=1}n^{\sigma}(l;t)=C_{l}(0;t)$ (dashed lines), which are identical to those of the corresponding spinless fermions.

As shown in an exact calculation antal99 , as well as using hydrodynamic arguments antal08 , $n(l;t)$ during the domain-wall melting obeys the scaling form

[TABLE]

when both $l,t\to\infty$, while their ratio is kept fixed within the region $-1\leq(l-l_{0})/(2t)\leq 1$. Here, $l_{0}$ denotes the initial domain-wall edge. Unless stated otherwise, in this section we set $l_{0}=N_{p}$. Outside the scaling region, site occupations simply equal $n(l;t)=1$ (to the left) and $n(l;t)=0$ (to the right). The thin solid lines in Figs. 2(c)–2(f) show $n(l;t)_{\infty}$. They fit very well the numerical results in finite systems. We identify the current-carrying steady state as that in the scaling regime, namely, as that in which local properties only depend on the ratio $(l-l_{0})/t$.

Figures 2(g) and 2(h) show the integrated number of particles of a given spin flavor $\sigma$ that have crossed $l_{0}$ ($l_{0}=N_{p}+1/2$ in Fig. 2) at time $t$, defined as $Q^{\sigma}(t)=\sum_{l>l_{0}}n^{\sigma}(l;t)$. $Q^{\sigma}(t)$ exhibits clear plateaus at integer values (better seen for $N=4$). Those values increase by one each time that a particle from the domain wall (with a given $\sigma$) crosses the initial edge. In contrast, a steady linear growth is seen for spinless fermions [defined as $\sum_{\sigma}Q^{\sigma}(t)$, shown as dashed lines in Figs. 2(g) and 2(h)] reflecting the ballistic nature of the charge current in the system.

III.1 Emergent eigenstate solution

One of the most remarkable features of the current-carrying states generated by the melting of domain walls such as the ones in Fig. 2, observed in studies involving spinless fermions antal99 ; eisler09 , hard-core bosons rigol04 ; rigol05a , Bose-Hubbard rodriguez_manmana_06 and Fermi-Hubbard hm08 models, and Bose gases denardis_panfil_18 , is the dynamical emergence of one-body correlations with ground-state-like behavior. We explore one-body correlations in our time-evolving states to see whether they exhibit a Gaussian decay like the one observed in the ground state of impenetrable SU($N$) fermions constrained to the same spin configurations zhang18 .

A theoretical framework that provides an understanding of the dynamical emergence of ground-state-like correlations during expansion dynamics was put forward in Ref. vidmar_iyer_17 . It was dubbed emergent eigenstate solution to quantum dynamics as it allows one to construct a local operator (called the emergent Hamiltonian) of which the time-evolving state is an eigenstate. For the initial domain wall of spinless fermions, the emergent local Hamiltonian is (up to a constant)

[TABLE]

where $\hat{n}_{l}=\hat{c}^{\dagger}_{l}\hat{c}_{l}$ is the particle-number operator at site $l$ vidmar_iyer_17 . Note that the expansion time $t$ is a parameter in $\hat{\cal{H}}_{\rm{SF}}(t)$. The time-evolving state is the ground state of $\hat{\cal{H}}_{\rm{SF}}(t)$ (up to corrections that vanish exponentially with system size) as long as particles (holes) do not reach the lattice boundary, i.e., as long as $n(1;t)=1$ and $n(L;t)=0$ vidmar_iyer_17 . Hence, in our setup, the ground state of $\hat{\cal{H}}_{\rm{SF}}(t)$ in Eq. (19), denoted by $|\Phi_{\rm{SF}}(t)\rangle$, provides an accurate (exact in the limit $N_{p}\rightarrow\infty$) description of the time-evolving state for the charge degrees of freedom of the impenetrable SU($N$) fermions $|\Psi_{\rm SF}(t)\rangle$. We then calculate one-body correlations $C^{\sigma}_{l}(x;t)$ via $|\Phi_{\rm{SF}}(t)\rangle$ using the spin projectors introduced in Sec. II.2.

While in previous applications of the emergent eigenstate solution $|\Phi_{\rm{SF}}(t)\rangle$ was evaluated exactly numerically vidmar_iyer_17 ; vidmar_xu_17 ; modak_vidmar_17 , here we use the local density approximation (LDA) in the scaling regime to simplify $\hat{\cal{H}}_{\rm{SF}}(t)$. This approach is similar to the generalized hydrodynamic approach bertini16 ; castroalvaredo_doyon_16 , in which a GGE is constructed for each rescaled position of the current-carrying state. We introduce, for $t>0$, an effective local chemical potential at the (continuous) rescaled position $s=(l-l_{0})/t$,

[TABLE]

where $\mu_{0}$ is the global chemical potential. For each $s$, the system is homogeneous with an effective chemical potential $\mu(s)$:

[TABLE]

We test the emergent eigenstate solution plus LDA by analytically calculating the site occupations $n(l;t)$ for spinless fermions. They follow $n(s)=\arccos[-\mu(s)/2]/\pi$ when $-2<\mu(s)<2$, and $n(s)=0$ [$n(s)=1$] for $\mu(s)\leq-2$ [$\mu(s)\geq 2$]. We fix $\mu_{0}$ using the total particle number, $N_{p}=\int n(s)ds$, which yields the final expression for site occupations

[TABLE]

when $-1\leq s/2\leq 1$. This is the result in Eq. (18).

III.2 Total one-body correlations

Next, we study the off-diagonal one-body correlations in the current-carrying state. They are exactly zero in the initial product state [Eq. (17)], and become nonzero as a result of the expansion dynamics. To minimize finite-size effects (see Appendix A), we focus on the total correlations $C_{l}(x;t)$ [see Eq. (II.3)]. As an example, Fig. 3 shows the absolute value of all matrix elements of the correlation matrix $C_{l,j}=|C_{l}(j-l;t)|$ at times $t=20$ and $50$, in a system with $N_{p}=200$. Those results reveal a very fast decay of the one-body correlations with $x=j-l$, in stark contrast with the power-law decay observed for spinless fermions (see, e.g., Fig. 2 in Ref. vidmar_iyer_17 ).

Figure 4 shows one-body correlations vs $x$ measured with respect to sites with a given particle occupation $n(l;t)=n_{0}$, for a system larger than the one in Fig. 3. We answer three questions about the current-carrying steady state:

(i) Does $|C_{l}(x;t)|$ exhibit a Gaussian decay with $x$ for different values of $l$ and $t$?

(ii) Do off-diagonal one-body correlations obey the scaling solution

[TABLE]

(iii) Do the exact numerical results for $|C_{l}(x;t)|$ match the predictions from the emergent eigenstate solution plus LDA?

Figure 4(a) shows one-body correlations with respect to a site with $n_{0}=0.5$ at different times [see legends in Fig. 4(b)]. It is apparent that one-body correlations exhibit a Gaussian decay that is independent of time. They only depend on the filling of the reference site, equivalently, $|C_{l}(x;t)|=|C_{s}(x)|$. The same conclusion applies, at long times, to the results for $n_{0}=0.2$ shown in Fig. 4(b). They converge to a Gaussian decay that is independent of time, i.e., $|C_{l}(x;t)|=|C_{s}(x)|$ at long times. The onset of the scaling behavior for the off-diagonal one-body correlations occurs at later times as the reference site occupation departs from half-filling.

In Figs. 4(a) and 4(b), we show (as dashed lines) numerical results for one-body correlations in the ground state of the emergent local Hamiltonian [Eq. (21)] at the target particle filling. Those correlations are in excellent agreement with the ones in the current-carrying state at long times. We should add that the absolute value of the one-body correlations in the ground state of the emergent local Hamiltonian are identical to the ground-state correlations in the physical Hamiltonian that governs the dynamics, Eq. (1) vidmar_iyer_17 . The latter is known to exhibit one-body correlations with a Gaussian decay

[TABLE]

where $x_{0}\propto n(s)^{-1}$ zhang18 .

III.3 Quasimomentum distribution function

To conclude the study of the current-carrying state, we calculate the total quasimomentum distribution function

[TABLE]

which is routinely measured in experiments with ultracold quantum gases bloch08 . (Spin-resolved results are shown in Appendix A.)

Figure 5(a) shows $m(k;t)$ at different times for a system with $N_{p}=400$ particles. In the initial product state, in which fermions are localized, all quasimomenta are equally occupied, i.e., $m(k;0)=\bar{m}=N_{p}/L=0.5$. Two main features are apparent in the dynamics: a peak emerges at $k=\pi/2$, and the height of the peak increases with time. Figure 5(a) also makes apparent that the results obtained for $m(k;t)$ during the exact time evolution are indistinguishable from those obtained in the ground state of the emergent local Hamiltonian $\hat{\cal H}_{\rm SF}(t)$ [Eq. (19)].

The emergence of the peak at $k=\pi/2$ can be understood from the phase $i=\exp(i\pi/2)$ in the hopping term of $\hat{\cal H}_{\rm SF}(t)$. As a result, the minimum of the energy dispersion of the kinetic part of $\hat{\cal H}_{\rm SF}(t)$ is located at $k=\pi/2$. For 1D hard-core bosons, which can also be mapped onto spinless fermions, the dynamical emergence of the peak at $k=\pi/2$ rigol04 has been experimentally observed in experiments with bosons in 1D optical lattices vidmar15 .

For the current-carrying state of spinless fermions, it was shown analytically that the evolution of the quasimomenta occupations relative to the initial occupations $m(k;t)-\bar{m}$ is linearly proportional to time $t$ vidmar_iyer_17 . In Fig. 5(b), we plot the rescaled quasimomentum distribution $[m(k;t)-N_{p}/L]/t$ for the impenetrable SU($N$) fermions. We observe a perfect data collapse, confirming the expectation from Fig. 5(a) that the height of the peak increases linearly with time.

IV Unveiling the rapidity distribution

As a second application of the approach introduced in Sec. II, we study the expansion of harmonically trapped impenetrable SU($N$) fermions in a lattice after suddenly turning off the trap. We take the initial state to be the ground state of the infinitely repulsive (generalized) Fermi-Hubbard model [Eq. (1)] in the presence of a harmonic confining potential

[TABLE]

where $l_{0}=L/2$ is the trap center, and $R$ is the characteristic length of the trap. Since the observables studied in this section are summed over all spin flavors, our results apply equally to the two-component Fermi-Hubbard model ($N=2$) as they do to arbitrary $N$-flavor SU($N$) models, for initial states in which every pair of consecutive fermions carries distinct spin flavors [Eq. (2)]. The corresponding spinless fermion Hamiltonian, to which we map the impenetrable SU($N$) fermions, has the form $\hat{H}^{\prime}_{\rm SF}=\hat{H}_{\rm SF}+R^{-2}\sum_{l}(l-l_{0})^{2}\hat{c}_{l}^{\dagger}\hat{c}_{l}$.

At time $t=0$, the harmonic confinement is turned off (we set $R^{-1}=0$) and the system evolves under $\hat{H}_{N}$. We are interested in dynamics that occurs before particles reach the boundaries of the lattice [$n(1;t)=n(L;t)=0$], which are equivalent to dynamics in an infinite lattice.

IV.1 Dynamical emergence of rapidities

Figures 6(a) and 6(b) show density plots of the expanding total site occupation $n(l;t)$ and the quasimomentum distribution function $m(k;t)$, respectively, as functions of time. Figure 6(b) shows that, after some transient time and despite the fact that $n(l;t)$ expands at all times, $m(k;t)$ becomes independent of time. [In this section, we normalize $m(k;t)$ using the characteristic length of the initial trap, $m(k;t)=\sum_{l,x}e^{-ikx}C_{l}(x;t)/R$.]

Our main focus next is to understand the stationary quasimomentum distribution function observed at long times. In various 1D bosonic systems including lattice hard-core bosons rigol05 ; vidmar_xu_17 , as well as Tonks-Girardeau minguzzi05 and Lieb-Liniger gases buljan_pezer_08 ; jukic_klajn_09 ; iyer12 ; campbell15 ; piroli_calabrese_17 , it has been shown that during expansion dynamics the quasimomentum distribution function of the bosons undergoes a dynamical fermionization. For lattice hard-core bosons rigol05 ; vidmar_xu_17 and Tonks-Girardeau gases minguzzi05 , this fermionization is nothing but a transformation of the quasimomentum distribution function of the physical particles (impenetrable bosons) into the quasimomentum (rapidity) distribution of the underlying noninteracting fermions into which the former can be mapped. This transformation of quasimomentum distribution functions of physical particles into rapidity distributions under expansion dynamics is expected to be a generic phenomenon in integrable models sutherland_98 ; bolech12 ; campbell15 ; mei16 .

Figure 7(a) shows $m(k;t)$ in the initial state and at time $t=375$, a time at which $m(k;t)$ has already become stationary (corresponding spin-resolved results are shown in Appendix A). We consider $N_{p}=100$ particles and an initial trap with $R=200$. Remarkably, the stationary result obtained for $m(k;t)$ is identical to the quasimomentum distribution function $m_{\rm SF}(k)$ of the spinless fermions to which we mapped the impenetrable SU($N$) fermions. $m_{\rm SF}(k)$ does not evolve in time because the quasimomenta occupations of the spinless fermions are conserved under expansion dynamics. $m_{\rm SF}(k)$ is the distribution of rapidities in our model, and our results show that it can be unveiled by allowing the impenetrable SU($N$) fermions freely expand up to about two to three times their initial size.

Figure 7(b), on the other hand, shows that despite the agreement between the quasimomentum distribution of impenetrable SU($N$) fermions and the quasimomentum distribution function of spinless fermions, at long times, the off-diagonal one-body correlations of both systems are not the same (this also occurs for lattice hard-core bosons rigol05 ; vidmar_xu_17 and Tonks-Girardeau gases minguzzi05 ). The absolute value of the one-body correlations exhibits a Gaussian decay for impenetrable SU($N$) fermions at all times and a power-law decay for spinless fermions.

IV.2 Emergent eigenstate solution

Figure 7(b) also shows the one-body correlations of the impenetrable SU($N$) fermions at $t=0$. They exhibit a Gaussian decay as the ground state of homogeneous systems zhang18 . The fact that in Fig. 7(b) the decay of correlations has the same behavior in the ground state and during the expansion indicates that the time-evolving states share properties with the initial ground state.

Like for the initial states considered in Sec. III, one can find an emergent eigenstate solution for the expansion of spinless fermions after turning off a harmonic trap vidmar_xu_17 . The emergent local Hamiltonian for the corresponding spinless fermion system vidmar_xu_17 ; modak_vidmar_17 is (up to a constant)

[TABLE]

where $A(l;t)=\sqrt{1+[(2t/R^{2})(l-l_{0}+1/2)]^{2}}$ is the nearest-neighbor hopping amplitude, and the phase is $\phi(l;t)=\arctan[2t(l-l_{0}+1/2)/R^{2}]$.

In Fig. 7(b) we show (as thin solid lines) the correlations $|C_{l}(x;t)|$ obtained in the ground state of $\hat{\cal H}_{\rm SF}^{\prime}(t)$, which perfectly overlap with the exact results from the nonequilibrium time evolution. Moreover, one can also combine the LDA with the emergent eigenstate solution, as in Sec. III, to describe the one-body correlations $|C_{l}(x;t)|$ in the nonequilibrium steady state in terms of those in a homogeneous system.

The appropriate homogeneous Hamiltonian can be obtained from $\mathcal{\hat{H}}_{\rm SF}^{\prime}(t)$ replacing the harmonic trap by an effective local chemical potential, and making the amplitude $A(l;t)$ and the phase $\phi(l;t)$ site independent:

[TABLE]

In Fig. 8 we compare $|C_{l}(x;t)|$ calculated during the exact time evolution (symbols) with the results in the ground state of $\hat{\cal{H}}_{\rm SF}^{\prime\rm LDA}(t)$ in Eq. (IV.2), when the particle occupation matches the target site occupation $n(l;t)$ (lines). We show result for correlations measured from the trap center ($l=l_{0}=L/2$), for which $A(l;t)\approx 1$ and $\phi(l;t)\approx 0$, for systems with $N_{p}=100$ impenetrable SU($N$) fermions and an initial $R=100$. The results shown in Fig. 8 from both approaches yield a perfect agreement, and follow the functional form in Eq. (24).

V Equilibration and Generalized thermalization

As a third, and final, application of the approach introduced in Sec. II, we study the equilibration of impenetrable SU($N$) fermions in a finite lattice with open boundary conditions (equivalent to a box trap with infinitely high walls) after suddenly turning off a harmonic trap in which the impenetrable SU($N$) fermions are initially confined. The setup is similar to the one in Sec. IV. The main difference is that the expanding particles now reach the boundaries of the lattice and equilibrate by moving back and forth in a finite system. The stationary quasimomentum distribution function in this case is not the rapidity distribution. Also, after particles reach the lattice boundaries, the emergent eigenstate solution is no longer valid vidmar_iyer_17 .

Figure 9(a) shows the dynamics of the total site occupations $n(l;t)$. At long times, as expected, $n(l;t)$ equilibrates and becomes (nearly) homogeneous, $n(l;t)\to N_{p}/L$. In the light of the findings in Sec. IV, the dynamics of the total quasimomentum distribution function $m(k;t)$, shown in Fig. 9(b), is more remarkable. After approaching the rapidity distribution during the expansion, when the propagating fronts reach the lattice boundaries for the first time, $m(k;t)$ begins to change and a revival occurs in the occupation of low quasimomentum modes. Ultimately, after some oscillations with decreasing amplitude, $m(k;t)$ equilibrates to a new distribution, which is different from that of the spinless fermions to which the impenetrable SU($N$) fermions are mapped.

In Fig. 10(a), we show the evolution of the total occupation of three quasimomentum modes ($k=0,\,\pi/5$, and $\pi/3$), as functions of time (normalized using the size of the final lattice $L$), for three values of $L$. We select $N_{p}$ and $R$ for the initial harmonic trap such that $N_{p}/L$ and $N_{p}/R$ are the same in all cases. Figure 10(a) shows that, as expected, the time scale for the oscillations in $m(k;t)$ is set by $L$ (notice the data collapse for different values of $L$). Figure 10(a) also makes apparent that, for $t/L\gtrsim 3$, the total occupations of the modes shown, and $m(k;t)$ in general, are nearly time independent, i.e., they have equilibrated.

In Fig. 10(b), we show the initial total quasimomentum distribution function, the long-time equilibrated $m(k;t)$ (corresponding spin-resolved results are shown in Appendix A), and the quasimomentum distribution of the spinless fermions to which we map the impenetrable SU($N$) fermions. They can all be seen to be different, with $m(k=0;t=0)$ being greater than the total occupation of the zero quasimomentum mode after equilibration, which in turn is greater than the occupation of the zero quasimomentum mode for spinless fermions.

In order to describe $m(k;t)$ after equilibration, we construct the generalized Gibbs ensemble (GGE) rigol_dunjko_07 ; vidmar16 , whose density matrix is defined as

[TABLE]

where $\{\hat{I}_{k}\}$ are the conserved quantities [the occupations of the single-particle eigenstates of the noninteracting spinless fermions to which we map the impenetrable SU($N$) fermions, see Eq. (7)], $Z_{\rm GGE}=\operatorname{Tr}[e^{-\sum_{k}\lambda_{k}\hat{I}_{k}}]$ is the partition function, and $\{\lambda_{k}\}$ are fixed by the initial values of the conserved quantities, $\operatorname{Tr}[\hat{\rho}_{\rm GGE}\hat{I}_{k}]=\langle\Psi^{I}_{\rm SF}|\hat{I}_{k}|\Psi^{I}_{\rm SF}\rangle\equiv\langle\hat{I}_{k}\rangle^{I}$. From these conditions, it follows that $\lambda_{k}=\ln[(1-\langle\hat{I}_{k}\rangle^{I})/\langle\hat{I}_{k}\rangle^{I}]$, and that $Z_{\rm GGE}=\prod_{k}(1-\langle\hat{I}_{k}\rangle^{I})^{-1}$ rigol_dunjko_07 ; vidmar16 .

Computing total one-body correlation functions of impenetrable SU($N$) fermions in the GGE is straightforward using the grand canonical ensemble approach introduced in Ref. zhang18 . The off-diagonal part of the total one-body correlations in the GGE is given by

[TABLE]

where $\bf I$ is the identity matrix, $\bf A$ is a matrix with elements $A_{ij}=\delta_{il}\delta_{j(l+x)}$, $\bf\Lambda$ is a diagonal matrix whose elements are the Lagrange multipliers $\lambda_{k}$, and $\mathbf{M}_{l,x}$ is the matrix form of the projection operator ${\cal\hat{M}}_{l,x}$ [see Eq. (9)], which is a diagonal matrix in which the elements between $l+1$ and $l+x-1$ are zero and the other ones are one. Moreover, $\bf U$ is the unitary matrix that diagonalizes $\hat{H}_{\rm SF}$ in Eq. (7), where $\mathbf{U}\mathbf{H}_{\rm SF}\mathbf{U}^{\dagger}=\mathbf{E}$ (here $\mathbf{E}$ is a diagonal matrix with the single-particle eigenenergies). The diagonal part of the total one-body correlations in the GGE is given by

[TABLE]

and is identical to that of the spinless fermions.

The results in Fig. 10(a) show that the total occupation of the quasimomentum modes of the impenetrable SU($N$) fermions equilibrate at the values predicted by the GGE. In Fig. 10(b), the long-time result for $m(k;t)$ is indistinguishable from the GGE one, and they are both clearly different from the prediction of the grand canonical ensemble (GE) zhang18 in which the temperature and the chemical potential are set by the energy and the number of particles in the system dalessio_kafri_16 . This makes apparent that these systems do not equilibrate at the traditional statistical mechanics prediction (do not thermalize).

In order to understand how the differences between the GGE and the exact dynamics predictions behave as one increases the system size, we compute the following time average of the differences between $m(k;t)$ after equilibration and the GGE prediction:

[TABLE]

Results for $\Delta m$ vs $L$ are reported in the inset in Fig. 10(b). They show that $\Delta m\propto L^{-1}$, consistent with findings in studies of 1D hard-core bosons gramsch_rigol_12 ; wright_rigol_14 . These results suggest that observables after long times and the GGE predictions become identical in the thermodynamic limit ($L\to\infty$, with $N_{p}/L$ and $N_{p}/R$ fixed).

Finally, we discuss the long-time behavior of the total one-body correlations for the largest system size considered in Fig. 10 ($L=1500$). For every $x$, we define the time average

[TABLE]

Results for $\overline{C}_{l}(x)$ are shown in Fig. 11 for $l=L/2$. They reveal that during the equilibration process the Gaussian decaying correlations present in the initial state evolve toward a slower, exponential-like, decay at large $x$. Those correlations are well described by the GGE, and are very different from the stretched exponential decay predicted by the grand canonical ensemble zhang18 .

VI Summary

We introduced an exact approach to study the time evolution of spin-resolved one-body observables after quantum quenches in the (generalized) SU($N$) Fermi-Hubbard model at infinite repulsion. This approach is tailored for initial states that exhibit generalized Neel order, namely, in which there is a periodic $N$-spin pattern with consecutive fermions carrying distinct spin flavors. Our approach is based on a compact representation of charge degrees of freedom as spinless fermions and spin degrees of freedom as nonlocal string of operators (projectors) expressed in terms of spinless fermions. It can be used to benchmark numerical and analytical studies of quantum quenches involving (generalized) Neel order in the (generalized) Fermi-Hubbard model in the limit of very strong repulsive interactions.

We studied three unique phenomena that occur during expansion dynamics of impenetrable SU($N$) fermions far from equilibrium. They are within reach in current optical-lattice experimental setups. The first one occurs during the melting of generalized Neel domain walls, in which we unveiled a dynamical emergence of Gaussian (ground-state-like zhang18 ) correlations in the resulting current-carrying steady state. We explained this phenomenon using an emergent eigenstate solution to quantum dynamics vidmar_iyer_17 . The second one occurs during the sudden expansion of harmonically trapped impenetrable SU($N$) fermions into an empty lattice after turning off the trap. We showed that the total quasimomentum distribution function of the impenetrable SU($N$) fermions evolves toward the quasimomentum distribution function of the spinless fermions onto which we map the former. Namely, as seen in other integrable systems rigol05 ; vidmar_xu_17 ; minguzzi05 ; buljan_pezer_08 ; jukic_klajn_09 ; iyer12 ; campbell15 ; piroli_calabrese_17 ; sutherland_98 ; bolech12 ; campbell15 ; mei16 , the quasimomentum distribution function of the real particles transforms into the rapidity distribution at long expansion times. Finally, we studied the equilibration dynamics of impenetrable SU($N$) fermions in a box trap after turning off a harmonic trap. We showed that observables (such as the quasimomentum distribution function) after equilibration can be described using a GGE.

VII Acknowledgments

We acknowledge discussions with B. Bertini and W. de Roeck. Y.Z. and M.R. acknowledge support from NSF Grant No. PHY-1707482. L.V. acknowledges support from the Slovenian Research Agency (ARRS), research core funding No. P1-0044.

Appendix A Spin-resolved quasimomentum distribution function

In the main text, in order to minimize finite-size effects, we focused on total one-body correlations [defined in Eq. (5)] and the corresponding total quasimomentum distribution functions. A detailed study of finite-size effects in spin-resolved one-body correlations and quasimomentum distribution functions of impenetrable SU($N$) fermions in equilibrium was presented in Ref. zhang18 . There we showed that those spin resolved observables approach the averaged ones (the total ones divided by the number of spin flavors) as one increases the number of particles while keeping the number of spin flavors fixed.

A similar study of finite-size effects in spin-resolved one-body correlations and quasimomentum distribution functions out of equilibrium is beyond the scope of this work. As for the equilibrium systems studied in Ref. zhang18 , we find that those observables out of equilibrium approach the averaged ones when increasing the number of particles (for a fixed number of spin flavors). Finite-size effects are small for the system sizes considered in the main text. In Fig. 12, we show the spin-resolved quasimomentum distribution functions $m^{\sigma}(k;t)$ for impenetrable SU(2) fermions ($\sigma=1,2$). Results for the melting of the generalized Neel domain wall (Sec. III), the sudden expansion after turning off the harmonic trap (Sec. IV), and equilibration in the box trap (Sec. V), are shown in Figs. 12(a), 12(b), and 12(c), respectively. The results for individual spin flavors are almost indistinguishable from the averaged ones, which are the total ones divided by 2.