The Magnetic Maze: A System With Tunable Scale Invariance
Tian-Gang Zhou Michael Winer Brian Swingle

Quasi-periodic X-ray eruptions years after a nearby tidal disruption event
    
M. Nicholl^10000-0002-2555-3192,
D. R. Pasham^20000-0003-1386-7861,  
A. Mummery^3,
M. Guolo^40000-0002-5063-0751, 

K. Gendreau^5,
G. C. Dewangan^6,
E. C. Ferrara^7,8,50000-0001-7828-7708,
R. Remillard^2, 

C. Bonnerot^9,10,
J. Chakraborty^20000-0002-0568-6000,
A. Hajela^11,
V. S. Dhillon^12,13,
A. F. Gillan^10000-0003-4094-9408,  

J. Greenwood^1,
M. E. Huber^14,
A. Janiuk^150000-0002-1622-3036,
G. Salvesen^160000-0002-9535-4914,
S. van Velzen^17, 

A. Aamer^1,
K. D. Alexander^18,
C. R. Angus^1,
Z. Arzoumanian^5,
K. Auchettl^19,20,  

E. Berger^21,
T. de Boer^14, 
Y. Cendes^21,22,
K. C. Chambers^14,
T.-W. Chen^230000-0002-1066-6098, 

R. Chornock^24,
M. D. Fulton^1,
H. Gao^14,
J. H. Gillanders^25,
S. Gomez^26, 

B. P. Gompertz^9,10,
A. C. Fabian^27,
J. Herman^14,
A. Ingram^28,
E. Kara^2, 
 
T. Laskar^29,30,
A. Lawrence^31,
C.-C. Lin^14,
T. B. Lowe^14,
E. A. Magnier^14, 

R. Margutti^24,
S. L. McGee^9,10,
P. Minguez^14,
T. Moore^1,
E. Nathan^32,  

S. R. Oates^33,
K. C. Patra^24,
P. Ramsden^1,9,100009-0009-2627-2884,
V. Ravi^32,
E. J. Ridley^9,10, 

X. Sheng^1,
S. J. Smartt^25,1,
K. W. Smith^1,
S. Srivastav^25,1,
R. Stein^340000-0003-2434-0387,  

H. F. Stevance^25,1,
S. G. D. Turner^350000-0002-8641-7231,
R. J. Wainscoat^14, 

J. Weston^1,
T. Wevers^26,
D. R. Young^1 

^ 1Astrophysics Research Centre, School of Mathematics and Physics, Queens University Belfast, 
 Belfast BT7 1NN, UK

^ 2Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 
 Cambridge, MA, USA

^ 3Oxford Theoretical Physics, Beecroft Building, Clarendon Laboratory, Parks Road, Oxford, 
 OX1 3PU, UK

^ 4Department of Physics and Astronomy, Johns Hopkins University, 3400 N. Charles St., 
Baltimore MD 21218, USA
 
^ 5NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA
 ^ 6Inter-University Centre for Astronomy and Astrophysics (IUCAA), PB No.4, Ganeshkhind, 
Pune-411007, India
 
^ 7Department of Astronomy, University of Maryland, College Park, MD, 20742, USA
 
^ 8Center for Research and Exploration in Space Science & Technology II (CRESST II), 
NASA/GSFC, Greenbelt, MD 20771, USA

^ 9School of Physics and Astronomy, University of Birmingham, Birmingham B15 2TT

^ 10Institute for Gravitational Wave Astronomy, University of Birmingham, Birmingham  
B15 2TT

^ 11DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 155, 2200 Copenhagen, 
Denmark
 
^ 12Department of Physics and Astronomy, University of Sheffield, Sheffield, S3 7RH, 
United Kingdom 

^ 13 Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain
 

^ 14Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu HI  
96822, USA

^ 15Center for Theoretical Physics, Polish Academy of Sciences, Al. Lotnikow 32/46,  
02–668, Warsaw, Poland

^ 16Center for Theoretical Astrophysics, Los Alamos National Laboratory, Los Alamos,  
NM 87545, USA

^ 17Leiden Observatory, Leiden University,  
Postbus 9513, 2300 RA Leiden, The Netherlands 

^ 18Department of Astronomy and Steward Observatory, University of Arizona, 933 North 
Cherry Avenue, Tucson, AZ 85721-0065, USA 

^ 19School of Physics, The University of Melbourne, VIC 3010, Australia 

^ 20Department of Astronomy and Astrophysics, University of California, Santa Cruz, 
CA 95064, USA 

^ 21Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, 
MA 02138-1516, USA
 
^ 22Department of Physics, University of Oregon, Eugene, OR 97403, USA

^ 23Graduate Institute of Astronomy, National Central University, 300 Jhongda Road,  
32001 Jhongli, Taiwan

^ 24Department of Astronomy, University of California, Berkeley, CA 94720-3411, USA

^ 25Astrophysics sub-Department, Department of Physics, University of Oxford, Keble Road, 
Oxford, OX1 3RH, UK
 ^ 26Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

^ 27Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge  
CB3 0HA, UK

^ 28School of Mathematics, Statistics and Physics, Newcastle University, Herschel Building, 
Newcastle upon Tyne, NE1 7RU, UK

^ 29Department of Physics & Astronomy, University of Utah, Salt Lake City, UT 84112, USA

^ 30Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, 6500 GL,  
Nijmegen, The Netherlands

^ 31Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, 
Edinburgh EH9 3HJ, UK
 
^ 32Cahill Center for Astronomy and Astrophysics, California Institute of Technology, 
Pasadena, CA 91125, USA

^ 33Department of Physics, Lancaster University, Lancaster LA1 4YB, UK

^ 34Division of Physics, Mathematics, and Astronomy, California Institute of Technology,  
Pasadena, CA 91125, USA

^ 35Department of Applied Mathematics and Theoretical Physics, University of  
Cambridge, Wilberforce Road, Cambridge, CB3 0WA, UK



    

===================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






§ INTRODUCTION


The fate of randomness in physical systems has long been an intriguing topic. There are many well-known examples that consider quantum particles moving in random scalar or vector potentials. Such randomness often leads to non-ergodic dynamics, manifesting in behaviors like glassiness and localization. One classic model is the p-spherical model <cit.>, in which a quantum particle moves in high dimensions in a random scalar potential. This is a standard solvable model of glassy physics with quenched disorder. Another related model is Parisi's hypercube model <cit.>, in which a particle jumps among the vertices of a high-dimensional hypercube with random magnetic flux through each face. The hypercube model was motivated by random Josephson junctions. 
Still another class of models is motivated by various two-dimensional devices with magnetic impurities that exhibit spatially random vector potentials. These models typically result in electron localization physics <cit.>, although these two-dimensional results are not easily generalized to higher dimensions. 

However, certain models can evade non-ergodic behavior and exhibit novel forms of low-energy dynamics. A well-known example is the Sachdev-Ye model, along with its extension, the Sachdev-Ye-Kitaev model <cit.>. Recently, a spin-based version of the SYK model was studied as part of a search for bosonic models with possible SYK-like behavior <cit.>. Although there is ample evidence that these systems can avoid a glass-like phase with replica symmetry breaking at low temperature, we wanted a system that would be manifestly replica symmetric. We then wanted to understand the potential for exotic quantum dynamics that might emerge at low energy.

Hence, in this paper we study a system—the “magnetic maze” (MM)—consisting of a single non-relativistic charged quantum particle moving in many spatial dimensions under the influence of a random magnetic vector potential. The vector potential is time-independent but varies randomly in space with translation-invariant spatial correlations. Focusing on a special case of power-law decaying vector potential correlations, we study the thermodynamics, real-time dynamics, and Lyapunov physics as a function of temperature and coupling. The large spatial dimensionality, N, provides a control parameter for our analysis. Our findings are summarized as follows; see also Figure <ref>.



At high temperature, the system is effectively classical, and we do not expect the physics to depend sensitively on the precise choice of vector potential correlations. We show that its thermodynamics are those of a free particle and that it exhibits chaotic dynamics in which the large-scale motion of the particle is diffusive, velocity auto-correlations decay exponentially in time, and there is a positive Lyapunov exponent[These results are obtained primarily for a specific choice of vector potential correlations discussed below. We consider a more general class of correlations in Appendix <ref>, primarily in the Euclidean context. Within this more general class, we have also studied a few examples of real-time dynamics in which we obtained the same signatures of chaotic dynamics.].

At low temperature, quantum effects manifest and the physical properties depend sensitively on the choice of vector potential correlations. We show that, for a special choice of power-law spatial correlations, the system develops an emergent scaling symmetry at low energy. Whereas the high-temperature thermodynamics is identical to that of a particle with energy, E, related to momentum, p, by the usual non-relativistic dispersion, E ∼ p^2, the low-temperature thermodynamics is instead determined by an effective dispersion, E ∼ p^z with a dynamical exponent z which is tunable via the strength of the magnetic field.

The low temperature dynamics is also chaotic, and for the same special choice of power-law spatial correlations, the emergent scaling symmetry determines the dynamical properties. In particular, there is a single timescale β = k_B T/ħ and a single length scale ℓ_β∼β^1/z which set the dynamical properties. The decay rate of the velocity auto-correlation function and the Lyapunov exponent determined from an out-of-time-order correlator are both of order β^-1, and the diffusivity is of order D∼ℓ_β^2/β∼β^2/z-1. Henceforth, we set k_B = ħ = 1; Appendix <ref> discusses the physical scales in the MM.

For other choices of vector potential correlations, the low-energy physics changes significantly. For example, when the spatial correlations decay more rapidly than the special case with scaling symmetry, the low temperature properties are still chaotic but more closely resemble those of a weakly perturbed free theory, with a relaxation rate and chaos exponent that vanish more rapidly than β^-1 at large β. We investigate this case and the opposite case of uniform magnetic field (corresponding to very long-ranged correlations) in Appendices <ref> (uniform field) and <ref> (more rapid decay of correlations).

These results are obtained via large N path integral methods <cit.>. We formulate the system on a thermal contour to access thermodynamics, on a Schwinger-Keldysh contour <cit.> to access dynamics, and on a contour with two time-folds to access out-of-time-order correlations and Lyapunov physics. Because we use large N as a control parameter, our results correspond to taking the limit N →∞ first before the low temperature limit. It is of course interesting to study 1/N corrections and the finite dimensional fate of the phenomena we find, but we do not pursue that here.

Central to our study is the result, discussed below, that the MM cannot have static correlations between different replicas and thus cannot have an equilibrium glass phase. The dynamics which emerges at low energy is instead that of a scale-invariant chaotic quantum system, one which features a tunable dynamical exponent. Moreover, at the lowest temperature and largest couplings we were able to access, we found a chaos exponent that comes within 15 % of the MSS bound <cit.>. We have not been able to analytically show that the model becomes maximally chaotic at infinite coupling, but we see no evidence of sub-maximal saturation in the chaos exponent. Hence, this system may provide a new example of maximal chaos. 










In the remainder of the introduction, we will introduce the model in more detail and give a guide to the structure of the paper. We will also throughout the paper compare and contrast our findings with the properties of the Sachdev-Ye-Kitaev model. 



 §.§ Setup and Overview





The Magnetic Maze consists of a single particle in N-dimensional space subject to a random vector potential A_i(x). In order to keep our partition functions finite, we will sometimes consider a confining potential V_conf=ϵ/2∑_i x_i^2, but most of our results are most interesting when ϵ is far smaller than any other scale. The Hamiltonian is thus

    H=∑_i=1^N [ 1/2m(p_i-A_i(x))^2+ϵ/2x_i^2 ].

The elements of A are picked from a Gaussian distribution with mean zero. The covariance of A falls off as A_i(x)A_j(y)=f(|x-y|^2/N) δ_ij. We are particularly interested in the case where f(x^2/N)∼1/x^2/N for large x, which is the special power-law form highlighted in the introduction. We take the function to be

    f(x^2/N)=J^2/ℓ^2+x^2/N,
 
where ℓ is some short-distance scale and J is a dimensionless measure of the field strength. 

Although our focus is (<ref>), we note that the model makes sense for other choices of f. We address some of these in Appendices: the case of a uniform field (Appendix <ref>), which is an integrable system with f(x^2/N)=C-1/2^2x^2/N, and the case of power-law correlations with different power-law exponent at large x  (Appendix <ref>). As we show, (<ref>) is the most interesting choice because it leads a host of interesting zero- and low-temperature properties including emergent scaling symmetry and near-maximal chaos. By contrast, the high temperature properties are less sensitive to the form of f.

This paper focuses on studying the classical and quantum behavior of (<ref>) with f given by (<ref>). Our main quantity of interest will be Δ(t_1,t_2)=1/N∑_i𝒯 (x_i(t_1)-x_i(t_2))^2, which can be interpreted as the squared distance traveled between t_1 and t_2, and we will study both the imaginary-time and real-time behavior of this quantity. The physics is conveniently discussed in terms of T/J and J as sketched in Figure <ref>. Both quantities are dimensionless in units where ℓ, m, and ħ are set to one (see Appendix <ref>). T≫ J is the classical regime and T ≪ J is the quantum regime. Moreover, J ≪ 1 is weak coupling and J ≫ 1 is strong coupling.  

At zero temperature, T=0, we find the imaginary time-dependence Δ(τ) ∼τ^α with α determined by J through

    2tanπα/2πα/(α+1)=J^-2

which takes us from the free value α=1 at J^2=0 to the confined α=0 as J^2 goes to infinity. In terms of the previously introduced dynamical exponent, we have z= 2/α. This identification arises because Δ scales like two powers of length and length scales like time to the 1/z power.

Moreover, we find that the system has reparameterization symmetry like the SYK model, and power-law decay in the velocity correlations. There is also a rescaling symmetry in space, completely separate from the rescaling in time included in the reparameterization mode. This is in contrast with the SYK model, where the fermion modes stretch in a way controlled by the reparameterization mode.

Turning on a non-zero but low temperature, T ≪ J, the thermodynamics is controlled by the modified dispersion with non-trivial dynamical exponent, leading to a heat capacity equal to N/z = N α/2. For the real time dynamics, we find that the system exhibits properties characteristic of a quantum chaotic system. The velocity auto-correlation now decays exponentially in time with a time-scale set by β = 1/T. Out-of-time-order correlators also feature an initial period of exponential growth with chaos exponent κ∼ 1/β. As J increases, κ gets closer to the chaos bound, κ^max = 2π/β. Both the relaxation time and the Lyapunov time depend on J, so we have in fact an infinite family of ensembles of non-relativistic scale-invariant quantum systems.

At higher temperatures, T ≫ J, the system becomes effectively classical and the thermodynamics reduces to that of a free particle, corresponding to z=2 and α=1. The real time dynamics is still chaotic and we compute the diffusivity, the velocity relaxation time, and the largest Lyapunov exponent.

Our detailed analysis of the Magnetic Maze proceeds as follows. In Section <ref> we use the large-N limit to derive a mean-field action in terms of collective displacement variables. We then derive a saddle-point equation. Throughout Section <ref> we assume a general correlation function f for the As.

In Section <ref> we specialize to the situation in equation (<ref>), and work at low temperature. We derive the power-law behavior of the correlation function, as well as equation (<ref>). We study the low temperature thermodynamics, showing that the low-temperature heat capacity is α N/2, which is consistent with the picture of power-law dispersion with exponent E∝ p^z with z= 2/α. Finally, we show that at low temperatures and long times, the system has a reparameterization symmetry and a rescaling symmetry.

Section <ref> covers the real time dynamics of the system. We derive Schwinger-Dyson mean-field equations, and numerically solve them. We show that a particle in the Magnetic Maze exhibits diffusive transport, with the squared displacement after time t growing in proportion to t. More precisely, the displacement grows at 2Dt, where the diffusion constant D scales as β^α-1 at low temperature. In the same temperature regime, the velocity auto-correlation decays exponentially with a time scale set by β. We also note that the system has a somewhat unusual structure on the Schwinger-Keldysh (SK) contour, as we discuss in detail. 

Section <ref> covers the Lyapunov physics. We setup a ladder calculation of an out-of-time-order correlator (OTOC) and carry out a detailed numerical analysis. We find that the low temperature chaos exponent is of order the bound, 2π/β, and comes close to it as J increases.

Section <ref> discusses directions for further work and open questions. Following it are several technical appendices. 

Appendix <ref> provides a crash course on the physics of magnetic fields in high dimensions, analyzes the uniform field problem (which is integrable and exactly solvable at finite N), and explains how the dimensionless parameters T/J and J arise.

The other appendices are as follows. Appendix <ref> reviews the exact solution for J=0. Appendix <ref> discusses other choices of f. Appendix <ref> presents details of derivation of the real-time equations of motion. Finally, Appendix <ref> discusses the numerical method used to solve the equations of motion.


       
























§ ACTION AND EQUATIONS OF MOTION



We are going to perform a mean-field calculation in which the collective variable

    G(_1,_2)=1/N𝒫∑ x_i(_1)x_i(_2)

will play a central role. Here, 𝒫 denotes an appropriate path- or time-ordering as specified below. Our goal is to find a mean-field expression for the free energy of the system, expressed as a path integral over G. 

We start with the partition function

    Z=∫exp(∫_0^β-m/2 v^2-iA· v-ϵ/2 x^2d)𝒟 x.

We then introduce the Lagrange multiplier Σ enforcing (<ref>), 





    Z=∫exp(∫_0^β-m/2 v^2-iA· v-ϵ/2 x^2d-1/2∫_0^βΣ(_1,_2){ N G(_1,_2) -x(_1)· x(_2)}d_1 d_2)𝒟 x 𝒟 G 𝒟Σ.


Next we perform the average over As to get

    Z̅   =∫exp(-S)𝒟 x 𝒟 G 𝒟Σ
    
        S   =∫_0^βm/2ẋ^2+ϵ/2 x^2d+1/2∫_0^βΣ(_1,_2){NG(_1,_2) -x(_1)· x(_2)}+f((x(_1)-x(_2))^2)v(_1)· v(_2)d_1 d_2.

Note this leads to an “annealed” free energy ∝lnZ, whereas the more physical “quenched” free energy is ∝ln Z. However, as we will see shortly, these free energies always coincide at large N in the Magnetic Maze.

Using the delta function setting G to (<ref>), we can replace (x(τ_1) - x(τ_2))^2/N = Δ(τ_1,τ_2) where

    Δ(_1,_2)=G(_1,_1)+G(_2,_2)-2G(_1,_2).

Similarly, we can replace v(τ_1) · v(τ_2) with ∂_τ_1∂_τ_2 G(τ_1,τ_2). The resulting functional integral is now quadratic in the xs, and performing the x integral gives

    S/N=1/2log(-m∂_^2+ϵ-Σ)+1/2∫Σ G+f(Δ(_1,_1))∂__1∂__2G(_1,_2)d_1d_2.

This is the sought-for mean-field action which can be used to compute the free energy.












The equations of motion can now be obtained as functional derivatives of (<ref>) with respect to G and Σ. Using the functional derivative

    δΔ(_1,_2)/δ G(_1',_2') =  δ(_1 - _1')δ(_1 - _2') + δ(_2 - _1')δ(_2 - _2') - 2 δ(_1 - _1')δ(_2 - _2')

and the identity -2∂__1∂__2 G = ∂__1∂__2Δ, we get the following equations:

    (-m∂_^2+ϵ-Σ)*G=δ(_1-_2)
    Σ(_1,_2)=∂__1∂__2(f(Δ(_1,_2))+f'(Δ)∂__1∂__2Δ-δ(_1-_2)∫_0^βf'(Δ(_1,'))∂__1∂_'Δd'.




It is the last term in the equation for Σ that is most unfamiliar. It is a consequence of the fact that the interaction term depends not just on G(_1,_2), but also on G(_1,_1) and G(_2,_2). This special dependence of the interaction on G only through its derivatives and through Δ means that adding an overall constant to G (corresponding to translating the entire path by a random vector) doesn't affect the magnetic part of the action, a consequence of the statistical translational symmetry of our magnetic field. Another consequence is that Σ_ω=0 is exactly zero.

It is the first equation in (<ref>) that illustrates why we included ϵ. Without it, the zero-frequency mode of -m∂_^2+ϵ-Σ would be zero, and the equation wouldn't be satisfiable. The weak potential from ϵ saves us from a diverging position by limiting the motion of the particle to a sphere of squared radius N/ϵβ. This also cuts out divergences in the partition function and entropy due to the infinite volume of space, instead replacing them with a contribution to the entropy 


    Entropy from volume of box=-N/2log(βϵ).


It is also useful to write out equations (<ref>) assuming time-translation invariance. In this case, cutting out G entirely, we have

    Δ_ω≠ 0=-2/mω^2+ϵ-Σ_ω≠0
    Σ()=-2∂_^2 Δ f'(Δ)-(∂_Δ)^2 f”(Δ)+δ()∫_0^βf'(Δ('))∂_'^2Δ(') d'.


Equations (<ref>) illustrate why the quenched and annealed free energies are equal in the Magnetic Maze. The quenched free energy can be obtained from a replica trick involving n replicas of the system in the limit n → 0. The mean-field variables would be generalized to G_ab, Δ_ab, and Σ_ab, where a,b=1,⋯,n are replica indices. Replica symmetry breaking requires correlations between the replicas, but the inter-replica correlations are typically time-independent owing to the separate time-translation symmetry in each replica. The generalization of (<ref>) to the multi-replica case would relate inter-replica terms in the self-energy, Σ_a ≠ b, to time derivatives of the inter-replica terms in the displacement, Δ_a ≠ b. Since the inter-replica correlators are time-independent, we can conclude from the generalization of the second equation of (<ref>) that Σ_a ≠ b = 0. Self-consistency from the multi-replica generalization of the first equation of (<ref>) then implies that Δ_a ≠ b=0. Hence, multi-replica saddle points are always replica diagonal in the Magnetic Maze.



§ LOW-TEMPERATURE LIMIT WITH F∼J^2/Δ



In this section, we study the thermodynamics of the Magnetic Maze with f(Δ) ∼J^2/Δ and reveal its fascinating low-temperature physics. In the IR limit, where Δ is large and ∂_ is small, we can simplify our action. Equation (<ref>) becomes

    S/N=1/2log(Σ)-1/4∫( ΣΔ +J^2∂__1∂__2Δ(_1,_2)/Δ(_1,_1))) d_1d_2.

This IR action behaves in a simple way under two important classes of transformations: time reparameterizations and scaling transformations.

For reparameterizations of time, →σ(), the transformations are

    G(_1,_2)    → G(σ(_1),σ(_2)) 
    Δ(_1,_2)    →Δ(σ(_1),σ(_2)) 
    Σ(_1,_2)    →Σ(σ(_1),σ(_2)) |σ'(_1) σ'(_2)|.

It is straightforward to verify that this transformation preserves the last two terms in (<ref>). The lnΣ term is not manifestly symmetric, but by writing the transformation as

    Σ(_1,_2) →∫ d_1' d_2' M(_1,_1') Σ(_1',_2') M(_2',_2)

for a “matrix” M(_1,_1') = δ(_1' - σ(_1)) | σ'(_1)|, we can argue that the determinant is invariant by arguing that the determinant of M is unity.[For example, if σ is a time shift, σ()=+a, then M = δ(_1' - _1-a). This is effectively a permutation, which has determinant one. For a general infinitestisimal transformation, σ =  + ϵ(), one can directly show that the determinant is one by using the single-valuedness of ϵ().] 

For scalings, the time is unchanged but we send

    G    →λ G 
    Δ   →λΔ
    Σ   →λ^-1Σ.

We can also extend this transformation to the underlying x() variables as x() →√(λ) x(). The last two terms in (<ref>) are again manifestly invariant, but the first term transforms by a field-independent shift that depends on λ. As a consequence, the equations of motion in the IR limit are invariant under the scaling transformation but the action is not. This situation is slightly unusual but it can be understood as follows. Going back to (<ref>) before the x variables were integrated out, we see that the action is fully invariant under the scaling transformation provided we drop the corresponding terms (the x kinetic term, the confining potential, and the short-distance part of f), but the 𝒟x measure is not invariant (while the combined 𝒟G 𝒟Σ measure is). The scaling symmetry is also explicitly broken by the aforementioned terms we neglected in the IR limit.

Both of these transformations are approximate symmetries of the action in the IR limit, (<ref>), with the caveat about scalings discussed just above. However, we will see in the next subsection that the ground state spontaneously breaks these symmetries to a special subgroup.










 §.§ Power-Law Behavior of Ground State Correlations


We now show that a power-law ansatz solves the equations of motion arising from (<ref>) in the zero temperature limit. We consider a translation invariant ansatz in which Δ(_1,_2) = Δ(_1-_2) and similarly for Σ.

Suppose the IR behavior of Δ(τ) is 

    Δ(τ)= k|τ|^α.

Then the IR value of Σ is given by 

    Σ(τ)=-J^2 2α/k |τ|^-α-2.
 

This gives 

    Δ_ω=-2sinπα/2Γ(α+1)k |ω|^-α-1
    Σ_ω=2sinπα/2Γ(-α-1)J^2 2α/k |ω|^α+1,

as well as

    Δ_ωΣ_ω=-8sin^2 πα/2Γ(-α-1)Γ(+α+1)α J^2=-2.

We can use the reflection formula,

    Γ(z)Γ(-z)=-π/zsin(π z),

to get

    2tanπα/2πα/(α+1)=J^-2.


It turns out that for any J we pick, there is a unique 0<α<1 which solves equation (<ref>). We see in figure <ref> that equation (<ref>) aligns perfectly with numerics for reasonably small J^2 at a certain large fixed β. For larger J^2, one would need correspondingly larger β to see this agreement.


From the explicit form Δ∼τ^α, we learn that time reparameterizations are broken down to translations. Moreover, the scaling symmetry is also broken. However, the combination of a time reparameterization that rescales time and a scaling transformation leaves the form of Δ invariant. In this sense, the symmetries discussed in the previous subsection are spontaneously broken in the ground state to a subgroup of translations and combined rescalings. Of course, these symmetries are also explicitly broken by the non-IR terms in the full action.






 §.§ Imaginary Time Numerics


We saw in the previous subsection that at zero temperature, Δ() is proportional to ^α with a specific α given by the solution to the transcendental equation (<ref>). By analogy with the SYK model, one might expect the finite-temperature analog to be Δ()∝(sinπ/β)^α. In fact, we see from the graphs in Figure <ref> that for small J^2 the solution looks like

    Δ()∝((β-))^α.

When J^2=0 and α=1, this ansatz becomes exact (see appendix <ref>). At large J^2, the magnetic length becomes much smaller than ℓ. This maps onto the constant field limit discussed in Appendix <ref>, where the large-β solution is found to be proportional to logβsinπ/β for large β. This is not the small-α limit of equation (<ref>), suggesting that the function takes the more general form

    Δ()=(β F_α(/β))^α

where F_α is some α dependent function satisfying F_α(x)=F_α(1-x) (by the KMS condition <cit.>) and F_α(x)∼ x for small x. For α∼ 1 F_α(x) seems to take the form x(1-x), while for α∼ 0 we have F_α(x)=sinπ x.







 §.§ The Free Energy at Low Temperature

To cap off our discussion of thermodynamics, we will calculate the free energy of our system at low temperatures.

As already discussed in equation (<ref>), there is a contribution to the entropy of -N/2log(βϵ), and a corresponding contribution to the free energy of N/2β^-1log(βϵ), essentially the log-volume of the "box" the particle is in.

The most physically meaningful quantities are the ones that don't depend on the box volume: energy and heat capacity. Even more meaningful is kinetic energy. The potential energy is an artifact of our mathematically convenient choice to confine our system with a quadratic potential; other shapes such as quartic or infinite well would give different potential energies, but the kinetic energy would be invariant so long as the box is large.

One can see by examination that Hamiltonian (<ref>) is always positive. Furthermore, at high temperatures the As become negligible compared to the ps, and our particle becomes essentially classical. A classical particle in a magnetic field will have kinetic energy Nβ^-1/2 by the equipartition theorem.



    












At low temperatures, there is an equally natural interpretation of the thermodynamics. The imaginary-time distance function Δ()∼^α corresponds to a particle with dispersion E∼|p|^2/α. Working out the thermodynamics of such a particle, we see that the entropy goes as S=αN/2logβ+const, and the heat capacity is precisely αN/2. We see in Figure <ref> that the intensive heat capacity is indeed α/2 at low temperatures.










































§ REAL-TIME DYNAMICS



In this section, we study the real-time dynamics of the magnetic maze. We firstly set up the Schwinger-Keldysh formalism and derive the saddle point equations. Secondly, we solve these equations for specific vector potential distribution Eq. (<ref>). The system shows a universal long-time limit in which the particle motion is diffusive.  

        

 §.§ Schwinger-Keldysh Formulism

        

  §.§.§ Derivation of the Saddle Point Equations


        To study the two-point correlation function in real-time, it is convenient to represent the time-evolved density matrix using the Keldysh contour [Kamenev, Altland]
	
    Z = (U_-∞,∞ U_∞,-∞ρ_0)=∫𝒟 x  exp(∑_j ∑_a (-1)^a∫_-∞^∞ t ( m/2ẋ_a,j^2 +  A(x_a,j)·ẋ_a,j -ϵ/2 x_a,j^2 ) ),

        where a = +, - indicates the forward and backward contours, respectively, with the phase (-1)^a representing the direction of unitary evolution. By convention, we assume a = 0, 1 corresponds to the contour labels +, -. Here, ρ_0 = e^-β H is the thermal density matrix. For the real-time dynamics, the corresponding Green's function is defined as follows:
	
    G_ab(t_1,t_2)=1/N𝒫∑_i x_a,i(t_1)x_b,i(t_2).

	
	We introduce the Lagrange multiplier Σ_ab to enforce Eq. (<ref>)

    Z   =∫𝒟 x  exp(∑_j ∑_a (-1)^a∫_-∞^∞ t ( m/2ẋ_a,j^2 + A(x_a,j)·ẋ_a,j -ϵ/2 x_a,j^2 ) 
        - ∑_ab (-1)^a+b1/2∫ t_1  t_2 [ NΣ_ab(t_1,t_2)G_ab(t_1,t_2) + ∑_j Σ_ab(t_1,t_2)x_a,j(t_1)· x_b,j(t_2) ] )

	Similar to the imaginary-time derivation, we average over the A fields to obtain
	
    Z̅=    ∫exp( S)𝒟 x
     S/N =   ∑_a (-1)^a∫ t_1  t_2  x_a(t_1) ( -m/2∂_t_1∂_t_2δ(t_1-t_2)δ_ab  -ϵ/2δ(t_1-t_2)δ_ab - ∑_b(-1)^bΣ_ab(t_1,t_2)/2) x_b(t_2) 
        - ∑_ab (-1)^a+b∫ t_1  t_2  ( 1/2Σ_ab(t_1,t_2)G_ab(t_1,t_2) + 1/2  f((x_a(t_1)-x_b(t_2))^2/N) ẋ_a(t_1)·ẋ_b(t_2)/N )

	
	
        We can perform the quadratic integral over the x_a, which leads to the effective action.
	
    S/N=     - 1/2log( (-1)^a(- m ∂_t_1∂_t_2-ϵ)δ(t_1-t_2)δ_ab - (-1)^a+bΣ_ab)  
       - ∑_ab (-1)^a+b∫ t_1  t_2  ( 1/2Σ_ab(t_1,t_2)G_ab(t_1,t_2) + 1/2 f(Δ_ab(t_1,t_2) ) ∂_t_1∂_t_2 G_ab(t_1, t_2) )

	where we defined 
	
    Δ_ab(t_1,t_2)=(G_aa(t_1,t_1)+G_bb(t_2,t_2)-G_ab(t_1,t_2) - G_ba(t_2,t_1))

	to be the squared distance between the positions at t_1 and t_2. In Eq. (<ref>), all fields are defined on the original Keldysh contour. It is convenient to obtain the relation between the self-energy and Green's function by using the first set of saddle point equations, ∂ S/∂ G_ab = 0, which reads
    
    Σ_ab(t_1,t_2) =-∂_t_1∂_t_2f(Δ_ab(t_1,t_2)) -  f'(Δ_ab(t_1,t_2))∂_t_1∂_t_2Δ_ab(t_1,t_2) + δ(t_12)δ_ab∫t_3(   
          f'(Δ_ab(t_1,t_3))∂_t_1∂_t_3Δ_ab(t_1,t_3) -1/2 f'(Δ_+-(t_1,t_3))∂_t_1∂_t_3Δ_+-(t_1,t_3) - 1/2 f'(Δ_-+(t_1,t_3))∂_t_1∂_t_3Δ_-+(t_1,t_3) )

    where G_ab≡[  G_T  G_<;  G_> G_T̃;      ]_ab and the same definition applies for Σ_ab. Here, T and T̃ denote time-ordering and anti-time-ordering, respectively, while <,> label the lesser and greater Green's functions, as the two time variables are on different contours.
    
    
    The Schwinger-Dyson equations are obtained by computing ∂ S/∂Σ_ab, with a,b either in the basis +,- or cl,q. In the basis x_cl, x_q after the Keldysh rotation, the Schwinger-Dyson equations read
    
    [                         0 -δ(t_12) (m∂_t_2^2+ϵ)-Σ_A; -δ(t_12) (m∂_t_2^2+ϵ)-Σ_R                      -Σ_K;                           ]∘[ G_K G_R; G_A   0;     ] = 𝕀

	Here, the ∘ symbol denotes convolution: A(t_1, t_2) ∘ B(t_2, t_3) ≡∫ t_2 A(t_1, t_2) B(t_2, t_3). It is also useful to express Eq. (<ref>) under the assumption of time-translation invariance and perform a Fourier transformation. The detailed derivation is provided in Appendix <ref>, and the self-energy is summarized in Table <ref>. We obtain the retarded component of the Schwinger-Dyson equation as follows: 

    (m(ω+ η)^2 - ϵ-(Σ̅_R(ω)-Σ̅_R(ω=0)))  G_R(ω)=1
       Σ̅_R(t)= Θ(t)(Σ_>(t)-Σ_<(t))

    The central result from the derivation is that Σ_R(ω) = Σ̅_R(ω) - Σ̅_R(ω=0). This exactly satisfies the condition Σ_R(ω=0)=0, as expected from the imaginary time calculation. The retarded Green's function generally requires a small cutoff η to retain the retarded causal structure. In principle, η should be taken as an infinitesimal number, η = 0^+, but a finite η is necessary in numerics to help with the convergence of the system, thus introducing artificial dissipation. In this case, we take ϵ to be zero while keeping the retarded cutoff η. Consequently, we focus only on the real-time Green's function where the time variable t ≪ 1/η to eliminate unphysical artifacts. This sets a lower bound for the low-temperature region, requiring the inverse temperature β≪ 1/η.
	
	To self-consistently solve the equation, we also require the fluctuation-dissipation theorem to determine the relation between G_R and G_≷,
    
    ρ_G(ω)    = -1/2π (G_R(ω) - G_R(ω)^†) 
    
    			G_<(ω)    = -2π n_B(ω) ρ(ω) 
    
    			G_>(ω)    = +2π n_B(-ω) ρ(ω),

    where n_B(ω) = 1/(e^βω - 1) is the bosonic distribution function and β is the inverse temperature.


Note that Green's function G is not well-defined without a trapping potential ϵ. An example of this is the analytical solution for free particles at J=0, as shown in Appendix <ref>. For the real-time Green's functions G_≷(t) and G_K(t), as well as the imaginary-time-ordered Green's function G(τ), all of them diverge in the limit ϵ→ 0. However, Δ_≷(t), Δ_K(t), or Δ(τ), which represent the squared displacement, are well-defined. Therefore, it is convenient to relate everything to the correlation function Δ.
From the definition of Δ_ab in Eq. (<ref>), we can assume time-translation invariance and transform to the frequency domain. Using the symmetry Δ_ab(t) = Δ_ba(-t), we obtain:

    Δ_ab(ω) =  (const.) δ(ω) -2  G_ab(ω).

Hence, up to a constant shift at zero frequency, Δ_ab and G_ab(ω) are related by -2. The constant shift can be ignored since we always enforce the condition Δ_ab(t=0)=0 by taking Δ_ab(ω = 0) = - ∑_ω≠ 0Δ_ab(ω) in the numerics.
Consequently, we can define ρ_Δ = - 1/2π(Δ_R(ω) - Δ_A(ω)) for ω≠ 0, and the corresponding relation between the retarded and advanced components is given by

    Δ_R^*(ω) ≡ -Δ_A(ω), Δ_R^*(t) ≡ -Δ_A(-t).

Here, the * symbol denotes taking the complex conjugate, and the interchange of time variables has been carried out explicitly. The spectral function ρ_Δ also satisfies the corresponding fluctuation-dissipation theorem for the correlator Δ. 

In summary, we have the full set of equations for Δ in Eq. (<ref>)-(<ref>). These can, in principle, be self-consistently solved to yield all real-time Green's functions by providing an initial guess for Δ_R(ω).
[0.98]

    ρ_Δ(ω)    = -1/2π (Δ_R(ω) - Δ_A(ω)),  ω≠ 0 
    Δ_<(ω)    = -2π n_B(ω) ρ_Δ(ω),  ω≠ 0 
    Δ_>(ω)    = +2π n_B(-ω) ρ_Δ(ω),  ω≠ 0 
    Σ_≷(t)    = ∂_t^2 f(Δ_≷(t)) +  f'(Δ_≷(t))∂_t^2 Δ_≷(t) 
    Σ̅_R(t)    = Θ(t) (Σ_>(t)- Σ_<(t)) 
    (m(ω+   η)^2 - ϵ-(Σ̅_R(ω)-Σ̅_R(ω=0) ) )  Δ^R(ω)=-2 ,  ω≠ 0





  §.§.§ Numerical Methods to Obtain the Solution

Here we explain how the real-time equations can be solved numerically. Since the value of Δ_≷(ω = 0) is undetermined in the self-consistent process, using the simple mixing method Δ_R^(n)=(1-ζ) Δ_R^(n-1)+ζΔ_R,new^(n-1) is numerically unstable. Here, n is the iteration step and Δ_R,new^(n-1) is calculated through the Schwinger-Dyson equation Eq. (<ref>) based on the solution Δ_R^(n-1). Instead, we choose to use a gradient descent protocol to achieve self-consistency.

If we only consider the retarded component, we will find that the gradient being zero is the same as the Schwinger-Dyson equation Eq. (<ref>). Leaving the details to the supplementary material, we choose the update program to be:

    Δ(ω)_R^(n)) = Δ(ω)_R^(n-1) + ζ/4( (-2)(Δ_R^(n-1)))^-1(ω) - (G_0,R^-1(ω) - Σ_R^(n-1)(ω))  ) (Δ_R^(n-1)(ω) )^2

Here, the second term /4( (-2)Δ_R^-1(ω) - (G_0,R^-1(ω) - Σ_R(ω))  ) = ∂ S[Δ]/∂Δ_R is the gradient of the action with respect to the variable Δ_R. We also introduce an extra Δ_R(ω) as a numerical trick, which stabilizes the iteration process. This is because Δ_R^-1, Σ_R, and G_0,R^-1 all approach zero as ω→ 0, and Δ_R^2 enhances the numerical difference around ω≈ 0. We have checked that convergence is reached when ||Δ_R^(n)(ω)-Δ_R^(n-1)(ω)||_2 < 10^-6 ||Δ_R^(n)(ω)||_2.


After obtaining the numerical solution, we can check the real-time dynamics by directly comparing them with the imaginary time dynamics as a benchmark. The correlator can be related to the spectral function by the spectral representation

    Δ_R(ω) = ∫ω'ρ_Δ(ω')/ω' - (ω + η).

In real-time dynamics, we can directly sample the imaginary time Green's function through the analytical continuation

    Δ(ω_n) = ∫ω'ρ_Δ(ω')/ω' - ω_n,

with the Fourier transformation

    Δ(τ) = ∑_ω_n=2π n/βΔ(ω_n) e^ω_n τ.













 §.§ Long-Time Limit of The Dynamics




Here we describe the structure of the real-time dynamics in equilibrium at inverse temperature β. We demonstrate two typical numerical results in Fig. <ref>, with (a1,a2,c1,c2) in the high temperature regime and (b1,b2,d1,d2) in the low temperature regime. From this data, it is apparent that the equilibrium real-time dynamics exhibits a rich structure. Here we show that, for all parameters we studied, the dynamics eventually enters into a late-time “hydrodynamic” regime in which the large scale motion is diffusive, Δ_<(t) →  2D t, with D the diffusivity. Moreover, the system relaxes to this diffusive behavior exponentially with a time constant t_rel. The time-scale for the onset of the late-time hydrodynamic regime is also of order t_rel. 

Before proceeding, we note that the dynamics at early time is quite complex. For example, the Euclidean scaling form, Eq. (<ref>), has an imprint on the early time dynamics at very low temperature (see (b1,b2) of Fig. <ref>). Nevertheless, we focus on the long-time dynamics as it is the most universal and it is closely related to existence of a non-vanishing chaos exponent (discussed in Sec. <ref>).

The main goal of the rest of this subsection is to explain our results for D and t_rel at very high and very low temperature. At low temperature, we find t_rel∼β and D ∼β^2/z-1. These low temperature forms are determined by the emergent scaling symmetry and the associated dynamical exponent z, but we have not been able to give an analytic derivation of the coefficients. At high temperature, we find t_rel∼β^-1/2 J^-2 and D ∼β^-3/2 J^-2. These we will derive analytically below.

The values for D and t_rel just quoted, and the remainder of the discussion in this subsection, all consider f ∼ J^2/Δ. For other choices of f, we expect the system will still exhibit exponential relaxation to diffusive dynamics at long time (apart from the special case of uniform magnetic field or zero coupling). It is also interesting to consider the possibility that the dynamics could localize in some fashion, but we have not seen evidence for this in our numerics.




  §.§.§ Theory of Low-Temperature Parameters




Here we discuss the low-temperature parameters, relying heavily on the emergent scaling symmetry at low temperature. In Fig. <ref>, we show the numerics of real-time dynamics in the low-temperature region. We focus on two quantities to describe the late-time behavior.

    
  * The diffusion constant D is obtained by fitting the slope of the late-time correlator Δ_<(t) ∼ 2Dt + const. The fitting region is chosen such that ^2/ t^2Δ_<(t) is small enough, ensuring that the system has entered the diffusion region.
    
  * The relaxation times  are obtained by fitting the slope of log| ^2/ t^2Δ_<(t) | ∼ t    (or log| ^2/ t^2Δ_<(t) | ∼ t), as illustrated in the insets in Fig. <ref>(c1, c2) and (d1, d2). We can exactly prove that real and imaginary part fitting leads to the same relaxation time  in the high-temperature case, as indicated by Eq. (<ref>). However, in low-temperature case, the numerical extraction of  for the real and imaginary parts is slightly different, which we attribute to numerical limitations. We find that both methods give a relaxation time that scales linearly with β. Therefore, we extract  from the slope fitting of log| ^2/ t^2Δ_<(t) | ∼ t.
    




Remarkably, our model exhibits power-law scaling with D ∼β^2/z-1, as shown in Fig. <ref>(a). Here we use the exponent 2/z-1 to label the slope on the D, β log-log plot, where z is the dynamical exponent. This scaling arises from dimensional analysis of the diffusion behavior Δ_<(t) ∼ 2Dt + ⋯, where [Δ_<(t)] ∼ [⟨ x^2 ⟩] ∼ [t]^2/z, and therefore [D] ∼ [t]^2/z-1∼ [β]^2/z-1. The brackets, [], indicate the unit of the physical quantity.
We observe that the dynamical exponent is continuously tunable and depends on J. In Fig. <ref>(b), the fitted exponent 2/z from the diffusion constant generally shows quantitative agreement with α_Theory from the IR analysis in Eq. (<ref>). When J^2 ≤ 0.1, the error increases. This artificial numerical effect, also encountered in the high-temperature region, occurs because smaller J leads to larger t_, causing the system to potentially not enter the diffusion region during the fitting range t ∈ [0, 1/η], which results in a larger numerical error. In conclusion, in low-temperature regions, the dynamical exponent is precisely determined by the IR theory α_Theory with z = 2/α_Theory.

We are also interested in the relaxation time  to the diffusion region in the low-temperature system. In the numerical data of Fig. <ref>(c),  is obtained by linearly fitting log| ^2/ t^2Δ_<(t)| in the time region t∈ [1.4β, 1.55β], with an R-square value R^2>0.999. The fitting regions ensure that the fast decaying mode has already vanished and only one late-time decaying mode exists. Furthermore, we verify the linear behavior of ∼β by checking the R-square value R^2>0.9998 for all different J^2 in Fig. <ref>(c).
We can extract the coefficient /β (and an analogous quantity from fitting the imaginary part) for different J^2. The relaxation time is an O(1) value as a function of J^2. Both /β and the analogous quantity extracted from the imaginary part show similar behavior regarding J^2, while they differ by a constant value. With a larger temporal fitting window, we expect the time constants of the exponential decay of both the real and imaginary parts to converge. 

We can unify the behavior of these two quantities,  and D, at low temperatures within one framework. In the low temperature regime, the largest characteristic time scale is β. Therefore, we expect the relaxation time to be proportional to β since it has the unit of time. By virtue of the dynamical exponent relating space and time scaling, the diffusion constant has the unit of [time]^2/z-1, and therefore we expect that D ∼β^α - 1 using the relation 2/z=α.













  §.§.§ Theory of High-Temperature Parameters


Here we discuss the high-temperature parameters, relying on some useful approximations to give an analytical treatment. At high temperature, the approximation for Δ_<(t) ≈Δ_>(t) ≈t^2/m β works at a relatively large time scale t∈ [β, ], where the imaginary part of Δ_≷(t) can be ignored due to the condition | t/m |≪t^2/m β when t≫β. We obtain Σ_K by referring to the self-energy formula in Table <ref>,

    Σ_K(t) = Σ_<(t) + Σ_>(t)  = 8 β  J^2 m (t^2-l^2 β  m)/(l^2β  m+t^2)^3.

After Fourier transformation, the corresponding Σ_K(ω) is

    Σ_K(ω) = -2 i π  J^2 e^-l ω√(β  m)(l ω(β  l m ω +√(β  m))+1)/l^3 √(β  m).

In the small ω region, it can be expanded as the following equation

    Σ_K(ω) = -2 i π  J^2/l^3 √(β  m)-i πβ J^2 m ω ^2/l √(β  m) + O(ω^3).

We aim to obtain the relaxation time, so we consider the time scale when | ω | ^2≪(1/l^2β m), or effectively t^2≫ l^2β m. We can approximate Σ_K as a constant by dropping the quadratic term in ω,

    Σ_K(ω) ≈ - 2 i π  J^2 /l^3√(β  m).


By the fluctuation-dissipation theorem in Table <ref>, we obtain the retarded self-energy term,

    Σ_R(ω) = 1/2 Σ_K(ω) tanhβω/2≈ - βω/2π  J^2 /l^3√(β  m),

where we approximate tanh(βω/2) as βω/2 in the high temperature limit. In principle, the real part of Σ_R should be determined from the Kramers-Kronig relation ReΣ_R(ω) = 1/π𝒫∫ImΣ_R(ω')/ω' - ωω', but the analytic continuation can be read off as

    Σ_R(ω) ≈ - ωβ/2π  J^2 /l^3√(β  m).

Notice that Σ_R(ω)∝ω also automatically satisfies the condition Σ_R(ω=0)=0 shown in Table <ref>.

With all the approximations above, the retarded component of the Schwinger-Dyson equation in Eq. (<ref>) can be further simplified to

    (m(ω+ η)^2 - ϵ)  Δ_R(ω) + ωπ J^2 √(β)/2l^3√(m)Δ_R(ω) =-2 ,  ω≠ 0 .

In the time domain, we set ϵ to zero and Eq. (<ref>) becomes

    -m ^2 Δ_R(t)/ t^2  - π J^2 √(β)/2l^3√(m)Δ_R(t)/ t = -2δ(t).

We solve the ordinary differential equation Eq. (<ref>) with the boundary condition Δ_R(0^+)=0, Δ_R'(0^+)=2/m, Δ_R(0^-)=0, Δ_R'(0^-)=0. Here Δ_R(0^+)=0, Δ_R'(0^+)=2/m can be directly obtained from the free theory J=0. The solution for the retarded Green's function reads

    Δ_R(t) = Θ (t) 2/m(1-e^-|t|/),

where the relaxation time is

    = 2 m^3/2 l^3/πβ^1/2 J^2.

Similarly, the advanced component can be obtained as

    Δ_A(t) = Θ (-t) 2/m(1-e^-|t|/).


Next, Δ_K(t) can be evaluated using the FDT,

    tanhβω/2Δ_K(ω) ≈βω/2Δ_K(ω) = (Δ_R(ω) - Δ_A(ω) ) .







Using the high-temperature approximation, we can further replace ω as ∂_t and obtain Δ_K(t) from a first-order ordinary differential equation, 

    Δ_K(t)/ t=2/β(Δ_R(t)-Δ_A(t)).

Consequently, we obtain Δ_K(t) and then the real-part Δ_≷, which leads to the final result

    Δ_K(t)    = 4(e^-|t|/-1 + |t|/)^2 /(mβ) 
    Δ_≷(t)    = 2(e^-|t|/-1 + |t|/)^2 /(mβ)  ±(t) /m(1-e^-|t|/).

From this result, the system displays diffusive behavior at late time with Δ_<(t) ≈ 2Dt + C, with constant C=-2^2/(mβ). The diffusion constant D can be simply related to the  by

    D = /mβ =  2 m^2 l^3/πβ^3/2 J^2 .




To benchmark our theory, we directly perform numerical calculations based on the Schwinger-Dyson equation. In Fig. <ref>, we provide concrete evidence for our high-temperature solution.

Firstly, we can check that the high-temperature theory, Eq. (<ref>), perfectly matches the numerics, as indicated by the overlap between the dashed lines and solid lines in Fig. <ref>. Additionally, in the derivation of the self-energy, Eq. (<ref>), we assume the real part of Δ_<(t) is significantly larger than the imaginary part by using the free theory solution, and the numerics show that this approximation remains true for the exact dynamics.

Secondly, after verifying the functional form of Δ_<(t), we focus on the characteristic relaxation time , or equivalently, the diffusion constant D in the high-temperature region. In Fig. <ref>(b,c,d), the diffusion constant (or ) is obtained by directly fitting the formula, Eq. (<ref>), in the reliable time region t ∈ [0, 1/η]. We choose different cutoffs, η = [0.007, 0.01, 0.015], in the numerics and extrapolate the data to η = 0.

In Fig. <ref>(b), we quantitatively verify the numerically obtained diffusion constant against the theoretically predicted value in Eq. (<ref>). In the small β→ 0 limit, the ratio D_Fit/D_Theory≈ 1 within an error of 5%. We note that the error increases as J approaches zero. This is an artificial numerical effect because smaller J leads to larger t_, and the system might not enter the diffusion region during the fitting region t ∈ [0, 1/η], which causes larger numerical error.

In Fig. <ref>(c) and Fig. <ref>(d), we illustrate that the diffusion constant follows a power-law dependence on β and J by observing the linearity in the log-log plot. Most importantly, we extract the scaling law by comparing with D ∼β^-3/2 and D ∼ J^-2, which again confirms the correctness of our prediction in the high-temperature region as given by Eq. (<ref>).




§ LYAPUNOV PHYSICS



In this section, we setup the calculation of an appropriate out-of-time-order correlator (OTOC) <cit.> in order to extract a quantum Lyapunov exponent. We then report the results of a numerical study of the derived equations. In particular, at low temperature we find a chaos exponent which is proportional to 1/β and grows with increasing J. We do not observe any satuaration of this growth and the chaos exponent comes close to 90% of maximal at the largest Js we studied.




We study the 4-point function of xs by performing a perturbation around the saddle point of the full effective action Eq. (<ref>)

    G(τ_1,τ_2) = G̃(τ_1,τ_2) + δ G(τ_1,τ_2)  
    Σ(τ_1,τ_2) = Σ̃(τ_1,τ_2) + δΣ(τ_1,τ_2).


Then the effective action becomes

    I_eff[δ G,δΣ]/N =   1/4∫τ_1τ_2τ_3τ_4δΣ(τ_1, τ_2) ( G̃(τ_1, τ_3) G̃(τ_2, τ_4)  ) δΣ(τ_3, τ_4) + 1/2∫τ_1τ_2
       [δ G(τ_1,τ_2)δΣ(τ_1,τ_2) 
            + (δ G(τ_1,τ_2) )^2 f”(Δ̃)/2 (2(δ(τ_1-τ_2)-1))^2 ∂_τ_1∂_τ_2G̃(τ_1,τ_2) 
       + (δ G(τ_1,τ_2) ) f'(Δ̃) (2(δ(τ_1-τ_2)-1)) ∂_τ_1∂_τ_2(δ G(τ_1,τ_2) ) ]

After integrating out δΣ, the effective action reads

    I_eff[δ G]/N = ∫τ_1τ_2τ_3τ_4δ G(τ_1,τ_2) (1/4K^-1 - 1/2(I∘Λ) ) δ G(τ_3,τ_4),

where the kernel K is

    K(τ_1,τ_2,τ_3,τ_4) = -  G̃(τ_1, τ_3) G̃(τ_2, τ_4) ,

the interaction 2-point vertex induced by the magnetic field is

    Λ(τ_1,τ_2)    =  [ f”(Δ̃)/2 (2(δ(τ_1-τ_2)-1))^2 ∂_τ_1∂_τ_2G̃(τ_1,τ_2) +  f'(Δ̃) (2(δ(τ_1-τ_2)-1)) (∂_τ_1∂_τ_2) ],

and the identity 4-point vertex for the boson is

    I(τ_1,τ_2,τ_3,τ_4) = 1/2( δ(τ_1-τ_3)δ(τ_2-τ_4) + δ(τ_1-τ_4)δ(τ_2-τ_3)).

We can check that ∫τ_1τ_2τ_3τ_4δ G(τ_1,τ_2) ( 1/2(I∘Λ) ) δ G(τ_3,τ_4) is the same magnetic interaction term in Eq. (<ref>).

The 4-point function of x can be regarded as the 2-point function of the perturbed Green's function in the effective action Eq. (<ref>). We introduce

    F(τ_1,τ_2,τ_3,τ_4)    = ⟨δ G(τ_1,τ_2)  δ G(τ_3,τ_4) ⟩
    
               = 1/2(I - 2K ∘ (I∘Λ)))^-1 K∘I1/N

Then we can directly prove that the 4-point correlation function is generated by a series of ladder diagrams. To simplify the above equation we notice that 

    (K ∘I)(τ_1,τ_2,τ_5,τ_6)    = ∫τ_3τ_4K(τ_1,τ_2,τ_3,τ_4) I(τ_3,τ_4,τ_5,τ_6) 
       = - 1/2G̃(τ_1, τ_5) G̃(τ_2, τ_6)  - 1/2G̃(τ_1, τ_6) G̃(τ_2, τ_5)

where we can further define the decoupled 4-point function

    F_0(τ_1,τ_2,τ_3,τ_4)=-1/N(G̃(τ_1,τ_4)G̃(τ_2,τ_3)+G̃(τ_1,τ_3)G̃(τ_2,τ_4) ).

Combining all the results in Eq. (<ref>), (<ref>), (<ref>) we get the final result

    F    = (I - 2 K ∘ (I∘Λ)))^-1 F_0 
       = ∑_n=0(2 K ∘ (I∘Λ)) )^n F_0

Therefore 2 K ∘ (I∘Λ)) is the effective ladder diagram in our model, which reads 


    K_lad(τ_1,τ_2,τ_3,τ_4)    ≡(2 K ∘ (I∘Λ)(τ_1,τ_2,τ_3,τ_4))  
    
               = -2G̃(τ_1, τ_3) G̃(τ_2, τ_4) 
            [(f”(Δ̃)/2 (2(δ(τ_3-τ_4)-1))^2 ∂_τ_3∂_τ_4G̃(τ_3,τ_4) 
               + (f'(Δ̃) (2(δ(τ_3-τ_4)-1)) (∂_τ_3∂_τ_4) ] .



Finally, we consider the OTOC as a special kind of 4-point correlator on the double Keldysh contour with the parameterization.

    OTOC(t_1,t_2,t_3,t_4) = 1/N F(τ_1,τ_2,τ_3,τ_4), 
       with τ_1=β/2+ t_1,τ_2= t_2,τ_3=β/4+t_3,τ_4=-β/4+t_4

It is known that any 4-point function satisfies the Bethe-Salpeter equation even for the deformed contour <cit.>,

    F(τ_1,τ_2,τ_3,τ_4) = ∫_𝒞τ_5τ_6 K_lad(τ_1,τ_2,τ_5,τ_6)F(τ_5,τ_6,τ_3,τ_4),


where we ignore the inhomogeneous term F_0, which means the ladder with no rungs, since it's much smaller than the OTOC contribution. In Eq. (<ref>), 𝒞 refers to the deformed double Keldysh contour. 
Following the same technique used in the analysis of SYK, the imaginary variable τ_5 can be represented as t_5 + β/2 and t_5, with each single t_5 having two points in the backward and forward contour. These points are denoted by the blue dots 1, 2, 3, and 4, respectively. For τ_1= t + β/2, the contributions of blue points 3 and 4 cancel out, but the contribution from points 1 and 2 is non-vanishing, leading to the contribution Θ(t_1-t_5)( G_<(t_1,t_5)- G_>(t_1,t_5))=- G_R(t_1, t_5) <cit.>. Similarly, we find that τ_6 only contributes from blue points 3 and 4, leading to - G_R(t_2, t_6). This automatically leads to the simplification that τ_5 and τ_6 differ by an imaginary time β/2, which means the delta function in Eq. (<ref>) has no contribution. Additionally, F(τ_5,τ_6,τ_3,τ_4) will be invariant when considering points 1 and 2, since the causal relation between the blue points and τ_3 and τ_4 is fixed with certain imaginary time. After considering these structures, the real-time kernel reads.

    K_lad,R(t_1,t_2,t_5,t_6)    = -2G̃_R(t_1,t_3) G̃_R(t_2, t_5) 
            [f”(Δ̃_W(t_5,t_6))  ∂_t_5∂_t_6G̃_W(t_5,t_6) + 2 f'(Δ̃_W(t_5,t_6))  (∂_t_5∂_t_6) ]

where the Wightman Green's functions are defined by

    G_W(t) ≡⟨ x(t-β/2) x(0) ⟩,  Δ_W(t) ≡⟨(x(t-β/2)- x(0) )^2 ⟩

and their relation to the spectral function can be found in the Appendix <ref>. 
Finally, we take an ansatz for OTOC

    OTOC(t_1, t_2, 0, 0) = e^ϰ (t_1+t_2)/2ℱ(t_1-t_2),

and plug it into the Bethe-Salpeter equation 

    OTOC(t_1,t_2,0,0) = ∫_ℝt_5t_6 K_lad,R(t_1,t_2,t_5,t_6)OTOC(t_5,t_6,0,0).

These manipulations yield at last an eigenvalue problem for the OTOC,

    ℱ(ω) = -2∫ω'/2π|G_R(ω+ϰ/2)|^2 ( Λ_W,1(ω-ω') + Λ_W,2(ω-ω') ( ϰ^2/4+(ω')^2)) ℱ(ω'),

where the Wightman rung diagrams are

    Λ_W,1(ω)     = ∫t f”(Δ̃_W(t)) (-∂_t^2 Δ̃_W(t)) e^ω t
    Λ_W,2(ω)    = 2 ∫t  f'(Δ̃_W(t)) e^ω t.

Additionally, the analytical continuation to G_R(ω+ϰ/2) can be represented using spectral function, which reads 

    G_R(ω + ϰ/2) =1/-2∫ω'ρ_Δ(ω')/ω' - (ω+ ϰ/2).










We can diagonalize the RHS of Eq. (<ref>) to find the largest eigenvalue. As a sanity check of the numerics, we find that all the eigenvalues are real. We first discuss the numerical results for the high-temperature Lyapunov physics in Fig. <ref>. We find the Lyapunov exponent ϰ/(2π/β) exhibits the same power law for all the J^2 in Fig. <ref>(a). Via linear fitting on a log-log plot, we extract the power-law behavior ϰ/(2π/β) ∼β^0.8. In the high-temperature region, the Lyapunov exponent is far from the chaos bound, i.e., ϰ/(2π/β) ≪ 1. We also provide the eigenvalue h_i and the corresponding eigenfunction in the frequency domain ℱ_i(ω) by solving Eq. (<ref>) for a typical high-temperature setting, where β=0.06, J^2=0.6. As shown in Fig. <ref>(b), the largest eigenfunction corresponds to a positive and symmetric function centered at ω=0. Other eigenfunctions are negative in some regions, and the cancellation leads to smaller eigenvalues for those modes.

We show the numerical results for the low-temperature Lyapunov physics in Fig. <ref>. 
In Fig. <ref>(a), the Lyapunov exponent ϰ is seen to be set by 1/β, albeit with some mild temperature dependence. The smaller the J, the weaker the temperature dependence of ϰ/(2π/β) in this range of β. The results are all consistent with the chaos bound, ϰ≤ 2π/β. In Fig. <ref>(b), the Lyapunov exponent increases monotonically when J^2 increases. Due to numerical limitations, we are only able to probe the parameter region β < 800, J^2 ∼ O(1) with controlled numerical error. Hence, it is still open whether the system approaches maximal chaotic behavior at the largest Js and βs. But we observe no saturation with increasing J. Additionally, we also plot the eigenfunction ℱ_i(ω) associated with the eigenvalue in the low-temperature region. We find that the eigenfunctions ℱ_1, ⋯, ℱ_5 are qualitatively the same in the high-temperature and low-temperature regions. In the intermediate large ω region, ℱ_1(ω) shows exponential decaying behavior. However, the high-temperature eigenfunction in the small ω region will significantly deviate from the exponential decay, as illustrated in Fig. <ref>(b). But this phenomenon is not significant in Fig. <ref>(c).








§ DISCUSSION



In this paper, we presented a systematic study of the Magnetic Maze (MM), including its Euclidean dynamics, real-time dynamics, and Lyapunov physics. The structure of the large-N path integral framework, specifically the presence of certain time derivatives in the equations of motion, guarantees the absence of a glassy equilibrium state in this model. We formulated all the governing equations for general vector potential correlations, but we focused our analysis on the case where those correlations have a special power-law falloff, f(Δ) ∼ J^2/Δ. We found a rather rich theory in the low-energy regime, one with tunable dynamical exponent, which gives a new kind of controllable scale invariant quantum system.

We conclude by highlighting some open questions and directions for further work:



  
Low Energy EFT There is a lot we know about the Magnetic Maze at low temperature. We know the low-energy equations of motion have two emergent symmetries: reparameterization and rescaling. However, it is unclear what signatures these symmetries have in Lorentzian time, or what form any Goldstone bosons for these symmetries might take.



  
Maximal chaos? We did not observe any signature that the chaos exponent saturated with increasing J, therefore it is possible that the Magnetic Maze approaches maximal chaos in the limit of large J and low temperature. It would be very interesting to better understand this regime, possibly by making use of the Low Energy EFT.




  
Gravity dual? It is interesting to ask if the low-energy physics of the Magnetic Maze could have a gravity dual. This question is particularly interesting if it can be shown that the Magnetic Maze is maximally chaotic at low temperature in the strong coupling limit. The dynamical exponent z is a key feature to match, along with the pattern of symmetries.



  
Generalizations There are interesting generalizations of the Magnetic Maze, including the possibility of adding random electric fields (which could be equivalent to adding random magnetic fields to the p-spherical model), the possibility of considering externally imposed time-dependent fields (analogous to a random circuit model), and various many-body generalizations. On the latter point, we note that the variable dynamical exponent allows one to tune the particle interactions to be relevant or marginal even in high dimensions. We can also generalize to other observables, such as the spectral form factor which probes statistical properties of the energy spectrum.



  
Finite N It is also interesting to explore finite N corrections to our analysis. A systematic 1/N expansion would be useful as well as an analysis at small finite N, even N=2. The problem of random magnetic fields is well-studied in 2d, but the special power-law correlations we studied may not have been considered before. Finally, we note that cosmic magnetic fields can have a random character <cit.>, and our work might provide an alternative perspective on such systems via a 1/N expansion.




 
 We thank Xiao Chen, Ping Gao, Martin Sasieta
and Pengfei Zhang for useful discussions.
 T.-G. Z. acknowledges support from the Tsinghua Visiting Doctoral Students
 Foundation. MW acknowledges support from the Joint Quantum Institute. BGS acknowledges support from the U.S. Department of Energy through DE-SC0009986.









§ UNIFORM FIELD PROBLEM AND PHYSICAL SCALES




In this appendix, we review the basic physical scales in the magnetic maze. Throughout the paper, we absorb the charge into the vector potential and set k_B= 1, but we keep ħ explicit in this appendix. It is useful to approach this analysis by first discussing the case of a uniform magnetic field. This clarifies the physical meaning of the dimensionless parameters in the more general problem.



In N dimensions, a random spatially uniform magnetic field is an antisymmetric tensor B_ij which we take to be an N × N antisymmetric matrix with matrix elements given by Gaussian random variables with zero mean and variance equal to ^2/N. Antisymmetric matrices don't have real eigenvalues, but they do have eigenplanes in which the matrix takes the form [   0   λ; - λ   0 ] for some field strength λ. It can further be shown that the field strengths, which vary from plane to plane, have a semi-circle distribution at large N,

    ρ(λ) = N/2π^2√(4 ^2 - λ^2),

for 0 < λ < 2. When N is odd, there is one special direction that does not experience any magnetic field.



 §.§ Review of the 2d Uniform Field Problem


The N-dimensional problem thus breaks up into N/2 copies of the 2-dimensional problem. In the 2d classical problem with a uniform field λ, a particle with speed v undergoes uniform circular motion with angular frequency ω_c = λ/m and radius r_c = v/ω_c. In the corresponding 2d quantum problem, we also have the magnetic length, ℓ_B^2 ∼ħ/λ, and the characteristic energy of the Landau levels, ħω_c ∼ħ^2/m ℓ_B^2.

For a particle with a thermal velocity, v ∼√(T/m), the ratio of the classical orbit size to the magnetic length is

    r_c/ℓ_B∼√(T/ħω_c).

Hence, at high temperature the problem is effectively classical and at low temperature quantum effects can appear, as signaled by the size of the classical orbit approaching the magnetic length.

We are primarily interested in the squared displacement. At the classical level, the trajectory is

    x(t) = r_c cos( ω_c t ),        y(t) = r_c sin( ω_c t) .

The squared displacement, averaged over a thermal distribution for the velocity at temperature T = 1/β, is

    ⟨  (x(t)-x(0))^2 + (y(t)-y(0))^2 ⟩ = ⟨ v^2 ⟩4/ω_c^2sin^2 ω_c t/2 = 8/ m βω_c^2sin^2 ω_c t/2.

Quantum effects modify this result when ħω_c ≳ T.

Working in Landau gauge, the full quantum result is obtained by expressing the position operator in terms of the Landau level ladder operator, a, and momentum in the y-direction, p_y, 

    x = √(ħ/2 m ω_c) (a + a^†) + p_y/m ω_c.

The corresponding Heisenberg operator is 

    x(t) = √(ħ/2 m ω_c) (a e^- i ω_c t + a^† e^i ω_c t) + p_y/m ω_c.

From these, we get the thermally averaged correlator,

    ⟨ x(t)x(0) ⟩ = ħ/2 m ω_c( [n_B +1] e^-iω_c t + n_B e^i ω_c t) + ∑_k_yħ^2 k_y^2/m^2 ω_c^2,

where we note that the final term does not depend on time.

The average square displacement is 

    ⟨ (x(t) - x(0))^2 ⟩ = 2 ⟨ x(0)^2 ⟩ - ⟨ x(t) x(0) ⟩ - ⟨ x(0) x(t) ⟩ = 2 ħ/ m ω_c (2 n_B + 1)sin^2 ω_c t /2.

In the limit that T ≫ħω_c, we have n_B ∼ (βħω_c)^-1 and (<ref>) approaches half of (<ref>); the other half comes from including the average of (y(t)-y(0))^2.




 §.§ Uniform Field in N Dimensions




Returning now the N-dimensional problem, the key new effect is that we have a distribution of the field strengths λ. The N-dimensional averaged square displacement is thus

    Δ(t) = 1/N∫^2 _0 dλρ(λ) 4 ħ/λ(2 n_B(λ/m) + 1) sin^2 λ t/2m .

At zero temperature, with n_B=0, we can obtain a slow logarithmic increase with t. 

At non-zero temperature, the dominant effect comes from the classical contribution, which can be estimated as follows. We adjust the upper limit to take into account the freezing out of the classical motion for planes with a high cyclotron frequency. Then at large t, all but the smallest frequencies contribute an erratic amount proportional to sin^2  λ t/2m. Replacing this by its value averaged over one cycle, we estimate

    Δ(t) ∼1/N∫^m T/ħ_2m/t dλρ(λ) 4 m/βλ^2∼4 m/βπ∫^m T/ħ_2m/t d λ1/λ^2 = 2/βπ t + ⋯.

The first ∼ is replacing sin^2 λ t/2m with 1/2 and putting the lower cutoff where the erratic contributions to Δ begin. The second ∼ is approximating the density of λ by its central value, which is appropriate since the low field planes are the ones which contribute significantly to Δ. The modes with λ < 2m/t give a contribution of the same order. The final result shows that the particle moves diffusively at large t (cutoff only when t ∼ m/λ_min∼ N). It is also possible to perform the integral exactly using Bessel functions.



 §.§ Imaginary Time Dynamics Within Mean Field Theory


Strictly speaking, the mean field theory developed in section <ref> doesn't apply to the case of constant field. This is because it is impossible to choose a gauge where the A-field has statistical translation symmetry for constant fields. This, however, is a matter of fine print. One can get arbitrarily close to constant field by choose f(Δ) to be a positive function whose slope is a constant f'(Δ)=-^2 in the range of Δs we care about. In this limit, equation
 (<ref>) for the imaginary-time dynamics becomes quadratic in the Fourier transform G_ω, for all ω≠ 0. We have the imaginary-time equation
 
    (mω^2 +^2ω^2 G_ω)G_ω=1.

 
 It can be solved exactly, giving us 

    G_ω=2/mω^2+√(m^2ω^4+4^2ω^2).

If we are interested in large βs and small ωs, we can simplify further to

    G_ω=1/ |ω|

 

    Δ() ∼2/π^2 logβsinπ/β+𝒪(1).





 §.§ Spatially Varying Field


Now consider the problem of a spatially varying field. With our conventions, the vector potential has units of momentum and its correlations are translation invariant and given by

    A_i(x) A_j(0) = δ_ij f(x^2/N).

We will focus on the case in which

    f(u) = ħ^2 J^2/u + ℓ^2,

although other choices of f can be similarly analyzed. J is a dimensionless coupling which dials the strength of the local magnetic and ℓ is a characteristic length scale which sets the correlation length of the random vector potential. 

The magnetic field is obtained from A_i(x) as

    B_ij = ∂_i A_j - ∂_j A_i,

and its correlations can be straightforwardly obtained from those of f. At large N, the leading contribution to the magnetic field correlation is 

    ∂_i A_j(x) ∂_k A_l(0) = - δ_ikδ_j l/N f'(x^2/N) + ⋯,

which in our case is

    ∂_i A_j(x) ∂_k A_l(0) = δ_ikδ_j l/Nħ^2 J^2/(x^2/N+ℓ^2)^2 + ⋯.




To analyze the dimensionless parameters in the problem, it is convenient to choose ℓ as a unit of length. Combining this length with the mass m and ħ we get a characteristic energy

    E_0 = ħ^2/mℓ^2.

The characteristic scale of the magnetic field is

    ∼ħ J/ℓ^2.

As discussed just above, the magnetic field is a tensor with a characteristic semi-circle spectrum of field strengths (which now vary in space), and  sets the overall scale of the spectrum. The corresponding magnetic length is 

    ℓ_B/ℓ∼1/√(J)

and the magnetic energy is 

    ħω_c/E_0∼ J .


Putting this all together, the physics can be classified in terms of two dimensionless ratios,

    r_c/ℓ_B∼√(T/ħω_c)∼√(1/JT/E_0)

and

    ℓ_B/ℓ∼1/√(J).

In units where ℓ, m, and ħ are all set to one, we see that the problem is classical for T ≫ J and quantum for T ≪ J. In addition to setting the quantum-classical boundary, J also plays the role of the coupling, with larger J leading to larger deviations from the physics in no magnetic field. 

From the analysis in this appendix, we can form some intuition as to how the physics depends on T and J. When T ≫ J, all the local eigen-planes of B_ij are effectively classical, whereas when T ≪ J, only the low field eigen-planes remain classical. It is this smaller number of “active” planes which should be responsible for the low-temperature dynamics. Moreover, according to the uniform field analysis, (<ref>), the average square displacement after a thermal time, t ∼β, is 1/J. Hence, the larger J is, the slower the particle is moving. In the main text, we indeed find that larger J leads to slower motion in many senses, including by increasing the dynamical exponent at low temperature. 





§ ANALYTICAL SOLUTION AT J=0

To simplify the formulas, we set m=1. We also include the confining potential as a convenient regulator and only set it to zero at the end of the calculation.



 §.§ Imaginary Time

At limit J→ 0, the imaginary time Green's function in freqeuncy domain is

    G(ω_n) = 1/ω_n^2 + ω_0^2,

where ω_0^2 ≡ϵ. 

We can perform the Matsubara frequency summation to obtain

    G(τ)    = 1/β∑_n G(ω_n)e^-ω_n τ
        = e^-τω_0((e^2 τω_0+1) n_B(ω_0)+1)/2 ω_0.

From Δ(τ) =2G(0) - 2G(τ) we get

    Δ(τ)    = 2 e^-βω_0csch(βω_0/2) sinh(τω_0/2) sinh(1/2ω_0 (β -τ ))/ω_0
        = τ  (β -τ )/β+ τ  (τ -β )ω_0 + O(ω_0^2),

at which point it is safe to take ω_0 → 0.



 §.§ Real Time


We start with the basic formulas 

    G_R(ω) = 1/(ω +  0^+)^2 - ω_0^2,


    ρ_G(z) = 1/2ω_0δ(z-ω_0) + 1/-2ω_0δ(z+ω_0),


    G_>(t) = -/2ω_0 n_B(ω_0) ( e^-ω_0 t e^βω_0 - e^ω_0 t),


    G_<(t) = -/2ω_0 n_B(ω_0)( e^-ω_0 t - e^ω_0 t e^βω_0).

The zero-temperature version is actually given by

    G_R(t)    = ∫ω/2π1/2 ω _0( 1/-ω _0+(ω +i η )-1/ω _0+(ω +i η )) 
       = Θ(t) e^i t ω _0-e^-i t ω _0/2 ω _0


Notice that G_R(t) is actually well defined in the ω_0 → 0 limit, where G_R(t)=Θ(t) t. But for G_>,G_<,G_K,G_W, the limit is not well defined.
After subtracting the appropriate zero point, Δ(t) becomes

    Δ_>(t) = 2 csch(βω_0/2) sin(t ω_0/2) sin(1/2ω_0 (t+i β ))/ω_0

We find that Δ_>(t) > 0. If we take ϵ→ 0, we obtain the same formula via analytic continuation,

    Δ_>(t)    = t (t+i β )/β
    Δ_<(t)    = t (t-i β )/β.







 §.§ Wightman Green's Function

The Wightman Green's function is defined by

    G_W(t) ≡⟨ x(t-β/2) x(0) ⟩ .

Therefore we can derive it using non-equilibrium Green's function

    G_W(t)    =  G_>(t-β/2) 
       = ∫ω/2π 2π n_B(-ω) ρ_G(ω) e^-βω/2 e^-ω t
    ⇒ G_W(ω)    = -π1/sinh(βω/2)ρ_G(ω) 
       = 1/2sinh(βω/2)(G_R(ω) - G_R^†(ω) ) .


For Δ, we have Δ_W(ω)= -2 G_W(ω) for ω≠ 0. The ω=0 zero point needs special treatment due to the condition

    Δ_>(t-β/2,0) = 2(G_>(t-β/2,t-β/2)+G_>(0,0) - 2 G_>(t-β/2,0)).

With translation invariance, the correlator is

    Δ_>(t-β/2) = 2(G_>(0) - G_>(t-β/2)),

which leads to the Wightman Green's function at t=0

    Δ_W(0) = 2(G_>(0) - G_>(-β/2)) = (-i)Δ_τ(β/2) ,

Δ_τ is imaginary time Green's function. This is explicitly true for the J=0 case, where from the analytic continuation we get

    Δ_W(t) = Δ_>(t-β/2) = t^2 + β^2/4/β.




§ OTHER CHOICES OF VECTOR POTENTIAL CORRELATIONS



It is an interesting question to ask what's special about the 1/Δ form of the vector potential covariance. To answer this question, we can generalize the vector potential correlations to

    ⟨ A_i(x_1) A_j(x_2)⟩=δ_ij J^2 ℓ^2(ξ-1) / ( ℓ^2+|x_1-x_2|^2/N)^ξ,

with arbitrary parameter ξ >0. We can examine the behavior by performing imaginary time calculations for different ξ and trying to compare to the scaling ansatz in Eq. (<ref>). We try to extract an effective value of α and coefficient c_α from Δ(τ) using the ansatz 

    Δ(τ)=c_α( τ(β-τ)/β)^α

In Fig. <ref>, we demonstrate the result for different ξ. Only ξ=1 gives a result consistent with the IR value α_Theory obeying Eq. <ref>. When ξ>1, we expect the system  will flow to the free theory, which corresponds to α=1 and a J dependent coefficient c_α. This is consistent with the data in  Fig. <ref>.





§ SCHWINGER-KELDYSH PATH INTEGRAL

   

 §.§ Definition of Keldysh Rotation

   To establish the Schwinger-Keldysh formalism, we introduce the double contour with x_+, x_- fields. In order to construct other useful Green's function, we perform a Keldysh rotation on the original basis and introduce the new basis x_cl, x_q:
   
    [ x_cl;  x_q;      ] = M
       		[ x_+; x_-;     ]

   Here the transformation matrix is M= 1/√(2)[  1  1;  1 -1;    ].
   Correspondingly, the new Green's function or self-energy is defined on such basis, which reads
   
    [ G_K G_R; G_A   0;     ]≡
       		-<[ x_cl;  x_q;      ][ x_cl  x_q;      ]> 
       		= M [ G_++ G_+-; G_-+  G_–;      ] M^T 
       [   0 Σ_A; Σ_R Σ_K;     ]≡ (M^T)^-1[ (-1)^0+0Σ_++ (-1)^0+1Σ_+-; (-1)^1+0Σ_-+  (-1)^1+1Σ_–;              ] (M)^-1.

   The rationale behind defining Σ_R, Σ_A, and Σ_K is to maintain the structure of the Schwinger-Dyson equation (G_0^-1 - Σ) ∘ G = 𝕀 invariant after the Keldysh rotation.
   
   
   It's noticed that the statistical translational symmetry of the magnetic field still appears in our real-time formalism, even though δ(t_12) term doesn't appear in the Σ_≷. This is because the analog of the Schwinger-Dyson equation in real-time involves the retarded (or advanced) component of G and Σ after the Keldysh rotation. The retarded component Σ_R should, in principle, be represented as Σ_+-, Σ_+-, Σ_++, Σ_++, where the delta function term could enter the self-energy through Σ_++ and Σ_–. 
   Similar to the imaginary time formalism, we should also expect Σ_R(ω=0)=0 in real-time dynamics. This is consistent with the fluctuation dissipation theorem 2Σ_R(ω) = tanh(βω/2)Σ_K(ω), which is always true in the Keldysh formalism.

   

 §.§ Calculation of the Rotated Self-Energy

    Since the effective action Eq. (<ref>) involves the correlator Δ_ab, it leads to the unusual delta function term in the diagonal component of the self-energy in Eq. (<ref>). We need to prove that it still leads to the correct structure of self-energy. The structure of self-energy after and before Kelydsh rotation is
    
    [   0 Σ_A; Σ_R Σ_K;     ]=1/2[ Σ _++-Σ _+--Σ _-++Σ _– Σ _+++Σ _+--Σ _-+-Σ _–; Σ _++-Σ _+-+Σ _-+-Σ _– Σ _+++Σ _+-+Σ _-++Σ _–;                        ]

    To begin with, we show all central results in Tab. <ref>
    
    In the following subsection, we will derive the first four results below. The lesser and greater self-energy are by definition. The fluctuation-dissipation theorem for the self-energy can be a consistent check for our special Schwinger-Keldysh structure, which requires Σ_R(ω=0) - Σ_A(ω=0)=0. 
    
    

  §.§.§ Constraint on Σ

    Before we discuss each component, we can obtain some useful relations. First, we recall the definition of time-ordered and anti-time-ordered correlator
    
    Δ_++(t_1,t_2)    = Θ(t_12)Δ_-+(t_1,t_2) + Θ(-t_12)Δ_+-(t_1,t_2) 
    Δ_–(t_1,t_2)    = Θ(t_12)Δ_+-(t_1,t_2) + Θ(-t_12)Δ_-+(t_1,t_2).

    Then we calculate the summation of time-ordered and anti-time-ordered self-energy
    
    Σ_+++Σ_–
    
        		=    -∂_t_1∂_t_2f(Δ_++(t_1,t_2)) -  f'(Δ_++(t_1,t_2))∂_t_1∂_t_2Δ_++(t_1,t_2)  -∂_t_1∂_t_2f(Δ_–(t_1,t_2)) -  f'(Δ_–(t_1,t_2))∂_t_1∂_t_2Δ_–(t_1,t_2) 
       + δ(t_12) ∫t_3(  f'(Δ_++(t_1,t_3))∂_t_1∂_t_3Δ_++(t_1,t_3) +  f'(Δ_–(t_1,t_3))∂_t_1∂_t_3Δ_–(t_1,t_3) 
                    -   f'(Δ_+-(t_1,t_3))∂_t_1∂_t_3Δ_+-(t_1,t_3) -  f'(Δ_-+(t_1,t_3))∂_t_1∂_t_3Δ_-+(t_1,t_3) )

    After we expand the time-ordered and anti-time-ordered correlator using Eq. (<ref>), we find the delta function term is exactly canceled, and the regular self-energy terms can be simplified as Σ_+-+Σ_-+. Therefore it immediately leads to the first result:
    
    
    1/2(Σ_+++Σ_–-Σ_+--Σ_-+)= 0.


    

  §.§.§ Keldysh Self-Energy Σ_K

    After using the constrain relation Eq. (<ref>), Keldysh Green's function reads
    
    Σ_K≡1/2(Σ_+++Σ_–+Σ_+-+Σ_-+) = Σ_+-+Σ_-+.


	

  §.§.§ Retarded Self-Energy Σ_R

    To deal with the retarded self-energy component, we further assume time translation invariance and simplify the formula Σ_++-Σ_–, which reads
    
    Σ_++-Σ_–
    
                =   ∂_t_1^2f(Δ_++(t_12)) +  f'(Δ_++(t_12))∂_t_1^2Δ_++(t_12)  -∂_t_1^2f(Δ_–(t_12)) -  f'(Δ_–(t_12))∂_t_1^2Δ_–(t_12) 
       + δ(t_12) ∫t_3( - f'(Δ_++(t_13))∂_t_3^2Δ_++(t_13) +  f'(Δ_–(t_13))∂_t_3^2Δ_–(t_13)  ) 
    
                =   (Σ̅_++ - Σ̅_–) + δΣ_1 +  δΣ_2.

    The result can be separated into a regular part Σ̅_++ - Σ̅_– and the δΣ_1,δΣ_2 part corresponding to the delta function term. 
    
    Firstly we consider δΣ_2, which is contributed by the explicit delta function term in the self-energy Σ_++,Σ_–, where
    
    δΣ_2(t_12) =δ(t_12)∫t_3(-  f'(Δ_++(t_13))(∂_t_3^2)Δ_++(t_13) +  f'(Δ_–(t_13))(∂_t_3^2)Δ_–(t_13)  )

    It can be further simplified using the following tricks
    
    δΣ_2(t_12)
                   =δ(t_12)∫t_3( - f'(Δ_++(-t_3))∂_t_3^2Δ_++(-t_3) +  f'(Δ_–(-t_3))∂_t_3^2Δ_–(-t_3)  ) 
       =δ(t_12)∫t_3( -Θ(-t_3) f'(Δ_-+(-t_3))∂_t_3^2Δ_-+(-t_3) - Θ(t_3) f'(Δ_+-(-t_3))∂_t_3^2Δ_+-(-t_3)  
                      +Θ(-t_3) f'(Δ_+-(-t_3))∂_t_3^2Δ_+-(-t_3) + Θ(t_3) f'(Δ_-+(-t_3))∂_t_3^2Δ_-+(-t_3) )  
       =-2 δ(t_12)∫t_3Θ(t_3)(   f'(Δ_-+(t_3))∂_t_3^2Δ_-+(t_3) -   f'(Δ_+-(t_3))∂_t_3^2Δ_+-(t_3) )

    In the first equation we absorb variable t_1 into the t_3 due to the integration is carried out from -∞→∞. In the second equation, we expand the time-order and anti-time-ordered correlator. In the third equation, we want to deal with the minus sign on the t_3 variable. For the terms with Θ(-t_3), we can change the integration variable t_3 → - t_3. Besides, for the terms with Θ(t_3), we can remove the minus sign in Δ by using the condition Δ_ab(t)=Δ_ba(-t), which origins from the fact that displacement operator x is a Hermitian operator.

    
    Secondly, the δΣ_1 denotes the delta function origins from the kinks in the time-ordered or anti-time-ordered correlator Δ_++,Δ_–. This contribution only appears in the combination of ∂_t_1^2f(Δ_++(t_12))-∂_t_1^2f(Δ_–(t_12)).
    By expanding Δ_++,Δ_– using Eq. (<ref>), we can obtain the kink contribution
    
    ∂_t_1^2  f(Δ_++(t_12))  - ∂_t_1^2  f(Δ_–(t_12)) 
    
                =   ∂_t_1^2 ( Θ(t_12) f(Δ_-+(t_12))+Θ(-t_12) f(Δ_+-(t_12)) ) - ∂_t_1^2 ( Θ(t_12) f(Δ_+-(t_12))+Θ(-t_12) f(Δ_-+(t_12)) ) 
    
                =     2δ(t_12)( ∂_t_12 f(Δ_-+(t_12))- ∂_t_12 f(Δ_+-(t_12)) ) + (t_12) (∂_t_1^2  f(Δ_-+(t_12))  - ∂_t_1^2  f(Δ_+-(t_12)) )

    From the second line to the third line, we have considered both the kink contribution using ∂_t Θ(t) = δ(t), and the regular part with a sign function. Therefore, the δΣ_1(t_12) can be defined as simplified as
    
    δΣ_1(t_12)    = 2δ(t_12)( ∂_t_12 f(Δ_-+(t_12))- ∂_t_12 f(Δ_+-(t_12)) ) 
       = -2δ(t_12) ∫t_3(Θ(t_3)∂_t_3^2  f(Δ_-+(t_3)) - Θ(t_3)∂_t_3^2  f(Δ_+-(t_3)) )

    Finally, the regular part is
    
    Σ̅_++-Σ̅_–   = (t_12) (∂_t_1^2  f(Δ_-+(t_12))  - ∂_t_1^2  f(Δ_+-(t_12)) +   f'(Δ_-+(t_12))∂_t_1^2Δ_-+(t_12) -   f'(Δ_+-(t_12))∂_t_1^2Δ_+-(t_12)) 
       =(t_12) (Σ_-+(t_12) -Σ_+-(t_12)),

    which immediately leads to the regular part contribution to the retarded self-energy
    
    Σ̅_R(t_12) ≡1/2(Σ̅_++ - Σ̅_– +Σ _-+ -Σ _+-) = Θ(t_12)(Σ_-+(t_12) -Σ_+-(t_12)).

    This is the structure that we expect for the regular retarded self-energy term. Remarkably, the delta function contribution to the retarded self-energy can also be greatly simplified to
    
    1/2(δΣ_1 +δΣ_2)= - δ(t_12) ∫t_3Θ(t_3) (Σ_-+-Σ_+-) = -δ(t_12) Σ̅_R(ω=0).

  This leads to the final result for retarded self-energy:
    
    Σ_R(t) = 1/2(Σ _++-Σ _+-+Σ _-+-Σ _–)    = Σ̅_R(t) + 1/2(δΣ_1 +δΣ_2) 
                =Σ̅_R(t) - δ(t) Σ̅_R(ω=0).

    
    

  §.§.§ Advanced Self-Energy Σ_A

    The advanced self-energy can be easily obtained in a similar way:
    
    Σ̅_A(t)    ≡Θ(t)(Σ_+-(t) -Σ_-+(t)) 
    Σ_A(t)    =  Σ̅_A(t) - δ(t) Σ̅_A(ω=0).














§ NUMERICAL METHODS FOR THE REAL-TIME SADDLE POINT EQUATIONS



As we discuss in the main text section <ref>, due to the statistical spatial translational invariance of the system, the real-time saddle point cannot be solved self-consistently using a simple mixing method since it's numerically unstable. Instead, we choose to use a generalized gradient descent protocol with a mask function.

From equation Eq. (<ref>), we can insert the saddle point equation and represent all self-energies via Δ,

    S/N=    - 1/2log 
            (G^-1)  -  ∫ t_1  t_2  ( 1/2Σ∘ G  + 1/2 I_0(G) ) 
    
            =     - 1/2log 
            (G^-1)  -  ∫ t_1  t_2  ( 1/2( G_0^-1∘ G - 𝕀δ(t_12) )  + 1/2 I_0(G) ) 
    
            =    - 1/2log 
            (-2Δ(ω)^-1)  -   (const+ 1/2∫ω/2π G_0(ω)^-1Δ(ω)/-2  + ∫ t_1  t_2 1/2 I_0(G) ).

The last step is to decompose ∫ t_1  t_2 I_0(G) and take its derivative:

    δ/δΔ(ω)∫ t_1  t_2  ( 1/2 I_0(G) ) 
    
             =   ∫ t δ G(t)/δΔ(ω)δ/δ G(t)∫ t_1  t_2 ( [∂_t_1∂_t_2 f(Δ_ab(t_12)) ] G(t_12) ) 
    
             =   ∫ t 1/∫ t' δΔ(t')/δ G(t)e^-ω t' (-Σ(t)) 
    
             =   ∫ t 1/∫ t' (-2δ(t'-t) - 2δ(t)) e^-ω t' (-Σ(t)) 
    
             =   ∫ t 1/-2 e^-ω t- 2δ(t)δ(ω) (-Σ(t)).


Since we're working in the ω≠ 0 sector, the δ(ω) term can be ignored. Summing all results we arrive at

    S/Δ(ω)/N  
            =   1/2Δ^-1(ω)  +  ( -1/2(  Δ_0^-1) + 1/2Σ(ω)/-2) 
    
            =   /4( (-2)Δ^-1(ω) - (G_0^-1 - Σ)  ).

If we only consider the retarded component, we will find the gradient to be zero is the same as the Schwinger-Dyson equation(<ref>). In the numerics, we choose the update rule to be

    Δ(ω)_new = Δ(ω)_old + ζ/4( (-2)Δ^-1_R(ω) - (G_0,R^-1(ω) - Σ_R(ω))  ) Δ_R(ω)^2.

Since the gradient can diverge near small ω, we introduce an extra Δ_R(ω)^2 factor as a mask function. This numerical trick stabilizes the iteration process. We consider convergence to be achieved when ||Δ_R,new(ω)-Δ_R,old(ω)|| < 10^-6.













JHEP


	
