Slow scrambling and hidden integrability in a random rotor model

TL;DR

This paper investigates the dynamics of quantum rotors with random interactions, revealing slow scrambling without exponential growth and uncovering integrability and connections to random-matrix theory.

Contribution

It demonstrates that in a large N, M limit, the model is integrable and exhibits slow scrambling, with a small exponential growth due to irrelevant terms.

Findings

Squared commutators do not grow exponentially in the large N, M limit.

The model is integrable and connected to random-matrix theory.

Irrelevant terms cause a small exponential growth in commutator dynamics.

Abstract

We analyze the out-of-time-order correlation functions of a solvable model of a large number, $N$, of $M$-component quantum rotors coupled by Gaussian-distributed random, infinite-range exchange interactions. We focus on the growth of commutators of operators at a temperature $T$ above the zero temperature quantum critical point separating the spin-glass and paramagnetic phases. In the large $N,~M$ limit, the squared commutators of the rotor fields do not display any exponential growth of commutators, in spite of the absence of any sharp quasiparticle-like excitations in the disorder-averaged theory. We show that in this limit, the problem is integrable and point out interesting connections to random-matrix theory. At leading order in $1/M$, there are no modifications to the critical behavior but an irrelevant term in the fixed-point action leads to a small exponential growth of the squared commutator. We also introduce and comment on a generalized model involving $p$-pair rotor interactions.

Equations

$$H=g\sum _i\frac{L_i^2}{2M}+\frac{M}{\sqrt[]{N}}\sum _{i<j}{J}_{ij}{\mathbf{n}}_i\cdot {\mathbf{n}}_j,$$

$$Q(i{\omega }_n)=\frac{g}{{\omega }_n^2+\lambda -g{J}^{2}Q(i{\omega }_n)},$$

$$A(\omega )=\text{sgn}(\omega )\frac{\sqrt{({\omega }^2-{\Delta }_s^2)({\Delta }_s^2+4Jg-{\omega }^2)}}{2g{J}^2}$$

$$Q_R(\omega )\approx -\frac{{\omega }^2}{2g{J}^2}+\frac{1}{J}+i\frac{\text{sgn}(\omega )}{J\sqrt{gJ}}|\omega |.$$

$${\Delta }_s(T)=\frac{2\pi k_{B}T}{\sqrt{3}{\log }^{1/2}(1/T)}.$$

$$C(t_1,t_2)=\text{Tr}\{\rho [V(t_1),W(0){]}^{†}[V(t_2),W(0)]\},$$

$$F(t_1,t_2)=-\sum _{i,j,a,b}\overline{\text{Tr}\{\sqrt{\rho }[n_i^a(t_1),n_j^b(0)]\sqrt{\rho }[n_i^a(t_2),n_j^b(0)]}\},$$

$$F({\omega }_1,{\omega }_2)\equiv \int \int F(t_1,t_2)e^{i({\omega }_1{t}_1+{\omega }_2{t}_2)}d{t}_{1}d{t}_2.$$

$$F({\omega }_1,{\omega }_2)=\frac{N}{M}\frac{Q_R({\omega }_1)Q_R({\omega }_2)}{1-J^2{Q}_R({\omega }_1)Q_R({\omega }_2)}.$$

$$Q(i{\omega }_n)$$

$${\Sigma }_R(\omega )\sim u^2|\omega {|}^5\text{sgn}(\omega ),$$

$$F(t_1,t_2)$$

$$Q_W(t){\delta }_{ij}{\delta }_{ab}$$

$$Q_W(\omega )=\frac{A(\omega )}{2\sinh (\beta \omega /2)}\approx \{\begin{array}{l}\frac{\sqrt{gJ}}{g{J}^2}\frac{\text{sgn}(\omega )\sqrt{{\omega }^2-{\Delta }_s^2}}{2\sinh (\beta \omega /2)},|\omega |>{\Delta }_s\\ 0,\text{otherwise}.\end{array}$$

$$\displaystyle\mathcal{S}=-\int d\tau\sum_{p=1}^{N}\frac{M}{2}\bigg{[}\frac{1}{g}[(\partial_{\tau}\boldsymbol{V}_{p})^{2}+\lambda(\boldsymbol{V}_{p})^{2}]+\frac{1}{\sqrt{N}}\lambda_{p}(\boldsymbol{V}_{p})^{2}+...\bigg{]},$$

$${\tilde{Q}}_p(i{\omega }_n)=\frac{1}{({\omega }_n^2+\lambda )/g+{\lambda }_p/\sqrt{N}}.$$

$${\sigma }_J(x)=\frac{1}{2\pi J}\sqrt{4N-(x/J{)}^2},$$

$$Q(i{\omega }_n)$$

$${\mathbf{n}}_i(\tau )=\sum _{p=1}^N{\phi }_i^{(p)}{\mathbf{V}}_p(\tau ).$$

$$H=g\sum _i\frac{{\widehat{L_i}}^2}{2M}+\frac{M^{\frac{p+1}{2}}}{\sqrt{p}{N}^{p-1/2}}\sum _{\{i,j\}}{J}_{i_1{j}_1...i_p{j}_p}\prod _{k=1}^p({\mathbf{n}}_{i_k}\cdot {\mathbf{n}}_{j_k}),$$

$$Q(i{\omega }_n)$$

$$\Sigma (\tau )$$

$$Q(\tau )$$

$$\Sigma (i\omega )$$

$$[L_{i\mu \nu },n_{j\sigma }]=i{\delta }_{ij}({\delta }_{\mu \sigma }{n}_{j\nu }-{\delta }_{\nu \sigma }{n}_{j\mu }).$$

$$Z_n$$

$$\displaystyle Z=\int\mathcal{D}n^{a}(\tau)~{}\mathcal{D}\lambda~{}\mathcal{D}Q^{ab}~{}\textnormal{exp}\bigg{[}-\frac{M}{2g}\int d\tau[(\partial_{\tau}n^{a})^{2}+\lambda((n^{a})^{2}-1)]$$

$$\displaystyle+\frac{MJ^{2}}{2}\int d\tau\int d\tau^{\prime}[Q^{ab}(\tau-\tau^{\prime})n^{a}(\tau)\cdot n^{b}(\tau^{\prime})-\frac{1}{2}Q^{ab}(\tau-\tau^{\prime})^{2}]\bigg{]}.$$

$$Q(i{\omega }_n)=\frac{{\omega }_n^2}{2g{J}^2}+\frac{\lambda }{2g{J}^2}-\frac{1}{2g{J}^2}\sqrt{({\omega }_n^2+\lambda -2gJ)({\omega }_n^2+\lambda +2gJ)}.$$

$$\displaystyle\int_{0}^{\infty}\frac{d\omega}{\pi}A(\omega)\coth\bigg{(}\frac{\beta\omega}{2}\bigg{)}=1.$$

$$gJ=\frac{9{\pi }^2{J}^2}{16}-3T^2.$$

$$F(\omega_{1},\omega_{2})=\frac{N}{M}\frac{Q_{R}(\omega_{1})Q_{R}(\omega_{2})}{1-J^{2}Q_{R}(\omega_{1})Q_{R}(\omega_{2})}=\frac{N}{MJ^{2}}\bigg{(}-1+\frac{1}{1-J^{2}Q_{R}(\omega_{1})Q_{R}(\omega_{2})}\bigg{)}.$$

Full text

Introduction.- The route to thermalization Deutsch (1991); Srednicki (1994); Tasaki (1998); Rigol et al. (2008), or lack thereof Anderson (1958); Fleishman and Anderson (1980), in closed quantum systems remains one of the central open questions in quantum statistical mechanics. The problem has been studied extensively over the past few decades and has received renewed attention in recent years Nandkishore and Huse (2015); Altman and Vosk (2015), partly as a result of great advances in experimental techniques that allow us to probe quantum dynamics in closed systems Kaufman et al. (2016); Langen et al. (2013); Jurcevic et al. (2014); Schreiber et al. (2015); Choi et al. (2016). The eventual global thermalization of closed quantum systems may, in general, proceed through multiple distinct stages. At the shortest timescales, there is a decay of local perturbations through the process of relaxation and this timescale affects the time dependence of simple two-point correlation functions of local observables Forster (1975). At intermediate timescales, a characteristically distinct phenomenon of scrambling Sekino and Susskind (2008) leads to the spreading of information across all of the degrees of freedom of the system. The connection between scrambling and thermalization remains a topic of active research.

It has been suggested that certain non-trivial out of time ordered correlation (OTOC) functions Larkin and Ovchinnikov (1969) of local operators (defined in Eq. 6 below) can probe the onset of scrambling Shenker and Stanford (2014a); Kitaev (2014). In recent years, these correlation functions have been used extensively to study scrambling in black holes Shenker and Stanford (2014a); Kitaev (2014); Shenker and Stanford (2014b), which are supposed to be the fastest known scramblers Maldacena et al. (2016) in nature. Moreover, since the memory of the initial state is effectively lost Hosur et al. (2016) due to scrambling, it can be heuristically related to an onset of ‘many-body quantum chaos’ Kitaev (2014) with a positive ‘Lyapunov’ exponent. A number of calculations in a variety of models at large $N$ or at weak-coupling have found an exponential growth of the OTOC (see e.g. Maldacena and Stanford (2016); Stanford (2016); Swingle and Chowdhury (2017); Patel and Sachdev (2017); Chowdhury and Swingle (2017); Aleiner et al. (2016); Gu et al. (2017); Banerjee and Altman (2017)) but a clear signature of such growth in numerical studies has been absent Bohrdt et al. (2017); Luitz and Bar Lev (2017); Heyl et al. (2018). Certain random-circuit models also explicitly point out the absence of such exponential growth of OTOC Nahum et al. (2018); von Keyserlingk et al. (2018) but it may be possible to reconcile these differences Xu and Swingle (2018). The precise relationship between scrambling and chaos remains unclear.

Inspired by the rapid developments in the study of OTOC in many-body systems, including in the family of Sachdev-Ye-Kitaev (SYK) models Sachdev and Ye (1992); Kitaev (2015); Kitaev and Suh (2018) which do not admit a description in terms of weakly interacting quasiparticles and scramble at the maximally allowed rate Maldacena and Stanford (2016), we address a number of questions in the remainder of this paper: (i) Are there $(0+1)-$dimensional models different from the SYK model that also scramble at the maximal rate? (ii) Does the absence of well-defined quasiparticle excitations ensure (near-)maximal chaos? (iii) In models with quenched disorder, is it possible to distinguish between the phenomenon of ‘dephasing’ arising as a result of inelastic processes due to strong interactions as opposed to coupling to disorder?

In this letter we study scrambling in a class of rotor models Ye et al. (1993); Read et al. (1995) that share some superficial resemblance to the SYK models but have remarkably distinct properties. We will begin with a large number, $N$, of $M$ component real rotor fields with $O(M)$ symmetric all-to-all interactions that are assumed to be quenched random variables with zero mean Ye et al. (1993). The system can be tuned through a quantum critical point (QCP) separating a paramagnet and a quantum spin glass. Upon disorder averaging, it is known that there are well-defined quasiparticle excitations on either side of the transition but the QCP itself lacks any quasiparticle-like excitations Ye et al. (1993). At low energies, the temperature ($T$) serves as the only relevant scale at the QCP Sachdev (2011), and one may then naturally expect the QC regime to scramble with a rate ($\lambda_{L}$) that is proportional to $(k_{B}T/\hbar)$ Patel and Sachdev (2017); Chowdhury and Swingle (2017).

Instead we find results that are surprisingly at odds with this natural expectation. In particular: (i) in the $N,~{}M\rightarrow\infty$ limit, the OTOC does not exhibit any exponential growth (i.e. $\lambda_{L}=0$) in the QC regime, (ii) the model can be solved (in the $N\rightarrow\infty$ limit) to all orders in $1/M$ Ye et al. (1993); Read et al. (1995); the leading (irrelevant) correction in $1/M$ leads to a small $\lambda_{L}>0$, (iii) the absence of an exponential growth of the OTOC in the former limit, in spite of the apparent non-quasiparticle character, is related to an underlying ‘hidden’ integrability. We demonstrate this by pointing out an interesting connection to a ‘random-matrix’ description Mehta (2004) for the saddle point of the $O(M)$ rotor model at large $N,~{}M$. Finally, inspired by these new insights, we construct a generalized rotor model with random exchange interactions between $p~{}(>2)$-pairs of rotors and point out the associated similarities with the generalized SYKq models with $q~{}(>4)$-fermion interactions Gross and Rosenhaus (2017).

Random $O(M)$ rotor model.- We consider a model of $O(M)$ rotors defined on $N$ sites with random interactions between any two sites. There is no notion of “space” in this model, which is effectively $(0+1)-$dimensional. The Hamiltonian of the model is given by,

[TABLE]

where $i,j\in\{1,...,N\}$ are site indices, ${\mathbf{n}}_{i}$ are $M$ component vectors of unit length, and $L_{i}$ represent angular momentum tensors with $M(M-1)/2$ components. The components of ${\mathbf{n}}_{i}$ are mutually commuting and satisfy standard commutation relations si . The $J_{ij}$ are uncorrelated random variables selected from a Gaussian distribution, $P(J_{ij})=e^{-J_{ij}^{2}/2J^{2}}/\sqrt{2\pi J^{2}}$.

We can obtain the saddle point action in the $N\rightarrow\infty,~{}M\rightarrow\infty$ limit, by integrating over disorder configurations Ye et al. (1993). In particular, the self-consistency condition for $Q(i\omega_{n})$, the Fourier transform of the correlator $Q(\tau)=\langle{\mathbf{n}}(\tau)\cdot{\mathbf{n}}(0)\rangle$ (assuming a replica-symmetric solution) is given by si

[TABLE]

where $\omega_{n}$ is a Bosonic Matsubara frequency and $\lambda$ is a Lagrangian multiplier that imposes the constraint of unit length for the rotor fields.

It is immediately clear that the spectral function, $A(\omega)=\textnormal{Im}~{}Q(i\omega_{n}\rightarrow\omega+i0^{+})$, is given by,

[TABLE]

where we have defined $\Delta_{s}^{2}=(\lambda-2gJ)$, such that $A(\omega)$ is finite if $\Delta_{s}^{2}<\omega^{2}<(\Delta_{s}^{2}+4Jg)$ and is zero otherwise.

It is evident then that $\Delta_{s}$ is the spin-gap to excitations in a paramagnetic phase, which vanishes at the $T=0$ QCP (when $\lambda=2Jg$) between the paramagnet and a spin-glass phase. At the QCP, for frequencies $\omega^{2}\ll gJ$, $Q_{R}(\omega)\equiv Q(i\omega_{n}\rightarrow\omega+i0^{+})$ is given by

[TABLE]

The non-analytic piece in $A(\omega)\sim|\omega|$ (or equivalently, $Q(\tau)\sim 1/\tau^{2}$ for large $\tau$) encodes the non-trivial character of the quantum critical regime. Similar behavior was also obtained for the transverse-field Ising model with random exchange interactions Miller and Huse (1993), which is expected to be in the same universality class as the above model with $M=1$.

The solution of the constraint equation, equivalently expressed as $Q(\tau=0)=1$, can be used to obtain the self-consistently generated ‘thermal-mass’, $\Delta_{s}(T)$ in the QC region at $g=g_{c}=9\pi^{2}J/16$. This leads to Ye et al. (1993); si ,

[TABLE]

OTO Correlation functions.- For any two generic Hermitian operators, $V$ and $W$, we can define an “unregulated” squared commutator,

[TABLE]

where $V(t)=e^{iHt}Ve^{-iHt}$ is the time-evolved Heisenberg operator corresponding to the Hamiltonian $H$. The system is taken to be in thermal equilibrium at a temperature $T=\beta^{-1}$ with the density matrix $\rho\propto e^{-\beta H}$. For the purpose of our computations in this paper, we will instead choose to work with a “regulated” squared commutator Maldacena et al. (2016) defined in terms of the rotor fields (i.e. $V,~{}W\equiv{\mathbf{n}}$),

[TABLE]

where $i,j=1,..,N$ are site indices, $a,b=1,..,M$ are the vector indices and “${}^{\overline{~{}~{}~{}}}$” represents disorder averaging. For systems with a large number of degrees of freedom, $N,~{}M$, the expectation 111Recent work Liao and Galitski (2018) has pointed out that the regulated vs. unregulated OTOC can display $\lambda_{L}$ that are significantly different from one another; the bound Maldacena et al. (2016) was proposed for the regulated OTOC. in a non-integrable chaotic system is that $F(t_{1},t_{2})\sim\epsilon~{}e^{\lambda_{L}(t_{1}+t_{2})/2}$. Here $\epsilon$ is some small parameter, which may depend on time and on the number of relevant degrees of freedom.

For the system described by the Hamiltonian in Eq. 1, the leading contribution to $F(t_{1},t_{2})$ arises from the series of diagrams shown in Fig.1(a) in the $N,~{}M\rightarrow\infty$ limit. The two horizontal lines represent the retarded propagators of the ${\mathbf{n}}$ fields along two real time branches respectively and the vertical dotted lines represent the disorder contraction over the random couplings, $J_{ij}$. Note that the ladder series is not simply a perturbative expansion in $J$, since the ${\mathbf{n}}$ propagators also include a self-consistent renormalization from the disorder averaging.

The ladder series in Fig.1(a) can be evaluated in a straightforward fashion and is given by a geometric series when expressed in terms of,

[TABLE]

The sum evaluates to,

[TABLE]

The exponential growth in time, if any, of the squared commutators in this case will arise from the (imaginary) poles of $F(\omega_{1},\omega_{2})$ in Eq. 9. However a simple exercise immediately shows that there are no such poles of $F(\omega_{1},\omega_{2})$ and as a result $F(t_{1},t_{2})$ does not exhibit any exponential growth with $\lambda_{L}>0$.

Corrections beyond large $M$.- In order to obtain the first non-trivial contribution to $\lambda_{L}$, we have to consider corrections beyond the $M\rightarrow\infty$ limit considered thus far. Previous work Ye et al. (1993); Read et al. (1995) has shown that the $1/M$ corrections in the disorder-averaged theory can be studied systematically within a soft-spin formulation of the problem. The leading $1/M$ correction arises from the four-body interaction term $\sim u({\mathbf{n}}^{2})^{2}$, which modifies the ladder series for $F(t_{1},t_{2})$ (see Fig.1(b)) and introduces a correction to the self-energy for $Q_{R}(\omega)$ (Fig.2). As has been noted earlier Stanford (2016); Swingle and Chowdhury (2017); Patel and Sachdev (2017); Chowdhury and Swingle (2017), the former can contribute to an exponential growth and the latter to an exponential decay for $F(t_{1},t_{2})$, which may lead to a net exponential growth with $\lambda_{L}>0$.

Including the $1/M$ correction, the propagator has the form

[TABLE]

which can be solved in a straightforward fashion and where the leading non-analytic correction to the self-energy (Fig.2), after analytically continuing to real frequencies si gives

[TABLE]

which is less singular than the $M\rightarrow\infty$ saddle point correction and hence does not affect the long time behavior of $Q(t)$. However, the presence of this self-energy correction immediately leads to an exponential decay for $F(t_{1},t_{2})\sim e^{-\Gamma(t_{1}+t_{2})}$ with a small decay rate $\Gamma=\alpha(u^{2}/M)\Delta_{s}^{4}(T)$, where $\alpha$ is a known function of $g,J$ si .

On the other hand, the ladder series for $F(t_{1},t_{2})$ itself is modified due to the addition of a new ‘rung’ (Fig.1 (b)), such that the self-consistent Bethe-Salpeter equation is given by

[TABLE]

where

[TABLE]

is the Wightman Green’s function. Note that the first line in the expression for $F(t_{1},t_{2})$ represents the contribution that we already evaluated earlier while the second line represents the contribution from the new rung. The exponential growth of $F(t_{1},t_{2})$, if any, will arise from the eigenvalues for the second line above, i.e. $\lambda_{L}\sim O(u^{2}/M)$. At weak coupling (i.e. small $u$), we are ignoring here another contribution to the ladder series at $O(u^{4}/M)$. In the $\omega\ll\sqrt{gJ}$ limit,

[TABLE]

(See Appendix C of Ref. Chowdhury and Swingle (2017)). An explicit numerical evaluation of the contribution to the small exponentially growing piece leads to $F(t_{1},t_{2})\sim e^{\Upsilon(t_{1}+t_{2})}$ with $\Upsilon=\eta(u^{2}/M)T^{3}\Delta_{s}(T)$, where $\eta$ is a known function of $g,~{}J$ si . Thus, at low temperatures $\Upsilon$ is parametrically larger than $\Gamma$ with $\lambda_{L}=(\Upsilon-\Gamma)>0$, resulting in an overall exponential growth of $F(t_{1},t_{2})$. Thus at low temperatures in the quantum critical regime, we find $\lambda_{L}\sim(u^{2}/M)[T^{4}/\log^{1/2}(1/T)]$ as shown in Fig. 3.

Integrability and random matrix model.- Let us now consider a collection of $M$ component vector fields, $\boldsymbol{V}_{p}$, with a unit constraint $\boldsymbol{V}_{p}^{2}=1$, where $p=1,..,N$. Let the ‘frequencies’ (or ‘mass’) terms for these fields be drawn from the distribution of eigenvalues of a large $N\times N$ random matrix $J_{\ell m}$, whose elements are independently distributed, such that $\sum_{m}J_{\ell m}\varphi_{m}^{(p)}=\lambda_{p}\varphi_{\ell}^{(p)}$. The $\lambda_{p}$ are the eigenvalues and $\varphi^{(p)}$ denote the eigenvectors of the random matrix $J_{\ell m}$. For any given realization of $J_{\ell m}$, the Lagrangian can be written as,

[TABLE]

where $\lambda$ is a Lagrangian multiplier that imposes the unit constraint and the action uses a specific normalization with respect to $M,~{}N$ in order to relate to the saddle-point action of the random rotor model. The ellipses denote additional terms beyond quadratic order that we leave unspecified for now.

The auto-correlator for the field $\boldsymbol{V}_{p}$ in any given realization is given by,

[TABLE]

It is then natural to ask how the disorder averaged auto-correlator behaves as a function of frequency, after averaging over the distribution of different $\lambda_{p}$. For the Gaussian orthogonal ensemble (GOE) of random matrices, the level density (or density of states), $\sigma_{J}(\lambda_{p})$, is given by,

[TABLE]

otherwise also known as the ‘Wigner semi-circle law’ Mehta (2004). The disorder averaged auto-correlator $Q(i\omega_{n})$ can be evaluated,

[TABLE]

It is not a coincidence that the above propagator is identical to the propagator for the rotor-fields, obtained upon solving Eqn. 2.

The mapping between the two seemingly unrelated models can be made precise by considering the following orthogonal transformation

[TABLE]

The problem is thus essentially free in the $N,~{}M\rightarrow\infty$ limit, described by a random-matrix, and is responsible for the absence of an exponential growth of $F(t_{1},t_{2})$. We can similarly express the quartic interaction between the rotor fields in terms of the fields introduced above and systematically compute the leading $1/M$ corrections to the self-energy si . Disorder averaging now requires a knowledge of the distributions of both the eigenvalues, $\lambda_{p}$ and the eigenvectors, $\varphi_{i}^{p}$, which are known to be independent for GOE random-matrices Anderson et al. (2010). Within the random-matrix formulation of the problem, the correlations among the different eigenvalues are also included in the $n-$level “cluster functions” Anderson et al. (2010). The replica-action for the rotor problem does not include these correlations and assumes that the eigenvalues are completely independent; the completely uncorrelated piece in the cluster function, which factorizes into a product of the $\sigma_{J}(\lambda_{i})$, reproduces the results obtained using the replica action si . The alternative formulation in terms of the random-matrix therefore allows us to clearly disentangle the features that arise from disorder averaging and interactions.

A generalized model.- Inspired by the structure of the saddle-point action of the $O(M)$ random rotor model, we propose here an extension of the model to instead include an interaction between $p$-pairs of rotors (instead of $p=1$)

[TABLE]

where $J_{i_{1}j_{1}...i_{p}j_{p}}$, symmetric under $i_{k}\leftrightarrow j_{k}$ and $(i_{k},j_{k})\leftrightarrow(i_{l},j_{l})$, are Gaussian random variables with mean zero and finite variance, $J^{2}$. After integrating over disorder and assuming a replica-symmetric solution, the saddle-point equations become,

[TABLE]

However the scaling behavior at the QCP between a paramagnet and a spin-glass, if it exists 222The above model while different from the quantum spherical $p-$spin glass model, bears some superficial resemblance at the level of the replica action. It has been argued Cugliandolo et al. (2001) that for the latter, there is no continuous QPT between the paramagnet and spin-glass for $p>2$, and the transition is first-order ($p=2$ reduces to the usual random rotor model Ye et al. (1993) considered earlier)., can be obtained using scaling (when $\Sigma(i\omega=0)=\lambda$). For $p>1$, this yields,

[TABLE]

where $\Delta(p)=1/(p+1)$.

At this point, the correspondence between the generalized rotor models with $p\geq 1$ and the SYKq models Gross and Rosenhaus (2017) should be clear. The QCP for the rotor model with $p=1$ bears striking similarity with SYK2, a random-matrix model, and explains the lack of any exponential growth of the OTOC. On the other hand, the QCP for the rotor models with $p>1$ have a scaling structure that is similar to the SYKq models. We leave a detailed analysis of our model for future work.

Discussion.- In this work, we have studied the OTOC for a rotor model interacting with random exchange interactions at the QCP between a paramagnet and a spin-glass. There are interesting connections between the large$-N,~{}M$ saddle point action and a purely random-matrix model supplemented with $O(M)$ invariant quartic interactions. It will be interesting to study the leading $1/N$ corrections to the saddle point solution using the properties of the asymptotic Tracy-Widom probability distribution Majumdar (2007) of the largest eigenvalues of a random matrix.

Acknowledgements.- DC is supported by a postdoctoral fellowship from the Gordon and Betty Moore Foundation, under the EPiQS initiative, Grant GBMF-4303 at MIT. TS is supported by a US Department of Energy grant DE-SC0008739, and in part by a Simons Investigator award from the Simons Foundation.

Note added.- While this manuscript was being finalized for submission, we became aware of related results Cheng and Swingle (2019).

I SUPPLEMENTARY MATERIAL for “Slow scrambling and hidden integrability in a random rotor model”

II Replica action for random $O(M)$ rotor model

In this subsection, we provide some additional details for the large $N,~{}M$ saddle point treatment for the rotor model. To begin, recall that the commutation relations for the rotor fields and the angular momenta are given by,

[TABLE]

In the $N\rightarrow\infty$ limit, the saddle point action can be obtained upon integrating over disorder configurations. The replicated partition function in Euclidean time is then given by,

[TABLE]

where $a,b\in\{1,...,n\}$ are replica indices, $i,j\in\{1,...,N\}$ are site indices, and $m,n\in\{1,...,M\}$ are vector indices. $Q_{mn}^{aa}$ and $P_{mn}^{ab}$ correspond to quadrupolar and spin glass order, respectively. Since we are only interested in the paramagnetic regimes in this work, both of the above quantities have zero expectation value.

Assuming the saddle point to be $O(M)$ invariant, the replica action can be simplified such that it describes the quantum mechanical action for a single rotor with multiple replica indices,

[TABLE]

In the large $M$ limit, we can define the imaginary time ($\tau$) auto-correlation functions $Q(\tau)=\langle{\mathbf{n}}(\tau)\cdot{\mathbf{n}}(0)\rangle$, where we are only interested in regimes where $Q(\tau)$ is replica diagonal. We can solve for $Q(i\omega_{n})$ (which leads to the above self-consistency condition Eq. 2) and obtain,

[TABLE]

At $T=0$, the spin-gap vanishes at the critical point between a paramagnet and a spin-glass phase when $\lambda=2Jg$ (see fig.4); the actual phase-boundary in the $g-T$ plane can be obtained by solving for the constraint equation $Q(\tau=0)=1$ with the above value of $\lambda$, which can be recast as,

[TABLE]

Near the $T=0$ quantum critical point, this is given by,

[TABLE]

III $F(t_{1},t_{2})$ in the $N,~{}M\rightarrow\infty$ limit

In this subsection, we obtain the time dependence of $F(t_{1},t_{2})$ in the large $N,~{}M$ limit, where there is no exponential growth. Recall that,

[TABLE]

The quantum critical region of interest to us here is when $\Delta_{s}^{2}\ll gJ$. Let us introduce the variables,

[TABLE]

such that the denominator can be simplified to,

[TABLE]

where we have assumed $x,y\ll 1$ in the first line above and ignored terms that are $O(xy^{2},~{}x^{2}y)$ and higher. By making this approximation, we only focus on the long time, long period behavior with a characteristic time scale $1/\sqrt{2gJ}$.

We can then approximate $F(x,y)$ as

[TABLE]

The leading time dependence of $F(t_{1},t_{2})$ comes from the second term above, which can be simplified as,

[TABLE]

Rescaling frequencies and time as $\tilde{\omega}_{1,2}=\omega_{1,2}/\Delta_{s},~{}\tilde{t}_{1,2}=\Delta_{s}t_{1,2}$, the inverse Fourier transform of Eqn. 36 is given by

[TABLE]

We can redefine $\tilde{t}=\frac{\tilde{t}_{1}-\tilde{t}_{2}}{2}$ and $\tilde{T}=\frac{\tilde{t}_{1}+\tilde{t}_{2}}{2}$. $I(\tilde{t}_{1},\tilde{t}_{2})=I(\tilde{T},\tilde{t})$. Taking partial derivative over Eqn.$\tilde{t}$ in 37 yields,

[TABLE]

We can integrate over $\tilde{\omega}_{2}$ in the first term in Eqn. 38, and integrate over $\tilde{\omega}_{1}$ in the second term in Eqn. 38.

[TABLE]

where $J_{1}$ is the Bessel function and we always let $\tilde{t}_{1,2}\geq 0$.

Let me explain in more detail how to perform the square-root integral,

[TABLE]

where $\theta(t)$ is the Heaviside function and $K_{1}$ is the Bessel function of the second kind.

Using Eqn.42, we can get the result in Eqn.39. Integrate over $\tilde{t}$ in Eqn.39, we can get,

[TABLE]

where $[_{1}F_{2}]$ is the hypergeometric function.

Next, we want to calculate $I(\tilde{T},0)$.

[TABLE]

Integrate over $\tilde{T}$, one can get,

[TABLE]

Note that $\lim_{t_{1}\to\infty,t_{2}\to\infty}F(t_{1},t_{2})=0$ by Riemann-Lebesgue lemma. Here, we approximate $I(\tilde{T},0)\approx 0$ in order to match the long time behavior of $F(t_{1},t_{2})$.

Finally, we obtain the approximate expression for the squared-commutator,

[TABLE]

IV Self-energy at $O(1/M)$ and decay rate for OTOC

As discussed in the main text, the leading correction to the self-energy at $O(u^{2}/M)$ from the quartic interaction term (Fig. 2) is given by,

[TABLE]

The above self-energy correction gives rise to an exponential decay of the Green’s function and to the OTOC. The decay rate can be estimated in a straightforward fashion by looking for the imaginary solutions for the following equation,

[TABLE]

This yields, $\omega=\pm\Delta_{s}-i\Gamma$, where

[TABLE]

as quoted in the main text. The negative imaginary part gives rise to an exponential decay of $F(t_{1},t_{2})$.

V Numerical analysis of $1/M$ corrections

In this subsection, we describe details of our numerical evaluation of the ladder sum for $F(t,t)$. The Bethe-Salpeter equation for the $1/M$ corrected squared-commutator $F_{u}(t_{1},t_{2})$ can be written as,

[TABLE]

where $F_{d}(t_{1},t_{2})$ represents the contribution to the ladder sum without the $O(u^{2}/M)$ correction to the rung, but includes the self-energy at $O(u^{2}/M)$ in the dressed propagators. As noted earlier, $F(t_{1},t_{2})$ in the large $N,~{}M$ limit does not have any exponential growth; an explicit analytical form for $F(t_{1},t_{2})$ appears below. Upon including the $O(u^{2}/M)$ self-energy correction into account, $F_{d}(t_{1},t_{2})=F(t_{1},t_{2})e^{-\Gamma(t_{1}+t_{2})}$, where we evaluated $\Gamma\sim(u^{2}/M)\Delta_{s}^{4}$. Finally, we have defined $P(t)=[Q_{W}(t)]^{2}$.

The above self-consistent integral equation can be viewed as solving a matrix inversion problem once we rewrite it as,

[TABLE]

In more explicit terms, the integral on the left hand side can be discretized as a summation on a fine (time) grid, $[A_{t_{1}t_{2};t_{3}t_{4}}F_{u;t_{3}t_{4}}]=F_{d;t_{1}t_{2}}$, such that computing the inverse $A_{t_{1}t_{2};t_{3}t_{4}}^{-1}$ will lead us to the required form of $F_{u}$.

In our numerical calculations, we focus exclusively on the exponential piece of $F_{u}(t,t)=F_{0}(t)~{}e^{\lambda_{L}t}$, where $F_{0}(t)$ is an undetermined function of time and $\lambda_{L}>0$ is the Lyapunov exponent. At late times, our numerical analysis is consistent with an exponent, $\lambda_{L}=a\frac{T^{4}}{\sqrt{\log(1/T)}}$, where $a$ is a temperature independent ‘fitting’ parameters. However, there is some uncertainty associated with observing the $1/\sqrt{\log(1/T)}$ piece in $\lambda_{L}$ in our numerical analysis.

The numerical result for $F_{u}(t,t)$ is summarized in Fig.5 for different temperature.

VI $1/M$ corrections to the random matrix model

The quartic interaction between the rotor fields can be written as,

[TABLE]

where $\varphi_{i}^{p}$ represents an orthogonal matrix.

Let us consider the $1/M$ correction to the Green’s function for the $\boldsymbol{V}_{p}$ field due to the above interaction term. At one-loop order, this leads to a contribution $O(u/M)$ which renormalizes the ‘mass’ with no frequency dependence and can be absorbed into a redefinition of $\lambda$. At two-loop order, the contribution of the interaction to the self-energy is,

[TABLE]

where $\widetilde{Q}(\lambda_{p})$ is the Green’s function for $\boldsymbol{V}_{p}$, as introduced earlier, and we have suppressed the frequency dependence above.

In order to carry out an averaging over the disorder distribution, we now need to know the distribution of both the eigenvalues, $\lambda_{p}$ and the eigenvectors, $\varphi_{i}^{p}$; for the GOE of random-matrices, these distributions are known to be independent Anderson et al. (2010). This allows us to integrate over the eigenvalues first. (Note that the bare Green’s function for one component of $\boldsymbol{V}_{p}$ has a factor $1/M$ such that the leading $1/M$ correction has a factor $1/M^{2}$.) Including the outer legs, we have,

[TABLE]

Averaging over the eigenvalue distributions, the $p$ dependence of eigenvalues drop out and one can perform the summation over eigenvectors,

[TABLE]

where $R_{5}(...)$ is the 5-level correlation function in GOE. We can express the 5-level correlation function more generally as,

[TABLE]

The first term above is a simple (independent) product of the density of states for the five eigenvalues. The correlation among the different eigenvalues is contained in the $2,~{}3,~{}4$ and $5-$level “cluster functions”, $T_{2},~{}T_{3},~{}T_{4},~{}T_{5}$, respectively Anderson et al. (2010). The summation $P$ is the permutation among $x_{1}$, $x_{2}$, $x_{3}$, $x_{4}$ and $x_{5}$. The first term in $R_{5}(x_{1},..,x_{5})$ above leads to the $\sim u^{2}|\omega|^{5}$ singular structure in the imaginary part of the self energy as in the main text since each integration $[\int~{}d\lambda_{i}~{}\sigma_{J}(\lambda_{i})\widetilde{Q}(\lambda_{i})]$ gives a factor of $Q(i\omega_{n})$ (the rest follows the discussion in the main text). The remaining terms in Eqn. 61, that take into account correlations among eigenvalues of the random matrix, are not included in the replica treatment of the large$-N,~{}M$ saddle-point action for the rotor-theory.

Within the random-matrix picture, the higher-order corrections to the “free” theory can be studied systematically by introducing the higher level correlation functions (i.e. $R_{n}(x_{1},...,x_{n})$) for GOE Anderson et al. (2010). As is clear from the above discussion, there is a piece for all $R_{n}(x_{1},...,x_{n})$ which corresponds simply to an independent product of the $\sigma_{J}(x_{i})$; these are the terms that are also included in the replica action for the rotor theory.