Two dimensional kicked quantum Ising model: dynamical phase transitions

Quantum dynamics of strongly interacting quantum systems is one of the most exciting fields of current physics research, in particular because many fundamental phenomena, such as thermalization and equilibration in large closed systems [1], are still lacking fundamental understanding. Moreover, as such models are very difficult to simulate by the best contemporary computers due to exponential growth of Hilbert space dimension [2], theoreticians have very little predictive power, as the theory of non-equilibrium quantum thermodynamics is only beginning to emerge [3, 4].

On the mathematical physics side, the key issue is a precise understanding of the notion of quantum ergodicity, and the ergodic to non-ergodic transitions in the thermodynamic limit. By ergodicity we mean that, for most observables, the time average of an observable and its time correlation will coincide with the canonical average of the observable alone. For periodically driven quantum systems, this means that starting from almost any initial state, all observables and correlations are determined by an infinite temperature Gibbs state. Lack of ergodicity thereof implies a sort of localization within a part of many-body Hilbert space. This type of many-body localization is essentially different from and perhaps more mysterious than the one studied in many-body systems with onsite disorder [5–8].

Some conceptual and computational attempts in this direction have been made in [9–15]. In particular, one showed [9–12] that kicked XX-Z chains and kicked Ising (KI) spin chains typically exhibit non-ergodic to ergodic transition with an increasing kicking period, despite the fact that the limiting (autonomous) system when the driving period goes to zero may be integrable. This transition was accompanied with the transition in quasi-energy level spacing distribution going from Poissonian statistics of uncorrelated levels for integrable regimes to Wigner-like level spacing distribution accurately describing Gaussian orthogonal (or unitary) ensembles of random matrix theory (GOE) for non-integrable regimes. The fact that in the semi-classical limit of the small effective Planck constant (if the latter can be meaningfully defined), both transitions (in the spectral correlations and the ergodicity of the system) are in one-to-one correspondence [16], has led to an intuitive belief that should hold in more general cases, such as many-body quantum systems without a small ℏ parameter. This, however, has not been carefully explored yet. A third measure of quantum chaos or quantum ergodicity in periodically driven quantum systems can be introduced, referred to as Loschmidt echo, survival amplitude, or equivalently, (Fourier transform of) the quasienergy spectral density. We shall refer to it as the spectral density, and for ergodic kicked quantum systems, one may expect it to be constant. Recently it has been suggested that for quantized chaotic single particle periodically kicked systems (e.g., for the so-called kicked top), one obtains dynamical instabilities signalled by discontinuous (sharp) transitions, or singularities in the spectral density. We shall argue later that spectral density is typically always constant in kicked locally interacting quantum spin $1/2$ chains, so there could be no transitions, but an interesting question arises for what happens in higher-dimensional spin lattices.

In this paper we study what is likely the simplest non-trivial dynamical many-body model in two dimensions, namely the two-dimensional (2-D) version of the KI spin $1/2$ model introduced in [12]. We numerically accurately calculate different dynamical and spectral properties of the model on a finite periodic rectangular ${{L}_{x}}\times {{L}_{y}}$ lattices implementing local one and two qubit operations on graphics processing units (GPUs) of a desktop computer, in a spirit similar to [17]. Namely we compute phase diagrams (dependencies on the modelʼs parameters) of dynamical susceptibilities, spectral densities, and quasi-energy level spacing distributions. We find that in most parts of parameter space these quantities turn out to be stable against increasing lattice sizes, so we formulate certain conjectures about the thermodynamic behavior. We find, first, similar to one-dimensional (1-D) chains, well-defined transitions between non-ergodic and ergodic dynamical susceptibilities of local observables. Second, we find non-trivial transitions between flat spectral density and spectral densities with shapes essentially determined by a single Fourier mode, $\rho (\phi )=\frac{1}{2\pi }+c{\rm cos} (k\phi )$, for some constant c and integer k. The Fourier coefficient c seems to decrease with increasing lattice size, but for any fixed Hilbert space dimension, one finds a dramatic difference of $|\rho (\phi )-\frac{1}{2\pi }|$ from a prediction of a flat spectrum (say of a random unitary matrix). However, remarkably and surprisingly, the transition points of dynamical susceptibilities do not correspond to the points of transitions of level density. Third, analysis of level spacing distributions of properly unfolded [16] quasienergy spectra reveals universal GOE (Wigner) statistics across most of the parameter space, whereas non-universal statistics are only observed on singular parameter regions with trivially integrable dynamics. This suggests there are no non-trivial integrable points. Therefore the transition in spectral correlations is instantaneous in thermodynamic limit and does not correspond with any other measure of quantum ergodicity, such as dynamical susceptibilities and spectral densities. We believe this is a remarkable observation that can question the versatility of level spacing distribution in quantum many-body systems. In particular, we believe other direct measures of quantum ergodicity should be used when discussing thermalization or equilibration.

Our paper is organized as follows. In section 2, we introduce the model and discuss its general properties, including its symmetries. In section 3, we comment on the spectral density and observe its properties for all model parameter values. A first picture of the model emerges. In section 4, we discuss the correlation properties of the spectra, using the nearest neighbour spacing distribution. The second picture emerges. In section 5, we study the dynamical susceptibilities (autocorrelation functions of certain observables) where a third picture of the model arises. Finally, we gather the results to come up with our conclusions in section 6.

We study a periodic 2-D lattice of kicked spin $1/2$ particles, inspired by the KI chain, proposed by one of the authors in [12]. The particle at site $m\in \{1,\ldots ,{{L}_{x}}\}$, $n\in \{1,\ldots ,{{L}_{y}}\}$ will be described by standard Pauli operators $\sigma _{m,n}^{\alpha }$, $\alpha \in \{x,y,z\}$.

Let us start by defining a 2-D Ising Hamiltonian

$$\begin{eqnarray}&&{{H}_{1}}=J{{H}_{{\rm I}}},\quad {{H}_{{\rm I}}}=\mathop{\sum }\limits_{m=0}^{{{L}_{x}}-1}\mathop{\sum }\limits_{n=0}^{{{L}_{y}}-1}(\sigma _{m,n}^{z}\sigma _{m+1,n}^{z}+\sigma _{m,n}^{z}\sigma _{m,n+1}^{z}),\end{eqnarray} \tag{ 1 }$$

with periodic boundary conditions $\sigma _{m,{{L}_{y}}}^{\alpha }\equiv \sigma _{m,0}^{\alpha }$, $\sigma _{{{L}_{x}},n}^{\alpha }\equiv \sigma _{0,n}^{\alpha }$. We now define a Zeeman Hamiltonian for a spatially homogeneous magnetic field $\vec{b}$

$$\begin{eqnarray}&&{{H}_{0}}=\mathop{\sum }\limits_{m=0}^{{{L}_{x}}-1}\mathop{\sum }\limits_{n=0}^{{{L}_{y}}-1}\vec{b}\cdot {{\vec{\sigma }}_{m,n}}=\vec{b}\cdot \vec{S},\quad \vec{S}=:\mathop{\sum }\limits_{m=0}^{{{L}_{x}}-1}\mathop{\sum }\limits_{n=0}^{{{L}_{y}}-1}{{\vec{\sigma }}_{m,n}}.\end{eqnarray} \tag{ 2 }$$

Notice that we can always choose the coordinate system such that $\vec{b}=({{b}_{x}},0,{{b}_{z}})$, so that both H₀ and H₁ are real. We will normally consider only a transverse field, that is, $\vec{b}=({{b}_{x}},0,0)$. The parameters J (interspin interaction) and b_x (transverse magnetic field) are independent dimensionless parameters that specify the model.

We consider a time-dependent Hamiltonian, where the magnetic field is modulated by periodic $\delta -$ pulses of period τ

$$\begin{eqnarray}&&H(t)={{H}_{1}}+{{H}_{0}}\mathop{\sum }\limits_{j\in \mathbb{Z}}\delta (t-j\tau ).\end{eqnarray} \tag{ 3 }$$

One-step quantum evolution propagator for the KI model over one period of the model $U(t)=\mathcal{T}{\rm exp} (-{\rm i}\int _{{{0}^{+}}}^{{{1}^{+}}}{\rm d}tH(t))$—the so-called Floquet map—reads, setting $\tau =1$ by a free choice of units:

$$\begin{eqnarray}&&{{U}_{{\rm KI}}}={{U}_{{\rm Ising}}}(J){{U}_{{\rm kick}}}(\vec{b}),\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&{{U}_{{\rm Ising}}}(J)={\rm exp} (-{\rm i}{{H}_{1}})={\rm exp} (-{\rm i}J{{H}_{{\rm I}}}),\quad {{U}_{{\rm kick}}}(\vec{b})={\rm exp} (-{\rm i}{{H}_{0}})={\rm exp} (-{\rm i}\vec{b}\cdot \vec{S}).\end{eqnarray} \tag{ 5 }$$

2.1. Symmetries

We shall now briefly discuss some obvious symmetries of the model that help reduce computational complexity of simulations.

Parameter space symmetries

Notice that the system is periodic in the parameters since the spectra of operators ${{H}_{{\rm I}}}$ and $\hat{b}\cdot \vec{S}$, where $\hat{b}=\vec{b}/|\vec{b}|$ is the unit vector in the direction of $\vec{b}$, form subsets of integers. In fact, in the spectrum of ${{H}_{{\rm I}}}$ there can be only integers with fixed remainder of division with 4, since flipping an arbitrary spin can change ${{H}_{{\rm I}}}$ only by $\pm 4$ or $\pm 8$. More precisely, the spectrum of ${{H}_{{\rm I}}}$ consists of points $\{2{{L}_{x}}{{L}_{y}},2{{L}_{x}}{{L}_{y}}-4,2{{L}_{x}}{{L}_{y}}-8,\ldots ,-2{{L}_{x}}{{L}_{y}}\}$, hence

$$\begin{eqnarray}&&{{U}_{{\rm Ising}}}(J+\pi /2)={{(-1)}^{{{L}_{x}}{{L}_{y}}}}{{U}_{{\rm Ising}}}(J).\end{eqnarray} \tag{ 6 }$$

Similarly, the spectrum of $\hat{b}\cdot \vec{S}$ consists of points $\{{{L}_{x}}{{L}_{y}},{{L}_{x}}{{L}_{y}}-2,{{L}_{x}}{{L}_{y}}-4,\ldots ,-{{L}_{x}}{{L}_{y}}\}$, hence

$$\begin{eqnarray}&&{{U}_{{\rm kick}}}(\vec{b}+\pi \hat{b})={{(-1)}^{{{L}_{x}}{{L}_{y}}}}{{U}_{{\rm kick}}}(\vec{b}).\end{eqnarray} \tag{ 7 }$$

Let us now further assume that the field is transverse $\hat{b}=(1,0,0)$ as will be the case for most of this paper. Then, performing a checkerboard canonical (unitary) transformation $\hat{C}$, namely: flipping the signs of $y,z$ components of spins ${{\vec{\sigma }}_{m,n}}$ for all even $m+n$ (i.e., rotating for angle πalong the x-axis), one finds that

$$\begin{eqnarray}&&{{U}_{{\rm Ising}}}(-J){{U}_{{\rm kick}}}(\vec{b})={{\hat{C}}^{\dagger }}{{U}_{{\rm Ising}}}(J){{U}_{{\rm kick}}}(\vec{b})\hat{C}.\end{eqnarray} \tag{ 8 }$$

Similarly, canonical transformation $\hat{D}$, which flips $x,y$ components of all spins (rotates around z axis for angle π), yields

$$\begin{eqnarray}&&{{U}_{{\rm Ising}}}(J){{U}_{{\rm kick}}}(-\vec{b})={{\hat{D}}^{\dagger }}{{U}_{{\rm Ising}}}(J){{U}_{{\rm kick}}}(\vec{b})\hat{D}.\end{eqnarray} \tag{ 9 }$$

Therefore, changing the sign of J or b_x leaves invariant all physical properties of the model, in particular the spectrum of ${{U}_{{\rm KI}}}$, so the principal domain of the phase diagram of the transverse field KI model only consists of a rectangle $(J,{{b}_{x}})\in [0,\pi /4]\times [0,\pi /2]$.

Symmetry reduction of the Hilbert space

In the general case, the system is symmetric under the following geometric operations, generating the symmetry group of the model: reflexion over the horizontal axis ${{R}_{x}}$, reflexion over the vertical axis ${{R}_{y}},$ horizontal translation ${{T}_{x}},$ and vertical translation ${{T}_{y}}$.

To illustrate these symmetries, let us number the sites in a $({{L}_{x}}=4)\times ({{L}_{y}}=3)$ grid from left to right and bottom to top, and consider a state of the computational basis $|{{\psi }_{0}}\rangle =|{{i}_{0}}{{i}_{1}}\cdots {{i}_{11}}\rangle$, with ${{i}_{j}}\in \{0,1\}$. The action of ${{R}_{x}}$ on the grid will be to transform it into

$$\begin{eqnarray}\;\begin{array}{ccccccccccccccc} 8 & 9 & 10 & 11 \\ 4 & 5 & 6 & 7 \\ 0 & 1 & 2 & 3 \\ \end{array}\qquad \to \qquad \begin{array}{ccccccccccccccc} 0 & 1 & 2 & 3 \\ 4 & 5 & 6 & 7 \\ 8 & 9 & 10 & 11 \\ \end{array},\end{eqnarray} \tag{ 10 }$$

so

$$\begin{eqnarray*}&&{{R}_{x}}\mid {{\psi }_{0}}\rangle =\mid {{i}_{8}}\cdots {{i}_{11}}{{i}_{4}}\cdots {{i}_{7}}{{i}_{0}}\cdots {{i}_{3}}\rangle .\end{eqnarray*}$$

The action of the reflection ${{R}_{y}}$ is similar:

$$\begin{eqnarray}\;\begin{array}{ccccccccccccccc} 8 & 9 & 10 & 11 \\ 4 & 5 & 6 & 7 \\ 0 & 1 & 2 & 3 \\ \end{array}\qquad \to \qquad \begin{array}{ccccccccccccccc} 11 & 10 & 9 & 8 \\ 7 & 6 & 5 & 4 \\ 3 & 2 & 1 & 0 \\ \end{array},\end{eqnarray} \tag{ 11 }$$

and the action over a member of the computational basis is

$$\begin{eqnarray*}&&{{R}_{y}}\mid {{\psi }_{0}}\rangle =\mid {{i}_{3}}\cdots {{i}_{0}}{{i}_{7}}\cdots {{i}_{4}}{{i}_{11}}\cdots {{i}_{8}}\rangle .\end{eqnarray*}$$

Similarly, we can picture the effects of the translations. The effect of vertical and horizontal translations on the original grid are

$$\begin{eqnarray}\;\begin{array}{ccccccccccccccc} 8 & 9 & 10 & 11 \\ 4 & 5 & 6 & 7 \\ 0 & 1 & 2 & 3 \\ \end{array}\to \begin{array}{ccccccccccccccc} 4 & 5 & 6 & 7 \\ 0 & 1 & 2 & 3 \\ 8 & 9 & 10 & 11 \\ \end{array},\qquad \begin{array}{ccccccccccccccc} 8 & 9 & 10 & 11 \\ 4 & 5 & 6 & 7 \\ 0 & 1 & 2 & 3 \\ \end{array}\to \begin{array}{ccccccccccccccc} 11 & 8 & 9 & 10 \\ 7 & 4 & 5 & 6 \\ 3 & 0 & 1 & 2 \\ \end{array},\end{eqnarray} \tag{ 12 }$$

respectively. Thus, the action of the symmetry is simply

$$\begin{eqnarray*}\begin{array}{rcl} {{T}_{x}}\mid {{\psi }_{0}}\rangle & = & \mid {{i}_{8}}\cdots {{i}_{11}}{{i}_{0}}\cdots {{i}_{7}}\rangle ,\;\;{\rm and} \\ {{T}_{y}}\mid {{\psi }_{0}}\rangle & = & \mid {{i}_{3}}{{i}_{0}}{{i}_{1}}{{i}_{2}}{{i}_{7}}{{i}_{4}}{{i}_{5}}{{i}_{6}}{{i}_{11}}{{i}_{8}}{{i}_{9}}{{i}_{10}}\rangle . \\ \end{array}\end{eqnarray*}$$

Notice that $R_{y}^{2}=R_{x}^{2}=T_{y}^{{{L}_{y}}}=T_{x}^{{{L}_{x}}}={\mathbb{1}}$. It can also be noted that ${{T}_{y}}={{T}^{{{L}_{x}}}}$, where T is the translation operator that acts as ${{T}_{|}}{{\psi }_{0}}\rangle =|{{i}_{11}}{{i}_{1}}{{i}_{2}}\cdots {{i}_{10}}\rangle$. The symmetry subspaces of the Hilbert space are therefore specified by two quasi-momenta ${{k}_{x}}\in \{0,\ldots ,{{L}_{x}}-1\},{{k}_{y}}\in \{0,\ldots ,{{L}_{y}}-1\}$, and for symmetric sectors with ${{k}_{x,y}}=0$ or ${{k}_{x,y}}={{L}_{x,y}}/2$ by additional reflection signs ${{\pi }_{x}},{{\pi }_{y}}\in \{\pm 1\}$, such that the states $|\psi \rangle$ from the subspace satisfy ${{R}_{x,y}}|\psi \rangle ={{\pi }_{x,y}}|\psi \rangle ,{{T}_{x,y}}|\psi \rangle ={{{\rm e}}^{-2\pi {\rm i}{{k}_{x,y}}/{{L}_{x,y}}}}|\psi \rangle$.

When we consider the special case of a transverse magnetic field, another symmetry arises. A parity operator ${{\prod }_{m,n}}\sigma _{m,n}^{x}$ commutes with the Floquet operator: $[{{\prod }_{m,n}}\sigma _{m,n}^{x},{{U}_{{\rm KI}}}]=0$.

Our KI model also has an anti-unitary symmetry, namely if K is complex conjugation in the standard Pauli basis, then ${{K}^{-1}}\sigma _{j,k}^{z}K=\sigma _{j,k}^{z}$, ${{K}^{-1}}\sigma _{j,k}^{x}K=\sigma _{j,k}^{x}$, and ${{K}^{-1}}\sigma _{j,k}^{y}K=-\sigma _{j,k}^{y}$. Writing a symmetrized Floquet propagator $U_{{\rm KI}}^{\prime }={{{\rm e}}^{-{\rm i}{{H}_{0}}/2}}{{U}_{{\rm KI}}}{{{\rm e}}^{{\rm i}{{H}_{0}}/2}}={{{\rm e}}^{-{\rm i}{{H}_{0}}/2}}{{{\rm e}}^{-{\rm i}{{H}_{1}}}}{{{\rm e}}^{-{\rm i}{{H}_{0}}/2}}$, we then have immediately

$$\begin{eqnarray}&&{{K}^{-1}}U_{{\rm KI}}^{\prime }K={{(U_{{\rm KI}}^{\prime })}^{-1}},\quad {\rm and}\quad U_{{\rm KI}}^{\prime }={{(U_{{\rm KI}}^{\prime })}^{T}}.\end{eqnarray} \tag{ 13 }$$

Using the standard wisdom [16], the model should then—if 'quantum chaotic'–correspond to circular orthogonal ensemble (COE) of random unitary symmetric matrices.

2.2. Steady field limit

Using the same machinery, we can study the time-independent limit, corresponding to keeping $\vec{b}/J$ fixed, while letting $|J|$ go to zero. Then, the scaled Hamiltonian is simply

$$\begin{eqnarray}&&H={{H}_{1}}+{{H}_{0}}.\end{eqnarray} \tag{ 14 }$$

To study this Hamiltonian, we also used compute unified device architecture (CUDA) machinery, and used both first- and second-order Trotter approximations. For the results presented here, we verified that the first-order Trotter evolution gives essentially the same results as the second-order, meaning they the changes in the figures presented are so small that cannot be noticed.

The spectrum of the Floquet map, $\mathcal{S}=\{{{\phi }_{n}};n=1,\ldots ,\mathcal{N}:={{2}^{{{L}_{x}}{{L}_{y}}}}\}$, defined by the unitary eigenvalue problem

$$\begin{eqnarray}&&{{U}_{{\rm KI}}}\mid {{\psi }_{n}}\rangle ={{{\rm e}}^{-{\rm i}{{\phi }_{n}}}}\mid {{\psi }_{n}}\rangle ,\end{eqnarray} \tag{ 15 }$$

entails the main dynamical features of the model. The statistical properties of $\mathcal{S}$ for systems with chaotic classical limit has been the central theme of quantum chaos [16]. However, very little is known about the distribution of $\{{{\phi }_{n}}\}$ for many-body quantum models, despite the fact that full many-body quantum dynamics are becoming experimentally accessible in recent years, in particular, in cold atom laboratories [18]. Even the behavior of the simplest spectral characteristic, the one-point function or the spectral density, defined as

$$\begin{eqnarray}&&\rho (\phi )=\frac{1}{\mathcal{N}}\mathop{\sum }\limits_{n=1}^{\mathcal{N}}\delta (\phi -{{\phi }_{n}}),\end{eqnarray} \tag{ 16 }$$

would be of great interest to know. In autonomous (time-independent) quantum many-body systems, the spectral density is predicted to go to a Gaussian in the thermodynamic limit [19, 20], while for periodically driven quantum systems, one may perhaps intuitively expect (and observe, in 1-D chains [11]) that the Floquet quasi-energy spectral density would be the constant (flat) function $\rho (\phi )=\frac{1}{2\pi }$, in a generic case.

The spectral density $\rho (\phi )$ is a $2\pi$-periodic function and therefore can be represented in terms of the Fourier modes as

$$\begin{eqnarray}&&\rho (\phi )=\frac{1}{2\pi }\left( 1+2\mathop{\sum }\limits_{k=1}^{\infty }{{\rho }_{k}}{\rm cos} (k\phi ) \right),\end{eqnarray} \tag{ 17 }$$

where the Fourier coefficients ${{\rho }_{k}}$ are given as traces of the k–step KI propagator

$$\begin{eqnarray}&&{{\rho }_{k}}:=\int _{0}^{2\pi }{\rm d}\phi \;\rho (\phi ){{{\rm e}}^{{\rm i}k\phi }}=\frac{1}{\mathcal{N}}{\rm tr}U_{{\rm KI}}^{k}.\end{eqnarray} \tag{ 18 }$$

The symmetry property ${{\rho }_{-k}}={{\rho }_{k}}$ has been used, following from the symmetry of the spectra of H₀ and H₁ around zero energy and cyclicity of the trace.

The sum in equation (16) can be carried in the whole Hilbert space, or in a single symmetry sector (with fixed quasi-momenta and/or parities). We calculated the spectra directly, using a basis that splits the evolution operator into different sectors, and numerically fully diagonalizing each sector independently. That way, we could calculate the entire spectrum for sizes of up to 5 × 4. The behavior over different symmetry sectors seems to be similar in all examples we considered, see figure 1 (right panel).

Zoom In Zoom Out Reset image size

For efficient numerical computation of the leading Fourier components ${{\rho }_{k}}$, one should instead use the expression in terms of traces of powers of the propagator directly (18). The computation can be further simplified by noting that the trace over the many-body Hilbert space is self-averaging and can be approximated using an expectation value in a single typical (random) state. Namely

$$\begin{eqnarray}&&\frac{1}{\mathcal{N}}{\rm tr}U_{{\rm KI}}^{t}=\int \langle \psi |U_{{\rm KI}}^{t}|\psi \rangle {\rm d}\mu (\psi )\approx \langle {{\psi }_{{\rm random}}}\mid U_{{\rm KI}}^{t}\mid {{\psi }_{{\rm random}}}\rangle ,\end{eqnarray} \tag{ 19 }$$

where ${\rm d}\mu (\psi )$ is the measure induced by the Haar measure over the unitary group, and $|{{\psi }_{{\rm random}}}\rangle$ is a state drawn at random with the Haar measure. Moreover, if one selects that state belonging to a given symmetry subspace, one then studies the spectral properties of that particular subspace. For practical computation, one may take all components of $|\psi \rangle$ as random Gaussian $c-$numbers with zero mean and equal variance and then normalize the state. In our case, calculating the kick-by-kick evolution of a state is quite efficient, so this form of calculating the Fourier transform of the spectral density is particularly convenient.

We now examine the behavior of the spectral density for several parameter values and several sizes. The left panel of figure 1 suggests a dominant Fourier component ${{\rho }_{k}}$ of the spectral density, whose magnitude varies with the size of the system. Examining each of the Fourier contributions as a function of the parameters proves very useful. Such analysis is carried out for all coefficients up to k = 12, varying the transverse component of the magnetic field. A similar behavior is obtained for different sizes, as can be seen in figure 2. The most outstanding fact is that there are clearly two different regions. One in which we have all Fourier coefficients magnitude close to the average random value, given by $|{{\rho }_{{\rm noise}}}|$, and another region, in which there is an ordered phase, manifested by $|{{\rho }_{k}}|\gg |{{\rho }_{{\rm noise}}}|$. We estimate $|{{\rho }_{{\rm noise}}}|$, replacing $U_{{\rm KI}}^{t}$ in equation (19) by a random unitary matrix of dimension ${{2}^{{{L}_{x}}{{L}_{y}}}}$. This gives rise to $|{{\rho }_{{\rm noise}}}|\propto {{2}^{-L/2}}$ where $L={{L}_{x}}{{L}_{y}}$, which is in good agreement with the tendency observed. We have also plotted the dominant Fourier component at b_x = 0.2. Indeed, there is a small decay of the oscillations for the b_x = 0.2, however, the comparative effect enhances with the system size, in the sense that $|{{\rho }_{3}}(b=0.2)/{{\rho }_{{\rm noise}}}|\underset{\scriptscriptstyle\thicksim}{\propto }{\rm exp} (0.27L)$, thus sharpening the transition from a disordered to an ordered phase. Therefore, even though the spectral density $\rho (\phi )$ seems to universally approach a constant $1/(2\pi )$ when one approaches the thermodynamic limit ${{L}_{x}},{{L}_{y}}\to \infty$, there are discontinuous transitions on the size scaling of the deviation $\rho (\phi )-\frac{1}{2\pi }$ when system parameters are changed.

Zoom In Zoom Out Reset image size

One can obtain an interesting global picture of the model by plotting a spectral density phase diagram. Namely, we determine and plot the leading nontrivial spectral component k for which ${{\rho }_{k}}$ is dominating, as a function of modelʼs parameters. We shall consider the spectrum to be flat, k = 0, if all Fourier coefficients ${{\rho }_{k}}$, for $k\gt 0$, are comparable to ${{\rho }_{{\rm noise}}}$. That is, the spectrum is declared flat, if

$$\begin{eqnarray}&&{{\rho }_{k}}\lt 20{{\rho }_{{\rm noise}}},\quad \forall k\gt 0.\end{eqnarray} \tag{ 20 }$$

On the other hand, we consider the system to be in the phase $k\gt 0$, if $|{{\rho }_{k}}|\gt |{{\rho }_{k^{\prime} }}|$, for all $k\ne k^{\prime} \ne 0$. We note that typically a single Fourier component is dominating others by several orders of magnitude. The gap between the Fourier components even increases when we increase the lattice size (see figure 2). See figure 3 for a comparison of phase diagrams for different lattice sizes that seems remarkable stable.

Zoom In Zoom Out Reset image size

To analyze the correlation properties of the spectrum, we used the commonly studied nearest neighbor spacing distribution. That is, we consider the distribution of level spacings

$$\begin{eqnarray}&&{{\tilde{s}}_{i}}={{\phi }_{i+1}}-{{\phi }_{i}}\end{eqnarray} \tag{ 21 }$$

where ${{\phi }_{i}}$ are the sorted eigenphases of the evolution operator. To remove the effect of non-uniform level density, we perform unfolding, i.e., a smooth non-linear scaling of the eigenvalues to obtain a uniform spectral density. Since the density is typically well described by few Fourier components, we shall take six of them to numerically perform the unfolding. Thus, we shall use the mapping

$$\begin{eqnarray}&&{{\varphi }_{i}}=2\pi \left( \frac{{{\phi }_{i}}}{2\pi }+\mathop{\sum }\limits_{k=1}^{6}{{\alpha }_{k}}\frac{{\rm sin} k{{\phi }_{i}}}{k} \right)\end{eqnarray} \tag{ 22 }$$

where ${{\alpha }_{k}}={{\rho }_{k}}/\pi$ is determined numerically using direct evolution. In this way, we can unfold the spectra to obtain fairly flat distributions of unfolded level spacings. Another very important aspect that must be considered is the fact that different symmetry sectors are not statistically correlated with each other, so we must rather look at the distribution of

$$\begin{eqnarray}&&s_{i}^{K}=\varphi _{i+1}^{K}-\varphi _{i}^{K}\end{eqnarray} \tag{ 23 }$$

where $K=\{{{k}_{x}},{{k}_{y}},{{\pi }_{x}},{{\pi }_{y}},\Pi \}$ denotes the set of quantum numbers that determines the irreducible quantum sector.

By construction, the mean spacing s_i^K equals one, so the probability density P(s) of $\{s_{i}^{K}\}$ is normalized such that $\int _{0}^{\infty }P(s){\rm d}s=\int _{0}^{\infty }sP(s){\rm d}s=1$. The famous quantum chaos conjecture states that P(s) behaves generically as the corresponding classical ensemble of random matrices [16], given that the classical limit is strongly chaotic (i.e., a hyperbolic dynamical system). Given the time reversal symmetry, the corresponding ensemble, in our case, would be the COE. To a very good approximation, the level spacing distribution is given in terms of 2 × 2 random real symmetric matrices, which is the so-called Wigner surmise ${{P}_{\,{\rm Wigner}}}(s)=\frac{\pi }{2}s{\rm exp} (-\pi {{s}^{2}}/4)$. A similar conjecture has been suggested for strongly non-integrable quantum many-body systems [21], but it has not yet been established precisely how non-integrability and level statistics are related in a given class of models.

We present in the left panel of figure 4 the P(s) for three typical cases of our KI model, each of which behaves differently with respect to the dynamical ergodicity measures discussed in this paper (spectral density and dynamical susceptibility). We see, however, that P(s) is in all three cases are excellently described by COE or Wigner surmise. In a finer scale (inset of the left panel of figure 4), even the difference between Wignerʼs surmise and the exact COE result can be resolved for the dynamical data. In the right panel of figure 4, we plot the Kolmogorov distance ${{\int }_{\mathcal{X}}}|f(x)-g(x)|{\rm d}x$ between the observed distribution of nearest neighbor spacings and the random matrix theory (RMT) prediction. There is good agreement with COE in the entire parameter space except for trivial integrable cases of zero field or zero spin interaction (both modulo $\pi /2$), or specially commensurate fields where the spectrum of ${{U}_{{\rm KI}}}$ can be explicitly computed in terms of regular or number-theoretic functions.

Zoom In Zoom Out Reset image size

So far we have analyzed dynamical properties of the 2-D KI model that depend solely on its spectrum. Now we shall focus on dynamical correlations of local observables, which are the key inputs to any linear response treatment of condensed matter theory [22–24].

Consider a traceless observable M (typically extensive and local), say a component of magnetization $M={{S}^{\nu }}$, $\nu \in \{x,y,z\}$. We define its time-autocorrelation with respect to the KI dynamics as

$$\begin{eqnarray*}&&C(t)=\frac{1}{{{L}_{x}}{{L}_{y}}\mathcal{N}}{\rm tr}(M(t)M)\approx \frac{1}{{{L}_{x}}{{L}_{y}}}\langle {{\psi }_{{\rm random}}}\mid M(t)M\mid {{\psi }_{{\rm random}}}\rangle .\end{eqnarray*}$$

If ${\rm tr}M\ne 0$, the corresponding constant has to be subtracted from C(t). Here M(t) denotes the time-dependent observable in the Heisenberg picture: $M(t)=U_{{\rm KI}}^{-t}MU_{{\rm KI}}^{t}$, $t\in \mathbb{Z}$.

One measures the ergodicity of an observable M by the so-called dynamical susceptibility, defined as the time average of C(t):

$$\begin{eqnarray}&&{{\chi }_{M}}=\mathop{{\rm lim} }\limits_{T\to \infty }\frac{1}{T}\mathop{\sum }\limits_{t=1}^{T}C(t).\end{eqnarray} \tag{ 24 }$$

By definition, observable M is ergodic with respect to dynamics ${{U}_{{\rm KI}}}(t)$, if ${{\chi }_{M}}=0$, and non-ergodic otherwise. Note that ${{\chi }_{M}}$ is always non-negative as it represents the spectral weight, i.e., the power spectrum of M(t), at frequency $\omega =0$.

For numerical investigations, to diminish the transient effects of relaxation, it is useful to define a finite time average between two, sufficiently large times ${{T}_{1}},{{T}_{2}}\in \mathbb{Z}$, as

$$\begin{eqnarray}&&\langle C(t)\rangle _{{{T}_{1}}}^{{{T}_{2}}}=\frac{1}{{{T}_{2}}-{{T}_{1}}+1}\mathop{\sum }\limits_{t={{T}_{1}}}^{{{T}_{2}}}C(t).\end{eqnarray} \tag{ 25 }$$

To illustrate the ergodic properties of the model for different parameters, we have analyzed a series of observables. We used both observables symmetric under particle permutation and non-symmetric observables, but restricted ourselves to sums of few-site local observables.

We studied the general case with $M=\hat{n}\cdot \vec{S}$ with $\hat{n}$ an arbitrary unit vector. Small system sizes revealed that when $\hat{b}$ points in any of the three axis directions (x, y, or z), the dynamical susceptibilities are exactly symmetric in parameter space with respect to the line ${{b}_{x}}=\pi /4$ in parameter space. However, this is the case for general observables, although there is a strong tendency to be exactly symmetric. Observables that are more general also have this tendency, which becomes increasingly more difficult to explore for moderate large systems.

We found three qualitatively different kinds of behavior of dynamical susceptibility, exemplified in figure 5. Here we compare the behavior for $M={{S}^{x}}$ for different sizes of the system and a fixed Ising interaction of J = 0.5. For b_x = 0.2, there is no decay for large times or dimensions. There seems to be a non-vanishing asymptotic value ${{\chi }_{M}}\ne 0$, which is also characteristic of integrable systems. For b_x = 0.3, C(t) seems to decay algebraically to an asymptotic value, which decays with increasing lattice size, so it is reasonable to conclude ${{\chi }_{M}}=0$. Finally, there are parameter values (say b_x = 0.5) for which fluctuations on top of an asymptotic value are reached exponentially fast. Again, the asymptotic value decreases with the Hilbert space dimension, as well as the fluctuations, suggesting they both vanish in the thermodynamic limit. In fact, the sources of data fluctuations at large times are two-fold: finite size effects and random initial state sampling (approximating the trace), whereas empirical evidence suggests that the latter (contributing to fluctuations as $\sim 1/\sqrt{{{2}^{{{L}_{x}}{{L}_{y}}}}}$) quickly becomes negligible.

Zoom In Zoom Out Reset image size

In left panel of figure 6, we observe the dynamical susceptibility as a function of both the transverse field b_x and J (here we set b_z = 0). There is clearly a set of parameters for which the model is not ergodic, that is, where ${{\chi }_{M}}\ne 0$. However, for J = 0.5, there seems to be a range of b_x where the correlations clearly vanish, namely for $0.4\lt {{b}_{x}}\lt 1.2$. We have observed also that as we were able to increase the number of particles, the transition was increasingly sharper. In the right panel of figure 6, we sketch the full three-dimensional (3-D) phase diagram (of order parameter ${{\chi }_{M}}$) in the parameter space $(J,{{b}_{x}},{{b}_{z}})$, clearly indicating distinct regions of ergodic and non-ergodic dynamics.

Zoom In Zoom Out Reset image size

Comparing these data to phase diagrams of level density figure 1, or nearest neighbor agreement with RMT figure 4 one finds that there is no point to point correspondence between different regimes in the parameter space. One can have a flat or non-flat level density in either an ergodic or non-ergodic regime for the dynamics of local observables. This speculative conclusion is certainly surprising and calls for a deeper understanding of the role of locality of observables in long time dynamics.

We have also used the same program a glimpse the behavior of the steady field model, and found hints that this is also a rich model, in which both situations of ergodic and non-ergodic dynamics are found for different parameter values, see figure 7.

Zoom In Zoom Out Reset image size

In this paper we describe a computational excursion into ergodic properties of 2-D periodically driven quantum spin systems. In the absence of efficient computational techniques, we implemented brute force simulation of the systemʼs dynamics. Speculating on the thermodynamic properties of the system by inspecting an increasingly large sequence of periodic lattices, our results suggest several rather intriguing conclusions. The spectral density of the Floquet operator displays phase transitions from regions of flat density to regions with non-trivial spectral densities dominated by non-zero Fourier components. Local observables display ergodic regimes with decaying correlations and non-ergodic regimes with non-decaying correlations, which do not correspond to regions of flat versus non-flat level densities. Moreover, the level spacing distribution is essentially given by Wigner surmise of random matrix theory over the entire parameter space, where the model is non-integrable, and therefore, surprisingly, does not provide any useful information on the systemʼs ergodicity. We believe that our numerical results generate a strong motivation for further theoretical investigations into dynamics of periodically driven (or discrete-time) interacting spin models on 2-D lattices.

TP acknowledges financial support from grant P1–0044 and J1–5439 of the Slovenian Research Agency. Support from CONACyT 153190 and UNAM-PAPIIT IA101713 is acknowledged by CP and EV.

The GPU implementation is made using Nvidia CUDA architecture for GPU parallel computing. We store the coefficients (in the computational basis) of the state one wishes to evolve on the global memory of the GPU and then apply the required quantum gates (4) on it. The parallelization is done, realizing that the application of a n-qubit gate can be decomposed in ${{2}^{L-n}}$ independent parallel operations. Each is completed by a single thread in the GPU.

Using the threads index threadIdx.x+blockIdx.x*blockDim.x and the computational base, we specify the entries on the state a specific thread will compute on; for example, when applying a one-qubit gate on the second qubit, each thread shall act on the coefficients of the components. In particular,

$$\begin{eqnarray*}\begin{array}{lc} {\rm thread}\ 0: & |\cdots 0\underline{0}00\rangle ,|\cdots 0\underline{1}00\rangle \\ {\rm thread}\ 1: & |\cdots 0\underline{0}01\rangle ,|\cdots 0\underline{1}01\rangle \\ {\rm thread}\ 2: & |\cdots 0\underline{0}10\rangle ,|\cdots 0\underline{1}10\rangle \\ \vdots & {} \\ {\rm thread}\;{{2}^{L-1}}: & |\cdots 1\underline{0}11\rangle ,|\cdots 1\underline{1}11\rangle . \\ \end{array}\end{eqnarray*}$$

The qubit over which the gate is acting is underlined, and is the one that 'couples' the computational states. By using this scheme, all gates be computed in parallel, regardless of the number of qubits worked on. The program that realizes these operations is publicly available in [25].

Let us compare the speed of the two setups, namely GPU and a usual CPU implementation. For that, we evaluate the average speed for the application of one time step of the KI model. That is, the application of the unitary operation (4) to a random state. We apply the operator several times for smaller system sizes, to obtain good statistics. The results are presented in figure A1 and show a more than satisfactory speed increase.

Zoom In Zoom Out Reset image size