Disconnecting a Traversable Wormhole: Universal Quench Dynamics in Random Spin Models

Understanding strongly interacting quantum matter and quantum gravity are both important open issues in theoretical physics, and the holographic duality between quantum field theory and gravity theory nicely brings these two topics together. Nevertheless, direct connections between gravity physics and experimental observations in quantum matter are still rare. Here we utilize the gravity physics picture to understand quench dynamics experimentally observed in a class of random spin models realized in several different quantum systems, where the dynamics of magnetization are measured after the external polarization field is suddenly turned off. Two universal features of the magnetization dynamics, namely, a slow decay described by a stretched exponential function and an oscillatory behavior, are respectively found in different parameter regimes across different systems. This work addresses the issues of generic conditions under which these two universal features can occur, and we find that a natural answer to this question emerges in the gravity picture. By the holographic duality bridged by a model proposed by Maldacena and Qi, the quench dynamics after suddenly turning off the external polarization field is mapped to disconnecting an eternal traversable wormhole. Our studies show that insight from gravity physics can help unifying different experiments in quantum systems.

The holographic duality between quantum field theory at boundary and gravity theory in bulk has shed new insights in understanding both quantum matter and gravity [1,2]. For example, on the gravity side, the wormhole in local Einstein gravity is not traversable, and recently mechanisms that can render a wormhole traversable have been studied by coupling quantum fields on boundaries [2][3][4]. On the quantum matter side, insights from gravity theory can help us understand strongly interacting quantum systems, such as non-Fermi liquids [2]. Nevertheless, so far only a few specifically designed quantum models, such as the Sachdev-Ye-Kitaev (SYK) model [1,6,7,9], have been explicitly shown to possess holographic duality to gravity theory. Such a model is seemingly different from realistic systems in quantum materials and is also hard to be realized by quantum simulations. Therefore, so far little connection between sights from gravity physics and realistic experimental observations in quantum matter has been established through holographic duality, except for few known examples such as bounds for viscosity and transport coefficients [10][11][12][13][14].
In this work we focus on a typical random spin model written aŝ This model has recently been realized by cold molecules in optical lattices [15,16], and NV centers [17], fermions in harmonic traps [18], Rydberg atoms [19], high spin atoms in optical lattices [20] and by solid state NMR [21,22]. J ij is random all-to-all spin interaction coefficients, which originates either from dipolar interaction or from the spin exchanging interaction, and the randomness comes from the random locations of spin carriers in these systems. We can denote J ij =J/N + δJ ij , wherē J/N is the averaged value of all J ij and N is the number of spins. δJ 2 ij = 4J 2 /N and J denotes the typical strength of random spin interactions. ∆ is the anisotropy between the longitudinal and the transverse spin-spin interactions. h is an external polarization field. Among these systems, /J varies from a few microseconds to a few milliseconds, and the values ofJ/J and ∆ are also different.
Taking the advantages of the controllability of these systems, these realizations enable experimental studies of far-from-equilibrium quantum dynamics, and the quench dynamics in these random spin models have also been reported in a number of recent experiments [15][16][17][18][19][20]. In the quench dynamics, spins are initially polarized in the transverse direction, say, alongx, by h-term in Eq. 1, and then the dynamics of the total magnetization alonĝ x are experimentally measured after h is suddenly turned off. In some of these experiments where the magnetization decay time due to the single-spin effect is shorter, a spin echo is applied to eliminate the single-spin effect [15][16][17]20], which ensures that the observed features are all due to spin-spin interactions described by Eq. (1). Two features of the quench spin dynamics have been observed in the time scale from a few to a few tens of /J. They take place in different parameter regimes, but both The correspondence between the Penrose diagram in the gravity theory (left) and the lesser Green's function G < (t), as well as the spectral function ρ(ω) (right) in the quantum theory, are shown for the BH phase (B) and the TW phase (C), respectively. In the Penrose diagrams, the colored regimes are the physical space in the geometry. Two yellow dashed lines represent the physical boundaries of two AdS2 spaces, denoted by the left (L) and the right (R) ones, where the quantum fields are defined. The black dotted lines with an arrow represent the null geodesic. The horizontal axis is the real-time axis.
are observed in different systems and across different parameter regimes. These two features are listed below.
• Slow decay of the total magnetization has been found when the magnetization decays to zero at long time, and the time dependence of the total magnetization is identified as a stretched exponential function. This feature has been observed in experiments on NV centers [17], Rydberg atom [19] and high spin atoms in optical lattices [20], • When the total magnetization saturates at finite value at long time, the total magnetization oscillates in time. This feature has been observed in experiments using cold molecules in optical lattices [15,16], and fermions in harmonic traps [18].
Here we will address the question of the general conditions under which these two features will occur in these random spin models. This work is to show that a scenario emerged from the dual gravity theory unifies these experiments. This work establishes a connection between the quench dynamics in these random spin models and the dynamics after suddenly turning off the coupling in a traversable wormhole. The connection is bridged by a model coupling two SYK systems proposed by Maldacena and Qi, as shown in Fig. 1(a). The Maldacena-Qi model is dual to a gravity theory hosting an eternal traversable wormhole [2], and a clear correspondence can be established between the behaviors of the spectral function in different phases of the Maldacena-Qi model and different space-time geometry in the gravity side with or without a traversable wormhole. We will show that the random spin models are closely related to the Maldacena-Qi model, and the spectral functions of these two models share similar behaviors. Hence, by studying the evolution of the spectral functions during the quench dynamics, we can establish a gravity interpretation of the quench dynamics.

Maldacena-Qi Model and its Gravity Dual.
Since the Maldacena-Qi model is a coupled SYK model, we need to first briefly introduce the SYK model. A single SYK model is written aŝ whereψ i (i = 1, . . . , N ) are N Majorana fermion operators, and J ijkl are random Gaussian variables with zero mean and the variance J 2 ijkl = 3!J 2 /N 3 . The SYK model is exactly solvable in the large-N limit. It can be shown explicitly that the low-energy physics of the SYK model has an emergent conformal symmetry, which can be dual to the Jackiw-Teitelboim gravity in the AdS 2 geometry with a black hole.
Various extensions of the SYK models have been considered, especially coupled SYK models [2, 4, 23-31, 33, 34]. Here we focus on the version proposed by Maldacena and Qi [2]. The Maldacena-Qi model is written aŝ whereψ L i andψ R i denote the left and the right Majorana fermions. If µ = 0, there is no coupling between two SYK models, and on the gravity side, two AdS 2 spaces are also independent. Therefore, the geodesic starting from the boundary of one AdS 2 space can never reach the boundary of the other AdS 2 space. The µ-term adds coupling between two boundaries of the AdS 2 spaces which distorts the boundary geometries. At finite temperature, when the distortion is weak, two boundaries still cannot be connected by a geodesic and this is referred to as the "black hole" (BH) phase. However, when the distortion is strong enough, the geodesic starting from the boundary of one AdS 2 space can reach the boundary of the other AdS 2 space, which renders these two black holes into a traversable wormhole. This is referred to as the "traversable wormhole" (TW) phase. The difference between these two geometries is represented by two different Penrose diagrams in Fig. 1

(B) and (C).
In the presence of the holographic duality, the behavior of correlation functions are closely related to the connectivity of the space-time geometry through the standard holographic dictionary. Let us consider the lesser Green's function G < aa (t) = i l ψa l (t)ψ a l (0) /N , where a, a = L, R. Roughly speaking, this Green's function can be physically interpreted as follows: first, a signal is created at t = 0 by an operator at boundary a , and then the signal propagates following the geodesic, and finally, the signal is detected at time t on boundary a. In the BH geometry, since there is no connection between two boundaries, it is natural to have G < aa (t) = 0 for a = a . For a = a , G < aa (t) decreases as t increases, because as one can see from the Penrose diagram in Fig. 1(B), the geodesic, denoted by the black dotted line, moves away from the boundary. However, the situation is different for the TW geometry. In the TW geometry, as also shown by the Penrose diagram in Fig. 1(C), the geodesic starting from the left boundary can reach the right boundary and can be further bounced back and forth between two boundaries for an eternal traversable wormhole [2]. Thus all G < aa (t) are always finite for a = a and a = a . At a time when the geodesic arrives at one of the boundaries, the corresponding Green's func-tion reaches a maximum. Therefore, G < aa (t) shows an oscillatory behavior in time. These two different behaviors of G < aa (t) are also illustrated in Fig. 1(B) and (C). Similar behaviors can also be found in the quantity i l [ψ a l (0),ψ a l (t)] /N = −G > aa (t) − G < aa (t), which more directly reflects whether and how the information gets transferred from one side to the other.
Knowing the Green's function G < aa (t), one can obtain G < aa (ω) by performing the Fourier transformation, and then derive the spectral function ρ aa (ω) = G < aa (ω)/2πin F (ω) by the fluctuation-dissipation theorem, where n F (ω) is the Fermi-Dirac distribution. In the BH phase, ρ aa (ω) shows a single broad peak near ω ∼ 0 for a = a , and is significantly smaller for a = a . In the TW phase, ρ aa (ω) displays multiple narrow peaks. ρ aa (ω) behaves similarly for a = a and a = a , except for their different parities in ω. This difference is also shown in Fig. 1(B) and (C). Hence, we have established correspondence between the behavior of the spectral function and the presence of a traversable wormhole. Equipped with this correspondence, we are now ready to move on to discuss the random spin models. First of all, the random spin model can be written in terms of fermion operatorsĉ i,s with spins s =↑, ↓, by rewriting the spin operators asŜ α i = 1 2ĉ † i,s (σ α ) ss ĉ i,s , where α = x, y, z and σ α denotes the corresponding Pauli matrices. Then, we define the imaginary-time Green's function of fermions G(τ ) as a 2 × 2 matrix, with the matrix elements defined as We first consider the case withJ = 0. The Fourier transformation of G(τ ) defines G(iω n ), and by the Schwinger-Dyson equation, we have with the Matsubara frequency ω n = (2n + 1)π/β and β = 1/(k B T ) being the inverse temperature, where Σ(iω n ) is the self-energy in the Matsubara frequency domain. With the large-N and the melon diagram approximations, the self-energy can be obtained as (see Methods for details) For the Maldacena-Qi model, we can similarly define a 2×2 matrix G(τ ) with G aa (τ ) = i T τψ a i (τ )ψ a i (0) /N , with a, a = L, R, and G(iω n ) in the Matsubara frequency domain also obeys the Schwinger-Dyson equation as It has been shown that the Green's function obeys the self-consistency equation [2] Σ aa (τ ) = J 2 G aa (τ ) 3 .
We first consider the case ∆ = 0, where direct evaluation of Eq. (6) yields wheres = s. First of all, to compare Eq. (5) with Eq. (7), we notice that the spin model has rotational symmetry alongẑ, and therefore one can replace σ x in Eq. (5) by σ y . Then note that the system is invariant when rotating π alongŷ:ĉ i → iσ yĉi , which gives G ↑↑ (τ ) = G ↓↓ (τ ) and G ↑↓ (τ ) = −G ↓↑ (τ ), and the system also has the particle-hole symmetryĉ i →ĉ † i that corresponds to (Ŝ These symmetry relations simplify Eq. (9) into Σ ss (τ ) = J 2 G ss (τ ) 3 , which takes the same form as Eq. (8). In this way, Eq. (5) and Eq. (9) for the random spin model become completely identical to Eq. (7) and Eq. (8) for the Maldacena-Qi model. Two spin components s, s =↑, ↓ play the role of a, a = L, R, and the Zeeman field h plays the role of the coupling term µ, which is responsible for the wormhole being traversable. The equivalence between these two models also holds approximately for non-zero but small ∆, where a non-zero small ∆ provides an extra channel to couple two sides, as one can see from the low energy effective action (see Supplementary Materials for details).
Next, we discuss how to deal with a non-zeroJ. Up to Here we take ∆ = −0.73, which is realized in the Rydberg atom experiment [19], and we setJ = J. a constant, the term associated withJ can be cast intō We implement the molecular field approximation by de- Then, Eq. (10) can be written as For equilibrium, M α needs to be determined selfconsistently. For dynamics, M α evolves in time. In this work, since we initially polarize spins alongx, we focus on the situation only M x = 0 and M y = M z = 0. Hence, this term can be combined with the h-term by replacing The quenching process we will consider is instantaneously turning off h from non-zero to zero, therefore, here we will first discuss the equilibrium phase diagram for both finite and zero h, respectively. The model parameters contain J, ∆, h andJ, as well as temperature T . Here we take J as the energy unit. When h is finite, it polarizes spins, and M x is non-zero, which further normalizes h tol through the effect ofJ. When the self-consistency is reached, there is a one-to-one correspondence between h and h tol , and only h tol enters the equation that determines the spectral functions. Hence, in Fig. 2(A), we fix h tol and plot the phase diagram in terms of the other two dimensionless parameters T /J and ∆. The phase diagram is obtained by numerically solving the self-consistency equations discussed above. Here the "BH" and "TW" are distinguished by their different behaviors of spectral functions, as we illustrate in Fig.  1 (B) and (C) and show in Fig. 2(C). Consistent with the Maldacena-Qi model [2], the random spin model with small ∆ displays the BH phase at high temperature and the TW phase at low temperature (see Supplementary Material for the discussion of the phase transition).
When h = 0, the system possesses a spin rotation symmetry of (Ŝ If this symmetry is respected, then M x is always zero and h tol is also zero. Nevertheless, the ground state can have non-zero M x due to the spontaneous symmetry breaking, which leads to non-zero h tol . Both a non-zeroJ and a non- zero ∆ can lead the system into a TW phase [4,33]. However, as shown in Fig. 2(B), the TW phase is significantly smaller compared with the h = 0 case shown in Fig. 2

(A).
Quench Dynamics. Now we discuss the two universal features observed in the quench experiments. As we have seen, the h-term is mapped to the µ-term in the Maldacena-Qi model, which couples two boundaries and makes the wormhole traversable. Hence, in the gravity picture, the quench dynamics of suddenly turning off h corresponds to turning off the coupling term, which is responsible for making a wormhole traversable. Below we will calculate the quench dynamics numerically by the Kadanoff-Baym formula, which gives how the spectral function evolves in time (see Method). By the correspondence between the spectral function and the space-time geometry discussed above, we can establish a gravity picture of the quench dynamics.
Firstly, we discuss the slow dynamics described by a stretched exponential. Initially, the spectral function is a typical TW-type one, as shown in Fig. 3(C1). We show in Fig. 3(A) that the dynamics contain two stages. In the first stage when tJ < t * J, with t * J ≈ 2 in the case shown in Fig. 3(A), M x decays quite fast. As shown in Fig. 3(C2), the spectral function still retains TWphase-like in this stage. It is known that the distribution function F(ω) at equilibrium should be tanh (βω/2), as in the initial case shown in Fig. 3(D1). As one can see in Fig. 3(D2), at the end of the first stage, the system strongly deviates from equilibrium.
In the second stage, with tJ ranging from t * J to a few tens, the dynamics is much slower. Note that if M x obeys a stretched exponential e −C(t/t0) η , it should manifest as a straight line in the log(| log M x |) − log t plot, and the slope determines η. We show this in Fig. 3(B), and the fitting yields η < 1, meaning that the dynamics is slower than usual exponential decay. In Fig. 3(C3), we show that the multiple peaks in the spectral function gradually merge into a single peak, and meanwhile, ρ ss with s = s is gradually suppressed. At the longtime limit, as shown in Fig. 3(C4), the spectral function becomes BH-phase-like, where ρ ss shows a single broad peak around ω ∼ 0 for s = s and is vanishing small for s = s . Hence, in the gravity picture, this process corresponds to the wormhole gradually closing. From Fig.  3(D3),(D4), we also see that the distribution functions fast reach quasi-equilibrium in this stage, and by fitting the distribution function with tanh(βω/2), we can obtain an effective temperature 1/β for the saturated state. By comparing Fig. 3(D4) with (D1), one can see that the temperature for the final state is close to the temperature of the initial state. By examining various different parameters (see Supplementary Material for more examples), we identify that the stretched exponential decay occurs generically when a traversable wormhole is disconnected and becomes two decoupled black holes, and the temperatures of the initial TW and the final BH phases are close.
Secondly, we discuss the oscillatory behavior. In Fig.  4 we compare two cases which only differ by the parameter βJ but show two different behaviors. For the one with a smaller βJ, M x monotonically decays to zero as discussed above. For the one with a larger βJ, M x saturates to a finite value in the long-time limit, and oscillates in time before saturation. Here we also want to find a generic condition to determine which behavior should take place, and it turns out that the answer is also quite clear in the gravity picture. As we show in Fig. 4(B) and (C), these two cases have similar initial states in the TW phase. However, their long-time saturation states are very different. ρ(ω) at long-time exhibits typical behavior of the BH phase for the one with smaller βJ, but exhibits typical behavior of the TW phase for the one with larger βJ. Note that in Fig. 2(B), even with h = 0, there still exists TW phase in the equilibrium phase diagram, which means that turning off h does not always disconnect the wormhole. If the long-time saturation phase is still retained in the TW phase, there will exhibit an oscillatory behavior of magnetization in time. This oscillatory behavior can also be found by using the effective action with gravity dual (see Supplementary Material). By examining different parameters (see Supplementary Material), we also find that whether the monotonical decay or oscillatory behavior occurs in the quench dynamics essentially depends on whether the final state is the BH phase or the TW phase.
Summary and Outlook. In summary, we show that the quench dynamics of randomly interacting quantum spins after turning off an external field can be understood in the dual gravity picture as turning off a coupling field for making the wormhole traversable. If this process disconnects the traversable wormhole and finally yields two decoupled black holes, the magnetization displays a slow decay described by the stretched exponential function toward zero magnetization. If this process does not disconnect the traversable wormhole, due to the existence of the residual couplings between the boundary fields, the magnetization saturates to a non-zero value at long time and oscillates around the saturation value. This gravity picture conceptually unifies these two universal phenomena observed in the random spin model realized in different physical systems.
The stratagem of our study is reminiscent of the concept of the"fixed point" in the renormalization group theory, where different microscopic models can share the same low-energy description given by a fixed point action. Here we use the Maldacena-Qi model, whose gravity dual has been shown explicitly, as the "holographic fixed point", and we use the holographic fixed point to define correspondence between quantum properties and the gravity properties. We further share this correspondence with different microscopic models in the neighborhood of the Maldacena-Qi model with similar quantum properties, and these models are much closer to reality. Therefore, the success of this approach can inspire more applications on realistic quantum models with insights from the gravity theory.

Methods.
Derivation of the Self-Consistent Equations. The original model given by Eq. (1) can be written by fermion representation withŜ α i = ss ĉ † i,s (σ α ) ss ĉ i,s /2 and the constrain sĉ † i,sĉ i,s = 1. We now promote this model by adding an additional M indices aŝ wherê Here m i = 1, 2, ...M and T γ are the generators of the SU (M ) group that satisfies The external field h only couples to the SU (2) part. The constrain is also promoted as s,mĉ † i,s,mĉ i,s,m = M . Firstly, we take the imaginary time approach in the large-N and large-M limits. The constrain is satisfied automatically due to the particle-hole symmetry, and the self-energy can be obtained by the melon diagrams as [6] Here we have defined By taking M → 1, we obtain Eq. (4) and Eq. (6) in the maintext with Σ(τ ) = Σ (1) (τ ) and G(τ ) = G (1) (τ ). Secondly, we also employ the real-time approach and we consider the M = 1 case. We begin with the greater and lesser Green's functions in real-time defined as as well as the retarded and advanced Green's function defined as where Θ (t) is the Heaviside step function and we have defined (16) In the thermal equilibrium, all Green's functions are only functions of t 12 due to the time-translational symmetry. To solve the real-time Green's functions self-consistently, we introduce the spectral function as which implies ρ(ω) = −ImG R (ω)/π. Other Green's functions are determined by the fluctuation-dissipation theorem as where n F (ω) is the Fermi-Dirac distribution function.
To derive the self-consistent equation for the spectral function, we apply the Fourier transformation on the Eq. (6) Σ(τ ) → Σ(iω n ) and perform the analytical continuation iω n → ω + i0 after the Matsubara frequency summation for the internal dummy frequency [1]. We Iteratively solving these equations gives the spectral function ρ(ω) and real-time Green's functions.

The Kadanoff-Baym Equation.
To study the real-time dynamics, we utilize the Kadanoff-Baym equation on the Keldysh contour [35], which describes the real-time evolution of G ≷ . Assuming the melon diagram approximation Eq. (6) and applying the Langreth rules [36] on the Schwinger-Dyson equation, we find that [37]: ).
(20) Here we have taken the contribution ofJ into account by The real-time self-energies in Eq. (20) are given by is given by the equilibrium solution, which serves as the initial conditions for the time dynamics. Evolving G ≷ χ (t 1 , t 2 ) using Eq. (20)- (22) gives the quantum dynamics.
To analyze the quench dynamics, we study the spectral function and the quantum distribution function at given time t. We firstly separate the relative time t r by defining the temporal Green's functionG ι (t, t r ) as Here ι ∈ {>, <, R, A} and we perform the Fourier transformation for the relative time. The temperal spectral function is then defined as Similar as Eq. (24), we could also defineΣ ι . At thermalequilibrium, there is no time dependence for bothG and Σ, and the fluctuation-dissipation theorem implies the relation [35] tanh Then for the non-equilibrium case, motivated by Eq. (26), we define the quantum distribution function at time t as Here we report the results from the exact diagonalization calculation of a finite-size spin system, which shows both the slow dynamics and the oscillatory dynamics in different parameter regimes. The results are consistent with those reported in the main text based on the large-N expansion and assuming the melon diagrams. This consistency provides supports to our approximations used in the main text. Nevertheless, we should also emphasize that, because the finite-size exact diagonalization and large-N theory are two different approximations, we do not expect a quantitative agreement between these two calculations.
Firstly, slow dynamics is shown in the Fig. S1(A) and (B). In the same parameter region as in our large-N theory discussed in the main text, we observe the stretched exponential decay behavior of M x . In the log(| log M x |) − log t coordinates, we find the decay exponent η < 1.
Secondly, we find some hints for the oscillatory behavior as shown in Fig. S1(C). Here we should first note that the symmetry breaking behaviors can hardly be achieved in such a small size spin system with N = 10 by the exact diagonalization method. In fact, we do not find oscillatory behavior in the exact diagonalization calculation when we quench the external field h to zero. This is actually consistent with our conclusion that the oscillatory behavior is associated with the final state in the TW phase, where M x is non-zero. Hence, to show the oscillatory behavior in the exact diagonalization calculations, we quench h to a very small value, say, h = 0.02, instead of h = 0. The results are shown in Fig. S1(C). It shows that for small βJ, at a long time the magnetization saturates to a constant, with very small variation. As βJ increases, the long time dynamics exhibits a larger and larger variation around its long-time averaged value, although the variation is not a perfect oscillation due to the finite-size effect.

A. Equilibrium solution
Here we will discuss the low-energy effective action for the random spin model in the TW phase. We can identify two parts of contributions in addition to the Schwarzian action. The first part corresponds to the transverse magnetic field h tol , the effect of which has already been considered in the Maldacena-Qi's work [2]. The second part represents the effect from a non-zero ∆, with the assumption that |∆| 1. Here we show how these two quantities affect the energy gap E gap in the TW phase, and the equilibrium magnetization M x .
To begin with, it is known that under the melon-diagram approximation, the effective action in the imaginary time can be written as (S1) Here we define the imaginary-time Green's function as: Here s, s =↑ or ↓ are the up or down spin indices respectively. From now on, we focus on the small h tol /J, where the dynamics can be described by the evolution of reparametrization modes. After imposing the spin rotational symmetries and the particle-hole symmetries discussed in the main text, direct calculation shows To the order of O(∆), Eq. (S3) can be truncated as When ∆ = 0, this effective action is identical to the effective action of the Maldacena-Qi model. After considering the last term in Eq. (S4) and rescaling the time variable, the effective action for the reparametrization modes becomes .
(S5) whereτ = J √ 2αs τ is the dimensionless imaginary coordinates time. For the SYK 4 model, α s ≈ 0.007 is a positive numerical factor in the action, and the scaling dimension D = 1/4. It's known that the Schwarzian action itself S3 describes the low energy fluctuations of the SYK model, and here we use {f, τ } to denote the Schwarzian derivative. Variable t l (τ ) and t r (τ ) in the low energy action are the time reparametrization modes on the left and right boundaries respectively. We have neglected other modes including the charge fluctuation modes, since they are not excited in our cases. The rescaled coupling constants arẽ These two parameters are related to the magnetic field h tol and the anisotropic paramter ∆. In the equilibrium case, the classical solution assumes the ansatz t l (τ ) = t l (τ ) = t τ , where t is a constant. The∆ term could be simplified as: . (S7) We have employed the ansatz and introduced the new coordinates x 1 = 1 2 (τ 1 +τ 2 ), x 2 = (τ 1 −τ 2 ). The dimensionless valueβ = βJ Adding up all the terms, the action reads It is also known that the classical saddle point dI dt = 0 determines the equilibrium solution as After the asymptotic expansion and assuming t 1, we can obtain equilibrium solution as The ∆ = 0 case reproduces the result in Ref. [2]. It is known that, in the gravity interpretation, the energy gap in the TW phase is E gap = J √ 2αs t D. Then after using (S10), we reveal different contributions of h tol and ∆ to the E gap term.
Second, with the value of t , the magnetization M x can also be calculated by the Green's function: plugging Eq.(S10) into the above relation between M x and t , the effects of h tol and ∆ are shown clearly.

B. Quench Dynamics with Molecular Field Correction
In the subsection A, we have studied the equilibrium properties by the imaginary time approach. In this subsection, we will study the real time evolution of this random spin model, starting from the TW phase. Therefore, we perform the analytically continuation on the action (S5) by substituting the coordinate timeτ → it + 0 (τ → it + 0) and reparametrization modes t l/r → it l/r . Here, for simplicity, we consider the case with ∆ = 0. Subsequently, Eq. (S5) becomes which is the same as the one in the Maldacena-Qi model. The action has SL(2) gauge symmetry, which means we should require the Noether charge vanishes Q 0 = 2N e −ϕ [−e 2ϕ − ϕ +hDe 2Dϕ ] = 0 with t l (t) = t r (t) = e ϕ(t) [2]. This determines the corresponding quench dynamics, during which the magnetization at coordinate time t can be expressed as (S13) Considering the non-zeroJ effect in the dynamical process, we can substitute h tol as h tol (t) = −JM x (t). Consequently, the equation of motion becomes with the proper initial conditions for this ordinary differential equation given by ϕ(0) = 2 log 2(2π) 1/4 √ α s M x (0) and ϕ (0) = 0. By numerically solving this differential equation, we find that forJ < 0 and h > 0, and in a parameter range ofJ/J, M x first decreases and then oscillates around certain non-zero value as time increases.

S3. PHASE TRANSITION
Here we introduce more details of two methods that we used to determine the phase diagram. The first one is the imaginary-time approach with which we can directly evaluate the thermodynamical quantities, such as the free energy f = − 1 βN log Z = 1 βN I (see Eq. (S1)) and the magnetization M x . The second one is the real-time approach, from which we can obtain the spectral function ρ ss (ω)(s, s =↑, ↓) and the spectral ratio we introduced as r ≡ ρ ↑↑ (0)/max ω ρ ↑↑ (ω).
With the imaginary-time approach, we start the self-consistent calculation from the high-temperature region, initialized with a proper high-temperature initial guess of G(τ ). When the self-consistent solution is reached at a certain temperature, we decrease the temperature and continue the self-consistent calculation for a slightly lower temperature, using the current self-consistent G(τ ) as an initial guess. In this way, we can obtain a series of M x and f evolving from the high-temperature side, denoted by M x H and f H . The calculation can be carried out similarly by evolving from the low-temperature side, which can also obtain M x and f denoted by M x L and f L . After obtaining the thermodynamics data by self-consistent calculation initiated from high-and low-temperature sides, we can determine the two phases and their boundary, as well as the order of phase transitions. For the imaginary time approach, we find both the crossover behavior (Fig. S2(A)) and the first-order transition behavior (Fig. S2(B,  C)). At relatively large |∆| region |∆| −0.5, both M x and f change continuously as lowering the temperature, and solutions initialized from high-temperature or low-temperature are identical. At relatively small |∆| region, say |∆| 0.5, thermodynamical properties behave as typical first-order transitions. Both M x and f show hysteresis effect for two different initializations from high-temperature and low-temperature sides, and a discontinuous jump exists at a certain temperature. Consequently, we can exactly determine the location of the first order transition temperature T * by the point f H (T * ) − f L (T * ) = 0. With the real-time approach, we can obtain the spectral function, and by definition, the spectral ratio r ∈ [0, 1]. r = 0 means a gap at ω = 0, describing the spectral function of the TW phase as discussed in the main text. r = 1 means that the spectral function is peaked at ω = 0, describing the spectral function of the BH phase. Here we can also find different transition types. For the first order transition region, r jumps discontinuously between 0 and 1. And in the crossover region, the gap closes gradually and r changes continuously.

S4. QUENCH DYNAMICS
In this section, we provide more examples to support the two main conclusions on the quench dynamics proposed in the main text. Fig. S3 shows the results of the quench dynamics, including the time dependence of the magnetization and the spectral functions in the initial and final states before and after the quench dynamics. In addition, we summarize the examples in the Table. S1. Clearly, in all these examples, we can see that • When initially r < 1 and finally r = 1, we observe a stretched exponential decay with exponent η < 1.
• When initially r < 1 and finally r is also < 1, we observe that the magnetization saturates at finite value at long time and oscillatory behavior at long time.  . For each given parameter, the results are described with a (decay exponent η/O, initial βJ → final βJ, initial r → final r). Where the symbol O stands for the oscillatory behavior around a non-zero saturation value. Here we recall that spectral ratio r ∈ [0, 1] and r → 0 means the TW phase and r → 1 means the BH phase. J/J 0 1 ∆ \ βJ