Many-Body Localization from Dynamical Gauge Fields

Here we consider two spatially overlapped one-dimensional chains, each of which consists of a series of double wells, as sketched in Fig. 1. In the first chain, each double well is formed by the $i$-th site and $(i+1)$-th site for even number $i$, and each double well contains exactly one $a$-atom and one $f$-atom. This means for each double well with a pair of sites $i$ and $i+1$, we can introduce a spin-$1/2$ operator $\hat{\bm{\tau}}_{i}$ for the $f$-atom

	$\displaystyle\hat{\tau}_{i}^{x}=\frac{1}{2}\left(\hat{f}^{\dagger}_{i}\hat{f}_{i+1}+\hat{f}^{\dagger}_{i+1}\hat{f}_{i}\right),$		(5)
	$\displaystyle\hat{\tau}_{i}^{y}=\frac{i}{2}\left(\hat{f}^{\dagger}_{i+1}\hat{f}_{i}-\hat{f}^{\dagger}_{i}\hat{f}_{i+1}\right),$		(6)
	$\displaystyle\hat{\tau}_{i}^{z}=\frac{1}{2}\left(\hat{f}^{\dagger}_{i}\hat{f}_{i}-\hat{f}^{\dagger}_{i+1}\hat{f}_{i+1}\right).$		(7)

Similarly, we can introduce another spin-$1/2$ operator $\hat{\bm{\sigma}}_{i}$ for the $a$-atom

	$\displaystyle\hat{\sigma}_{i}^{x}=\frac{1}{2}\left(\hat{a}^{\dagger}_{i}\hat{a}_{i+1}+\hat{a}^{\dagger}_{i+1}\hat{a}_{i}\right),$		(8)
	$\displaystyle\hat{\sigma}_{i}^{y}=\frac{i}{2}\left(\hat{a}^{\dagger}_{i+1}\hat{a}_{i}-\hat{a}^{\dagger}_{i}\hat{a}_{i+1}\right),$		(9)
	$\displaystyle\hat{\sigma}_{i}^{z}=\frac{1}{2}\left(\hat{a}^{\dagger}_{i}\hat{a}_{i}-\hat{a}^{\dagger}_{i+1}\hat{a}_{i+1}\right).$		(10)

In this way, we have introduced $\hat{\bm{\tau}}_{i}$ and $\hat{\bm{\sigma}}_{i}$ for every even number $i$. The Hamiltonian of this chain is simply a sum of that of each individual double well, Eq. (4), and takes the following form

$$\hat{H}=\sum\limits_{i=\text{even}}-4g\hat{\tau}_{i}^{z}\hat{\sigma}_{i}^{x}-2h\hat{\tau}_{i}^{x},$$

(11)

in terms of previously defined spin-$1/2$ operators. This Hamiltonian has included the effect of lattice shaking, the hopping of both $a$- and $f$-atoms inside the double wells, and the interactions between $a$- and $f$-atoms Bloch_Z2_gauge ; Bloch_Z2_gauge_theory . It is exactly the same model as that has been realized in the recent experiment Bloch_Z2_gauge .

Next, we consider the second chain constituted with $b$- and $d$-atoms that play the same role as $a$- and $f$-atoms in the first chain respectively. In this chain, each double well is formed by the $i$-th and $(i+1)$-th sites with odd $i$. Similar to the first chain, these double wells are disconnected and each double well contains exactly one $b$-atom and one $d$-atom. We can, therefore, introduce two spin-$1/2$ operators $\hat{\bm{\tau}}_{i}$ and $\hat{\bm{\sigma}}_{i}$ for every odd number $i$. Their definitions are given by Eqs. (5)–(7) upon substituting $f$-atom operators with $d$-atom operators and by Eqs. (8)–(10) after replacing $a$-atom operators with $b$-atom operators. In this way, we have introduced $\hat{\bm{\tau}}_{i}$ and $\hat{\bm{\sigma}}_{i}$ for every odd number $i$. Moreover, the same shaking protocol is applied to this chain, which generates the same Hamiltonian as Eq. (11) except that the summation is over odd number of $i$. At this point, we have defined two sets of spin operators $\hat{\bm{\sigma}}_{i}$ and $\hat{\bm{\tau}}_{i}$ for all links in this one-dimensional system.

Now we turn on interactions between the two chains. We emphasize again that these two chains are arranged to be spatially overlapping with each other such that atoms at the same site can interact with each other. The interactions between $a$- and $f$-atoms and between $b$- and $d$-atoms have already been taken into account when writing down the effective Hamiltonian Eq. (11). Among the rest interactions, we consider the case where the dominate contribution comes from the interactions between $a$- and $b$-atoms, $U\sum_{i}\hat{a}^{\dagger}_{i}\hat{a}_{i}\hat{b}^{\dagger}_{i}\hat{b}_{i}$. Noting that for even number $i$, we have

	$\displaystyle\hat{a}^{\dagger}_{i}\hat{a}_{i}$	$\displaystyle=1/2+\hat{\sigma}_{i}^{z};$		(12)
	$\displaystyle\hat{b}^{\dagger}_{i}\hat{b}_{i}$	$\displaystyle=1/2-\hat{\sigma}_{i-1}^{z}.$		(13)

And for odd number $i$, we have

	$\displaystyle\hat{a}^{\dagger}_{i}\hat{a}_{i}$	$\displaystyle=1/2-\hat{\sigma}_{i-1}^{z};$		(14)
	$\displaystyle\hat{b}^{\dagger}_{i}\hat{b}_{i}$	$\displaystyle=1/2+\hat{\sigma}_{i}^{z}.$		(15)

The interactions between $a$- and $b$-atoms can then be cast into the following form

$$U\sum_{i}\hat{a}^{\dagger}_{i}\hat{a}_{i}\hat{b}^{\dagger}_{i}\hat{b}_{i}=-U\sum\limits_{i}(\hat{\sigma}^{z}_{i}\hat{\sigma}^{z}_{i+1}-1/4)\,.$$

(16)

Hence, the total Hamiltonian of this system reads

$$\hat{H}=\sum\limits_{i}\left(-4g\hat{\tau}_{i}^{z}\hat{\sigma}_{i}^{x}-2h\hat{\tau}_{i}^{x}-4\hat{\sigma}_{i}^{z}\hat{\sigma}_{i+1}^{z}\right).$$

(17)

Here and hereafter, we set $U/4=1$ as our energy unit, and $g$ and $h$ here are actually dimensionless numbers $4g/U$ and $4h/U$ in terms of the original parameters. Obviously, this Hamiltonian is invariant under the transformation $\hat{\sigma}_{i}^{x}\rightarrow-\hat{\sigma}_{i}^{x}$ and $\hat{\tau}_{i}^{z}\rightarrow-\hat{\tau}_{i}^{z}$, and possesses a local $\mathbb{Z}_{2}$ gauge symmetry.

For later convenience, we shall rotate $\hat{\bm{\tau}}_{i}$ as follows:

$$\hat{\tau}^{x}_{i}\rightarrow\hat{\tau}^{z}_{i},\ \ \hat{\tau}^{y}_{i}\rightarrow-\hat{\tau}^{y}_{i},\text{ and }\hat{\tau}^{z}_{i}\rightarrow\hat{\tau}^{x}_{i}\,.$$

(18)

The Hamiltonian then becomes

$$\hat{H}=\sum\limits_{i}\left(-4g\hat{\tau}_{i}^{x}\hat{\sigma}_{i}^{x}-2h\hat{\tau}_{i}^{z}-4\hat{\sigma}_{i}^{z}\hat{\sigma}_{i+1}^{z}\right)\,.$$

(19)

In this new basis, one can easily see that $\hat{Q}_{i}=4\hat{\sigma}_{i}^{z}\hat{\tau}_{i}^{z}$ commutes with $\hat{H}$ and its eigenvalue $q_{i}$ takes the value of $+1$ or $-1$. For any eigenstate $|\rangle$, we have $\hat{Q}_{i}|\rangle=q_{i}|\rangle$ and $q_{i}$ can be understood as the conserved local gauge charge associated with the local gauge symmetries. The four local bases $|\sigma_{i}^{z}\tau_{i}^{z}\rangle$ of site $i$ fall into two classes: $|\uparrow\uparrow\rangle$ and $|\downarrow\downarrow\rangle$ with $q_{i}=1$, and $|\uparrow\downarrow\rangle$ and $|\downarrow\uparrow\rangle$ with $q_{i}=-1$. To make use of the reduced local Hilbert space for given $q_{i}$, we introduce another spin-$1/2$ operator $\hat{\bm{\Gamma}}_{i}$ defined as

$$\hat{\Gamma}_{i}^{x}=2\hat{\tau}^{x}_{i}\hat{\sigma}_{i}^{x},\quad\hat{\Gamma}_{i}^{y}=2\hat{\tau}^{x}_{i}\hat{\sigma}^{y}_{i},\quad\text{and}\quad\hat{\Gamma}_{i}^{z}=\hat{\sigma}_{i}^{z}.$$

(20)

Indeed, these operators satisfy the $SU(2)$ algebra: $[\hat{\Gamma}_{i}^{x},\hat{\Gamma}_{i}^{y}]=i\hat{\Gamma}_{i}^{z}$, $[\hat{\Gamma}_{i}^{y},\hat{\Gamma}_{i}^{z}]=i\hat{\Gamma}_{i}^{x}$ and $[\hat{\Gamma}_{i}^{z},\hat{\Gamma}_{i}^{x}]=i\hat{\Gamma}_{i}^{y}$ (note $\hbar=1$). It can be verified that the spin-up and spin-down states, denoted as $|\Uparrow\rangle$ and $|\Downarrow\rangle$, of the $\hat{\bm{\Gamma}}_{i}$ operator, are $|\uparrow\uparrow\rangle$ and $|\downarrow\downarrow\rangle$ respectively for subspace with $q_{i}=1$, and are $|\uparrow\downarrow\rangle$ and $|\downarrow\uparrow\rangle$ respectively for subspace with $q_{i}=-1$. Noticing that

$$\hat{\tau}^{z}_{i}=4\hat{\tau}^{z}_{i}\hat{\sigma}^{z}_{i}\hat{\sigma}^{z}_{i}=\hat{Q}_{i}\hat{\Gamma}^{z}_{i}\,,$$

(21)

we can always replace $\hat{\tau}_{i}$ by $q_{i}\hat{\Gamma}_{i}$ when acting the Hamiltonian on eigenstates. With the help of the definition of the $\hat{\bm{\Gamma}}_{i}$ operators, the Hamiltonian can be simplified as

$$\hat{H}=\sum\limits_{i}\left(-2g\hat{\Gamma}_{i}^{x}-2hq_{i}\hat{\Gamma}_{i}^{z}-4\hat{\Gamma}_{i}^{z}\hat{\Gamma}_{i+1}^{z}\right).$$

(22)

Because $q_{i}$ takes $\pm 1$ in different symmetry sectors, $hq_{i}$ can be either $+h$ or $-h$ and acts as a $\mathbb{Z}_{2}$ random field. In this way, we map our model to a transverse Ising model in the presence of a disordered longitudinal field. This mapping clearly reveals that the local conserved gauge charges play the role of disordered potentials. Since the gauge charge is always quantized, the difference between this model and usual models of MBL is the distribution of the random field. In the usual models of MBL, the field takes a continuous distribution, say, uniformly distributed in $[-h,h]$, whereas in our effective model, the field takes two discrete values with equal probabilities.

Nevertheless, this difference does not exclude the possibility of MBL. In the weak disorder limit $h\ll g$, the system is ergodic and obeys the eigenstate thermalization hypothesis (ETH) ETH1 ; ETH2 ; ETH3 ; Rigol:AIP . In the opposite strong disorder limit $h\gg g$, ETH fails and the system is expected to be in the MBL regime MBL1 ; MBL2 . Below we will demonstrate the presence of both the ETH (thermal) regime and the MBL regime in this model from several different metrics.

We first focus on the eigenstate properties. It is known that in the ETH regime the energy level spacing obeys the Wigner–Dyson distribution and in the MBL regime the energy level spacing obeys the Poisson distribution Rigol:AIP . To determine the specific distribution, a dimensionless quantity $r$ has been introduced r

$$r=\langle r_{n}\rangle=\left\langle\frac{\min(\delta_{n},\delta_{n-1})}{\max(\delta_{n},\delta_{n-1})}\right\rangle,$$

(23)

where $n$ is the ascending energy level index, $\delta_{n}=E_{n+1}-E_{n}$ is the energy level spacing between two consecutive energy levels, and the average, in our case, is taken over all energy levels, including a first average over all eigenstates in a given symmetry sector and then a second average over different symmetry sectors. It can be shown that $r\approx 0.53$ for the Wigner–Dyson distribution and $r\approx 0.39$ for the Poisson distribution r .

We use the exact diagonalization method with the open boundary condition to compute all eigenvalues and determine the value of $r$. As can been seen from Fig. 3, the general trend is that when $g$ becomes larger, the value of $r$ approaches $0.53$, and when $g$ is small, the value of $r$ is around $0.39$. By comparing Fig. 3 (a) and (b), one can see that when $h$ is larger, it also requires a larger $g$ to reach the ETH regime where $r\approx 0.53$. This agrees with our expectation that the system is driven from the ETH regime to the MBL regime as the relative strength of disorder increases.

How the von Neumann entanglement entropy $S^{v}$ of a generic eigenstate with non-zero energy density scales with the size of the subsystem, $L_{\text{A}}$, is another way to differentiate the ETH regime from the MBL regime. The ETH ansatz implies the reduced density matrix of a small subsystem $A$ approaches its thermal density matrix at the temperature $T$ fixed by the energy of the eigenstate, leading to the volume law for $S^{v}$ Abanin:RMP ; Huse:ARCMP . In contrast, a MBL system fails to thermalize, and the entanglement entropy obeys the area law Huse:ARCMP ; Abanin:RMP .

The above difference is schematically shown in Fig. 2, and the behavior of $S^{v}$ is also verified in our numerical studies. As can be seen from Fig. 4(a), when $h$ is small compared with $g$, the entanglement entropy increases almost linearly with $L_{\text{A}}$ until $L_{\text{A}}$ approaches half of the system size, confirming the volume law. In contrast, as shown in Fig. 4(b), when $h$ is large compared with $g$, the entanglement entropy tends to quickly saturates as $L_{\text{A}}$ increases, demonstrating the area law.

Aside from eigenstate properties, the differences between the ETH regime and the MBL regime also manifest themselves in quench dynamics. Starting with a product state of different sites, the entanglement entropy for subsystem $A$ will increase from zero. In the ETH regime, the entanglement entropy increases linearly until it saturates linear_growth . In the MBL regime, following a short linear growth, the entanglement entropy will undergo a slow logarithmically increase in an extended time region log_growth1 ; log_growth2 ; log_growth3 ; log_growth4 . This logarithmic increase is often seen as the hallmark of the MBL phase and can be explained by the phenomenological theory of MBL phenomenology_MBL1 ; phenomenology_MBL2 . Yet this linear versus logarithmic growth of entropy is also directly related to the exponential versus power-law behavior of the out-of-time-ordered correlation Zhai_OTOC ; OTOC1 ; OTOC2 ; OTOC3 ; OTOC4 ; OTOC5 .

Our initial product state is prepared with all $\bm{\tau}_{i}$ spins polarized along the $x$-direction and all $\bm{\sigma}_{i}$ spins aligned upwards or downwards in a domain-wall fashion, as illustrated in the inset of Fig. 5(b) and Fig. 6(b). In terms of the original $\bm{\tau}_{i}$ spins before rotation, Eq. (18), all $\bm{\tau}_{i}$ spins are actually polarized along the $z$-direction. That is to say all the $f$- and $d$-atoms are prepared in the left wells, and all the $a$- and $b$-atoms are localized either in the left ($\sigma_{i}^{z}$ upwards) or right wells ($\sigma_{i}^{z}$ downwards). The advantage of choosing such an initial state is twofold. Firstly, the wave function is an equal weight superposition of states from all symmetry sectors. Secondly, the initial state has a relative high energy which mitigates the problem of limited accessible system sizes in numerical simulations.

We numerically study the growth of both the von Neumann entanglement entropy $S^{v}$ and the second-order Rényi entropy $S^{R}$ of the subsystem $A$ that is chosen as half of the entire system. The presence of both the ETH regime and the MBL regime is demonstrated using two representative parameter sets, $(g,h)=(1.0,0.2)$ and $(g,h)=(0.2,1.0)$, as shown in Fig. 5 and Fig. 6 respectively. From Fig. 5(a), one can see that for the former case where the disorder ($h$-term) is relatively weak, two entanglement entropies both increase linearly in time, as indicated by the red dashed line. In contrast, for the latter case where the disorder is relatively strong, a logarithmic growth of both entanglement entropies is revealed, as demonstrated by the red dashed lines in Fig. 6(a). The characteristic saturation time is also much longer than the former case.

Finally, we study the time evolution of physical observables following a quench. If a system is thermal, all local physical observables will evolve toward their thermal equilibrium values. For a MBL system, local physical observables, however, keep part of the memory of the initial state and do not necessarily approach the thermal equilibrium values. This criterion has been applied in identifying MBL phases in cold atom experiments Bloch_MBL .

Here we consider a special observable, the staggered magnetization, whose definition depends on the initial state,

$$\hat{m}_{s}=\frac{2}{L}\sum_{i}\xi_{i}\hat{\sigma}_{i}^{z},$$

(24)

where $\xi_{i}=1$ if initially $\sigma_{i}^{z}=1/2$ and $\xi_{i}=-1$ if initially $\sigma_{i}^{z}=-1/2$. The advantage of defining such an observable is also twofold. First, initially, the expectation value of our observable is always $m_{s}=1$. Secondly, it is easy to show that the thermal average of $\hat{m}_{s}$ is always zero because of the $\hat{\sigma}_{i}^{z}\rightarrow-\hat{\sigma}_{i}^{z},q_{i}\rightarrow-q_{i}$ symmetry of the Hamiltonian, under which $m_{s}$ changes its sign. By comparing Fig. 5(b) and Fig. 6(b), it is clear that the system thermalizes in the former case but fails to thermalize in the latter case.