Two dimensional kicked quantum Ising model: dynamical phase transitions

Using an efficient one and two qubit gate simulator, operating on graphical processing units, we investigate ergodic properties of a quantum Ising spin 1/2 model on a two dimensional lattice, which is periodically driven by a $\delta$-pulsed transverse magnetic field. We consider three different dynamical properties: (i) level density and (ii) level spacing distribution of the Floquet quasienergy spectrum, as well as (iii) time-averaged autocorrelation function of components of the magnetization. Varying the parameters of the model, we found transitions between ordered (non ergodic) and quantum chaotic (ergodic) phases, but the transitions between flat and non-flat spectral density {\em do not} correspond to transitions between ergodic and non-ergodic local observables. Even more surprisingly, we found nice agreement of level spacing distribution with the Wigner surmise of random matrix theory for almost all values of parameters except where the model is essentially noninteracting, even in the regions where local observables are not ergodic or where spectral density is non-flat. These findings put in question the versatility of the interpretation of level spacing distribution in many-body systems and stress the importance of the concept of locality.


Introduction
Quantum dynamics of strongly interacting quantum systems is one of the most exciting fields of current physics research, in particular due to the fact that 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 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 more mathematical physics side, the key question is a precise understanding of the notion of quantum ergodicity, and the ergodic to non-ergodic transitions in 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 and perhaps more mysterious from the one studied in many-body systems with on-site disorder [5,6,7,8].
Some conceptual and computational attempts in this direction have been made in Refs. [9,10,11,12,13,14,15]. In particular, one of us showed [9,10,11,12] that kicked XX − Z chains and kicked Ising spin chains typically exhibit nonergodic to ergodic transition with increasing period of the kicking, 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 nonintegrable regimes. The fact that in the semi-classical limit of 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 lead to an intuitive belief that 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. Yet a third measure of quantum chaos or quantum ergodicity in periodically driven quantum systems can be introduced, going under different names as Loschmidt echo, survival amplitude, or equivalently, (Fourier transform of) the quasienergy spectral density. We shell simply 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 opens of what happens in higher dimensional spin lattices.
In this paper we study probably the simplest non-trivial dynamical many-body model in two dimensions, namely the two-dimensional version of kicked Ising spin 1/2 model introduced in Ref. [12]. We numerically accurately calculate different dynamical and spectral properties of the model on a finite periodic rectangular L x × L y lattices implementing local one and two qubit operations on graphical processing units (GPUs) Two dimensional kicked quantum Ising model: dynamical phase transitions 3 of a desktop computer. Namely we compute phase diagrams (dependencies on the model's parameters) of dynamical susceptibilities, spectral densities and quasienergy level spacing distributions. We find that in most parts of parameter space these quantities turn out to be stable against increasing the lattice sizes, so we formulate certain conjectures about the thermodynamic behavior. We find, firstly, similarly as in one-dimensional chains, well defined transitions between non-ergodic and ergodic dynamical susceptibilities of local observables. Secondly, we find non-trivial transitions between flat spectral density and spectral densities with the shapes essentially determined by a single Fourier mode, ρ(φ) = 1 2π + c cos(kφ), for some constant c and integer k. The Fourier coefficient c seems to be decreasing with increasing the lattice size, but for any fixed Hilbert space dimension one finds a dramatic difference of |ρ(φ) − 1 2π | from a prediction of a flat spectrum (say of a random unitary matrix). However, remarkably and surprisingly the points of transitions of dynamical susceptibilities do not correspond to the points of transitions of level density. Thirdly, 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 that there are no nontrivial 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 which can put under 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 values of the parameters of the model. A first picture of the model emerges. In the next part (section 4), we discuss the correlation properties of the spectra, using the nearest neighbour spacing distribution. The second picture emerges. We next (section 5) 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.

The two dimensional quantum kicked Ising model
We study a periodic 2 dimensional lattice of kicked spin 1/2 particles, inspired by the Kicked Ising (KI) chain, proposed by one of the authors in [12]. The particle at site m ∈ {1, . . . , L x }, n ∈ {1, . . . , L y } will be described by standard Pauli operators σ α m,n , α ∈ {x, y, z}.
Let us start by defining a 2 dimensional Ising Hamiltonian with periodic boundary conditions σ α m,Ly ≡ σ α m,0 , σ α Lx,n ≡ σ α 0,n . We now define a b · σ m,n = b · S, S =: Notice that we can always choose the coordinate system such that b = (b x , 0, b z ), so that both H 0 and H 1 are real. We will normally consider only a transverse field, that is, b = (b x , 0, 0). The parameters, J (inter-spin 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 δ−pulses of period τ One-step quantum evolution propagator for the KI model over one period of the model -the so-called Floquet map -reads, setting τ = 1 by a free choice of units: where

Symmetries
We shall now briefly discuss some obvious symmetries of the model which help in reducing computational complexity of simulations.
Parameter space symmetries. Notice that the system is periodic in the parameters since the spectra of operators H I andb · S, whereb = b/| b| is the unit vector in the direction of b, form subsets of integers. In fact, in the spectrum of H I there can be only integers with fixed remainder of division with 4, since flipping an arbitrary spin can change H I only by ±4 or ±8. More precisely, the spectrum of H I consists of points Similarly, the spectrum ofb · S consists of points Let us now further assume that the field is transverseb = (1, 0, 0) as will be the case for most of this paper. Then, performing a checkerboard canonical (unitary) transformationĈ, namely: flipping the signs of y, z components of spins σ m,n for all even m + n (i.e. rotating for angle π along the x-axis), one finds that Two dimensional kicked quantum Ising model: dynamical phase transitions   5 Similarly, canonical transformationD, which flips x, y components of all spins (rotates around z axis for angle π) yields Therefore, changing the sign of J or b x leaves invariant all physical properties of the model, in particular the spectrum of U KI , so the principal domain of the phase diagram of the transverse field KI model only consists of a rectangle ( 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 vertical axis R y , reflexion over the horizontal axis R x , vertical translation T y , and horizontal translation T x . To illustrate these symmetries, let us number the sites in a (L x = 4) × (L y = 3) grid from left to right, and bottom to top, and consider a state of the computational basis |ψ 0 = |i 0 i 1 · · · i 11 , with i j ∈ {0, 1}. The action of R y on the grid will be to transform it into so The action of the horizontal reflection R x is similar: and the action over a member of the computational basis is Similarly we can picture the effects of the translations. The effect of vertical and horizontal translations on the original grid are respectively. Thus, the action of the symmetry is simply and It can also be noted that T y = T Lx , where T is the translation operator that acts as T | ψ 0 = |i 11 i 1 i 2 · · · i 10 . The symmetry subspaces of the Hilbert space are therefore specified by two quasi-momenta k x ∈ {0, . . . , L x − 1}, k y ∈ {0, . . . , L y − 1}, and for symmetric sectors with k x,y = 0 or k x,y = L x,y /2 by additional reflection signs π x , π y ∈ {±1}, such that the states |ψ from the subspace satisfy R x,y |ψ = π x,y |ψ , T x,y |ψ = e −2πikx,y/Lx,y |ψ .
When we consider the special case of a transverse magnetic field, another symmetry arises. A parity operator m,n σ x m,n commutes with the Floquet operator: [ m,n σ x m,n , U KI ] = 0. Our KI model also has an anti-unitary symmetry namely if K is a complex conjugation in the standard Pauli basis, then Using the standard wisdom [16], the model should then -if 'quantum chaotic'correspond to Circular orthogonal ensemble (COE) of random unitary symmetric matrices.

Steady field limit
With the same machinery we can study the time independent limit, corresponding to keeping b/J fixed, while letting |J| go to zero. Then, the scaled Hamiltonian is simply In order to study this Hamiltonian we also used the 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 that the changes in the figures presented are so small that cannot be noticed.

Spectral density
The spectrum of the Floquet map, S = {φ n ; n = 1, . . . , N := 2 LxLy }, defined by the unitary eigenvalue problem entails the main dynamical features of the model. The statistical properties of 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 {φ n } for many-body quantum models, despite the fact that full many-body quantum dynamics is becoming experimentally accessible in recent years, in particular in cold atom laboratories [17]. Even the behaviour of the simplest spectral characteristic, the 1-point function or the spectral density, defined as would be of great interest to know. In autonomous (time-independent) quantum manybody systems the spectral density is predicted to go to a Gaussian in thermodynamic limit [18,19], while for periodically driven quantum systems one may perhaps intuitively expect (and observe, in 1D chains [11]) that the Floquet quasienergy spectral density would be the constant (flat) function ρ(φ) = 1 2π , in a generic case. The spectral density ρ(φ) is a 2π-periodic function and therefore can be represented in terms of the Fourier modes as where the Fourier coefficients ρ k are given as traces of the k−step KI propagator The symmetry property ρ −k = ρ k has been used, following from the symmetry of the spectra of H 0 and H 1 around zero energy and cyclicity of the trace. The sum in (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 whole spectrum for sizes of up to 5 × 4. The behavior over different symmetry sectors seems to be similar, in all examples that we considered, see Fig. 1 (right panel).
For efficient numerical computation of the leading Fourier components ρ 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 where dµ(ψ) is the measure induced by the Haar measure over the unitary group, and |ψ random 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 |ψ 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 behaviour of the spectral density, for several parameter values, and several sizes. The left panel of fig. 1 suggests a dominant Fourier component ρ 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 proofs very useful. Such analysis is carried out for all coefficients up to k = 12, varying the transverse component of the magnetic field. A similar behaviour is obtained for different sizes as can be appreciated 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 |ρ noise |, and another region, in which there is an ordered phase, manifested by |ρ k | |ρ noise |. We estimate |ρ noise |, replacing  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 |ρ 3 (b = 0.2)/ρ noise | ∝ ∼ exp(0.27L) thus sharpening the transition from a disordered to an ordered phase. Therefore, even though the spectral density ρ(φ) seems to universally approach a constant 1/(2π) when one approaches the thermodynamic limit L x , L y → ∞, there are discontinuous transitions on the size scaling of the deviation ρ(φ) − 1 2π with changing the system parameters.
One can get 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 ρ k is dominating, as a function of model's parameters. We shall consider the spectrum to be flat, k = 0, if all Fourier coefficients ρ k , for k > 0,  On the left panel, we observe the nearest neighbour spacing distribution for three different transverse fields, bx = 0.2, 0.3 and 0.5 in red, green and yellow respectively, J = 0.5, and we consider a 5 × 4 lattice. In all cases, we are considering sx = ±1, kx ∈ {1, 2} and ky = 1. The thick black curve correspond to the Wigner surmise. In the inset, we show the average of these three curves, minus the Wigner surmise, together with the theoretical prediction. On the right panel, we consider the Kolmogorov distance between the unfolded P (s), and the Wigner surmise, for all the parameters of the model, and a 4 × 3 lattice. Very good agreement with the RMT prediction is observed except when J or bx are zero, or J = bx = π/4. are comparable to ρ noise . That is, the spectrum is declared flat, if On the other hand, we consider the system to be in the phase k > 0, if |ρ k | > |ρ k |, for all k = k = 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 which seems remarkable stable.

Level spacing distribution
In order to analyze the correlation properties of the spectrum we used the commonly studied nearest neighbour spacing distribution. That is, we consider the distribution of level spacingss where φ i are the sorted eigenphases of the evolution operator. In order to remove the effect of non-uniform level density we perform unfolding, i.e., a smooth non-linear scaling of the eigenvalues in order to get a uniform spectral density. Since the density is typically well described by few Fourier components, we shall take 6 of them to numerically perform the unfolding. Thus, we shall use the mapping where α k = ρ k /π is determined numerically using direct evolution. In this way we can unfold the spectra to obtain fairly flat distributions of unfolded level spacings.
Two dimensional kicked quantum Ising model: dynamical phase transitions

11
Another very important aspect that must be taken into account is the fact that different symmetry sectors are not statistically correlated with each other, so we must rather look at the distribution of where K = {k x , k y , π x , π y , Π} denotes the set of quantum numbers that determines the irreducible quantum sector. By construction, the mean spacing s K i equals one, so the probability density P (s) of {s K i } is normalized such that ∞ 0 P (s)ds = ∞ 0 sP (s)ds = 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., hyperbolic dynamical system). Given the time reversal symmetry, the corresponding ensemble, in our case, would be the circular orthogonal ensemble (COE). To a very good approximation, the level spacing distribution is given in terms of 2 × 2 random real symmetric matrices, the so-called Wigner surmise P Wigner (s) = π 2 exp(−πs 2 /4). A similar conjecture has been suggested for strongly non-integrable quantum many-body systems [20], but it has not been established precisely (yet), 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 completely 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 excellently described by COE or Wigner surmise. In 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, X |f (x) − g(x)|dx, between the observed distribution of nearest neighbour spacings and the 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 π/2), or specially commensurate fields where the spectrum of U KI can be explicitly computed in terms of regular or number-theoretic functions.

Dynamical susceptibilities and non-ergodicity to ergodicity transition
So far we have analyzed dynamical properties of the 2D kicked Ising model which depend solely on its spectrum. Now we shall focus on dynamical correlations of local observables, which are the key input to any linear response treatment of condensed matter theory [21,22,23].
Consider a traceless observable M (typically extensive and local), say a component of magnetization M = S ν , ν ∈ {x, y, z}. We define its time-autocorrelation with respect to the kicked Ising dynamics as If tr M = 0, the corresponding constant has to be subtracted from C(t). Here M (t) denotes the time-dependent observable in the Heisenberg picture: One measures the ergodicity of an observable M by the so-called dynamical susceptibility, defined as the time-average of C(t): By definition, observable M is ergodic with respect to dynamics U KI (t), if χ M = 0, and non-ergodic otherwise. Note that χ M is always nonnegative as it represents the spectral weight, i.e. the power spectrum of M (t), at frequency ω = 0. For numerical investigations, in order 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 ∈ Z, as In order to illustrate the ergodic properties of the model for different parameters, we have analysed 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 =n · S withn an arbitrary unit vector. Small system sizes revealed that whenb 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 = π/4 in parameter space. However, this is the case for general observables, although there is a strong tendency to be exactly symmetric. More general observables also have this tendency, that becomes increasingly more difficult to explore for moderate large systems. We found three qualitatively different kinds of behaviour of dynamical susceptibility, exemplified in figure 5. Here we compare the behaviour 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 neither for large times nor dimensions. There seems to be a nonvanishing asymptotic value χ M = 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 χ M = 0. Finally, there are parameter values (say b x = 0.5) for which the asymptotic value, is reached exponentially fast and again the asymptotic value decreases with the Hilbert space dimension, suggesting that it vanishes in the thermodynamic limit.
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 χ M = 0. However, for J = 0.5 there seems to be a range of b x where the correlations clearly vanish, namely for 0.4 < b x < 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 phase diagram (of order parameter χ M ) in the parameter space (J, b x , b z ), clearly indicating distinct regions of ergodic and non-ergodic dynamics.
Comparing these data to phase diagrams of level density figure 1, or nearest neighbour 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 ergodic or non-ergodic regime for the dynamics of local observables. This speculative conclusions 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 to have a glimpse into 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.

Conclusions
In this paper we describe a computational excursion into ergodic properties of twodimensional 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 nontrivial spectral densities dominated by nonzero Fourier components. Local observables display ergodic regimes with decaying correlations and non-ergodic regimes with non-decaying correlations, which however, 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 nonintegrable, and therefore, surprisingly, does not provide any useful information on system's ergodicity. We believe that our numerical results generate a strong motivation for further theoretical investigations into dynamics of periodically driven (or discretetime) interacting spin models on 2D lattices. TP acknowledges financial support by the grant P1-0044 and J1-5439 of the Slovenian Research Agency. Support by the projects CONACyT 153190 and UNAM-PAPIIT IA101713 is acknowledged by CP and EV. computational base we specify the entries on the state a specific thread will compute on, for example when applying a 1-qubit gate on the second qubit, each thread shall act on the coefficients of the components. In particular, thread 0: | · · · 0000 , | · · · 0100 thread 1: | · · · 0001 , | · · · 0101 thread 2: | · · · 0010 , | · · · 0110 . . . thread 2 L−1 : | · · · 1011 , | · · · 1111 .
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 con be computed in parallel regardless of the number of qubits it works on. The program that realizes these operations is publicly available in [24]. 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 Kicked Ising model. That is, the application of the unitary operation (4) to a random state. We apply the operator several times for smaller system sizes, so as to get good statistics. The results are presented in figure A1, and show a more than satisfactory speed increase. . We can see that around 13 qubits already we have an important speed factor increase, which stabilizes around 240x. We are limited by the size of the memory of the card to 25 qubits.