Inhomogeneous quasi-adiabatic driving of quantum critical dynamics in weakly disordered spin chains

We introduce an inhomogeneous protocol to drive a weakly disordered quantum spin chain quasi-adiabatically across a quantum phase transition and minimize the residual energy of the final state. The number of spins that simultaneously reach the critical point is controlled by the length scale in which the magnetic field is modulated, introducing an effective size that favors adiabatic dynamics. The dependence of the residual energy on this length scale and the velocity at which the magnetic field sweeps out the chain is shown to be nonmonotonic. We determine the conditions for an optimal suppression of the residual energy of the final state and show that inhomogeneous driving can outperform conventional adiabatic schemes based on homogeneous control fields by several orders of magnitude.


Introduction
Techniques to control or assist adiabatic dynamics are of broad interest in quantum technologies, including quantum simulation and quantum computation [1,2]. The breakdown of adiabatic dynamics in quantum critical systems is conveniently described using the Kibble-Zurek mechanism (KZM) [3][4][5][6][7]. This is the paradigmatic theory to account for the nonequilibrium dynamics across a continuous quantum phase transition. It exploits the divergence of the relaxation time in the neighborhood of the critical point in combination with scaling theory to predict the density of excitations in the final state that results from crossing the transition at a finite rate. The KZM can also be used to estimate the mean energy of the final state, known as residual energy, when measured with respect to the corresponding ground state energy [8]. Both the density of excitations and the residual energy are shown to scale as a universal power law of the quench rate, where the power-law exponent is fixed by the correlation length and dynamic critical exponents, ν and z, as well as the dimensionality of the system.
Strategies that have been developed to mimic adiabatic dynamics in quantum critical systems, in the absence of disorder, often boil down to finding ways out of the KZM. This is challenging as the KZM is broadly applicable and it holds even in strongly-coupled systems [9,10]. Yet, a variety of protocols have been put forward. Prominent examples include the driving of finite many-body systems with a nonzero energy gap [11], coupling the system of interest to a thermal bath [12], and engineering the time variation of the driving field using scaling theory [13,14] or optimal control [15,16]. Further approaches include the design of counterdiabatic fields to implement shortcuts to adiabaticity [18][19][20][21][22], or the modulation of multiple control parameters in time [23], see [24] for a review. All these strategies are particularly crucial-and should be additionally scrutinized-in the presence of noisy control fields, that preclude adiabatic dynamics [17].
Recently, it has been shown that the KZM should be modified when spatial or temporal inhomogeneities affect the critical behavior of classical and quantum systems. In the absence of disorder, the residual energy dependence on the quench rate can be enhanced in classical inhomogeneous systems [25][26][27][28][29][30][31]. This is the case when the critical point is first reached at a local front that subsequently spreads throughout the system, completing the crossing of the phase transition. A number of recent experiments are consistent with the fact that the interplay between the velocity of the critical front and the speed of information paves the way to suppressing defect formation [7]. The role of causality has also been established in quantum spin chains with no disorder [34][35][36].
By contrast, the development of techniques to approach the adiabatic limit in disordered systems is much more limited. This is however a pressing issue for boosting the performance of quantum annealers that encode combinatorial optimization problems and thus are inherently disordered [1]. It has been shown that, for disordered Ising spin chains, the KZM predictions are severely modified and the residual energy of the state upon completion of the driving exhibits only a weak dependence on the quench rate [37,38]: it no longer follows a power-law and vanishes only with the inverse of the logarithm of the quench rate. One main outstanding challenge, that we address in this work, is to explore the effect of inhomogeneous driving across a quantum critical point in the presence of disorder. The main goal is to spatially coordinate symmetry breaking events among neighboring regions by finding the appropriate degree of inhomogeneity and the speed of critical front to reduce the number of topological defects.
In this work, we introduce an inhomogeneous driving strategy for a weakly disordered Ising spin chain by a transverse magnetic field with a smooth step-like spatial profile that sweeps out the chain from side to side, as illustrated in figure 1. The length scale in which the field is modulated controls the number of spins that simultaneously cross the critical point. We study the dependence on the shape (slope) of the profile and the velocity at which it is moved throughout the chain to minimize the residual energy of the final state, effectively approaching adiabatic dynamics. We show that our local driving protocol can outperform conventional quantum annealing schedules based on homogeneous fields, reducing the relative residual energy by several orders of magnitude.
The paper is organized as follows: in section 2 we introduce the weakly disordered transverse Ising model and briefly review previous relevant results. In section 3 we use the adiabatic theorem to show the existence of a threshold velocity of the critical front below which the evolution is adiabatic. By an analysis in the spirit of the inhomogeneous KZM, we show that this velocity has a universal character for smooth fronts. Subsequently, in section 4 we present simulations of full dynamics to find the inhomogeneous protocols that optimize the residual energy. While for fast driving a homogeneous control field proves advantageous, in slower schemes for which the critical front does not exceed the threshold velocity the inhomogeneous protocol with a smooth front extended over several spins turns out to be optimal, effectively reaching the adiabatic limit. We close with a discussion and conclusions in section 5.

Model
We consider a chain of N spins described by the random Ising Hamiltonian with quenched disorder (constant in time) in the nearest-neighbors couplings J n . We set   = = 1 in subsequent calculations, that is equivalent to use   as a unit of time. Realizations of disorder are drawn from a uniform distribution The equilibrium properties of the model are well understood as it is solvable using the strong disorder renormalization group approach [39], see [40] for a review. In particular, the critical point satisfies ( ) (| |) = g J log log c n for uniform g n = g c ( ¼ 0.9558 ), and belongs to the universality class of the infiniterandomess fixed point.
We consider a driving protocol where at an initial time t i the system is prepared in the ground state for the magnetic field ( ) = g t g n i i deeply in the paramagnetic phase (in the simulations we use g i = 3) and then driven at a finite rate to the final value ( ) = g t g n f f , deeply in the ferromagnetic phase. We set g f =0 for which the Hamiltonian in equation (1) can be considered classical, with a non-vanishing energy gap (notice that  J 0.5 n which we refer to as weak disorder) and is outside of the Griffiths phase surrounding the critical point. Under homogeneous driving and in the absence of disorder the nonequilibrium dynamics is well described by the KZM [3][4][5][6][7]. The mechanism exploits the divergence of the equilibrium relaxation time [ ] | | t e t e = n z 0 , known as critical slowing down, as a function of the dimensionless distance to the critical point In the proximity of g c the modulation of the magnetic field can be linearized in the form ( ) where t Q sets the quench time. The critical point is reached at t=0 around which the relaxation time diverges. As the system is driven through the phase transition, the evolution can be roughly split in three sequential stages where the dynamics is approximately adiabatic, frozen and adiabatic again. The time scale in which the system leaves the frozen stage to evolve adiabatically in the broken-symmetry side of the transition is known as the freeze-out timeˆ( ) t t = . In the Ising model without disorder, n = = z 1, and the KZM prediction for the density of excitations reads [41][42][43], as corroborated by the exact solution [41]. Similar power-laws can be derived for other observables using scaling theory [8,44,45]. This is the case for the residual mean energy defined as the difference between the mean energy of the system following the completion of the protocol and the corresponding ground state energy, e.g.,= á ñ - 1 effectively diverges as the system approaches the critical point (  ¥ z as e  0) [39]. As a result, the scaling of the density of excitations is no longer described by a power-law and a more careful analysis based on KZM predicts t d 1 ln Q 2 in a slow transition (for large t Q ) [37,38]. The dependence on the quench time becomes then much weaker than in the absence of disorder. We note that a similar result can also arise in a clean model as a result of a particular decoherence mechanism [46].
The existence of these logarithmic scaling laws signifies that one is forced to consider exponentially long evolution times to reduce the residual energy of the final state. As an alternative to a global homogeneous modulation of the magnetic field, we introduce an inhomogeneous protocol where the linear front interpolating between g i and g f travels through the chain with velocity v and gradually drives the system from the paramagnetic to the ferromagnetic phase, by sweeping out the chain from end to end. We denote by = n vt f the position of the spin for which the magnetic field equals the arithmetic mean of g i and g f . The slope α sets the effective number of spins driven at a given instant. The resulting protocol is illustrated in figure 1 and interpolates between the following two limiting cases: (i) homogeneous driving, which is recovered in the limit of a  0 and  ¥ v while keeping a t =v Q 1 fixed, and (ii) driving of one spin at a time, when a » 1 [47].
In absence of disorder, the KZM can be extended to account for an inhomogeneous scenario [27][28][29][30][31]. The central prediction is the existence of threshold velocity v t that determines the relevance of the driving scheme. When the velocity of the front widely surpasses this threshold value,  v v t , the effect of the local modulation of the control magnetic field is negligible and the nonadiabatic critical dynamics resembles that under homogeneous driving, well described by the standard KZM. By contrast, in the case  v v t the length scale in which the front is smoothed out becomes relevant and determines the number of spins that simultaneously experience criticality. The smooth front of the inhomogeneous field opens up an energy gap, that allows for adiabatic evolution.
For a smooth front with  a 1, KZM predicts that the penetration depth across the critical point, i.e., the size of the critical region separating the phases to both sides of the inhomogeneous front, varies as [31][32][33] ( ) When the gap at the critical point vanishes polynomially with the system size this leads to the opening of instantaneous gap that scales asˆx a D~-n n+ i i z z 1 . By combining the characteristic time and length scales in the problem one can then estimate the threshold velocity as In the Ising model without disorder z=1 and v t is a constant independent of α (  a 1). The analytical solution [34] shows that when J n = 1 in equation (1), v t =2 and equals the sound velocity at the critical point. What is more, the transition between the adiabatic regime for < v v t and the effectively homogeneous regime for > v v t is actually sharp.
The presence of disorder changes the universality class of the model and the assumptions leading to equation (5) do not longer hold. Naively setting  ¥ z in that equation leads to a vanishing threshold velocity and would indicate that inhomogeneous driving does not favor adiabatic dynamics in disordered systems. In this article we show this not to be the case. Indeed, the analysis presented in the next section predicts that the threshold velocity v t acquires a finite value, that admits a universal form for small α. This paves the way to implement adiabatic dynamics by inhomogeneous driving.
All simulations shown below are done using the Jordan-Wigner transformation that maps the Hamiltonian in equation (1) onto the system of free-fermions where it can be solved numerically in a standard way. For details, we refer the reader to appendix B.

Threshold velocity at the adiabatic limit
We next provide a quantitative prediction of the threshold velocity under quasi-adiabatic dynamics when diabatic transitions occur within the manifold spanned by the ground and the first-excited states. The formation of excitations is proportional to the mixing matrix element between these two states, . The instantaneous energy gap D n f can be parameterized both as a function of the front position n f  and the time of evolution t, are the energies of the instantaneous ground state | ñ t 0, and the first excited state | ñ t 1, of Hamiltonian (1), respectively. Since the parity operatorˆs =  = P i N n z 1 is conserved during time evolution, we consider only the subspace with the same parity as the initial ground state.
According to the adiabatic theorem, the evolution follows the instantaneous ground state with high fidelity provided that which allow us to estimate the value of threshold sweep velocity v t below which the dynamics of the quench is effectively adiabatic as The factor of 4 in the definition of v t above is introduced so that ( ) v n t f matches the exact known value in the case without disorder, when v t = 2, see appendix A for a detail discussion of this case.
The matrix element reads This analysis is restricted to quasi-adiabatic dynamics governed by adiabatic following of the ground-state and | -ñ 0 | ñ 1 transitions. As a result, it is expected to fail when transitions to higher excited states are dominant, e.g., for large α (» 1), when the inhomogeneous front approaches a step function.
In figure 2(a) we show the instantaneous gap during the evolution for a single realization of disorder and fixed value of the slope a = 1 32. While the gap fluctuates as the front travels through the chain, it remains finite. We define D min as the minimal gap for a given realization of disorder. This definition involves averaging over a finite-length chain (N = 512 in figures 2 and 3) and in principle depends on N. However, since only a fixed fraction of the system is being driven at given instant, we expect the dependence to be weak for finite a > 0. In that case, the critical front can be thought of as probing different local realizations of disorder where the effective size of the systemx i (set by α) is finite. As a result, fluctuations of the instantaneous gap are limited and there is a negligible probability of having a gap smaller than a given threshold. This can be qualitatively seen in figure 3(a) where we consider many realizations of disorder and show different quantiles of D min , as we shall discuss further at the end of this section where we provide a quantitative scaling argument. By contrast, in the homogeneous case (a = 0) the typical energy gap at the critical point is expected to vanish as ( ) -C N exp [38,39,48,49], where C is a nonuniversal constant.
In figure 2(b) we show the mixing term ( ) W n f for the same realization of the couplings, and by combining it with the gap, from equation (9) we estimate the local value of threshold velocity of the front ( ) v n t f below which the evolution is expected to stay adiabatic. We define v t as the minimum of ( ) v n t f for a given realization; see figure 2(c), where its value is of order unity for a = 1 32. Since ( ) v n t f is widely fluctuating one could envision a fine-tuned inhomogeneous driving protocol where the velocity of the front is adjusted with respect to the local value of the threshold velocity. However, the design of the corresponding driving protocol g n (t) would require quite a specific knowledge of the system. Here, to explore the broad applicability and universality of the inhomogeneous scheme, we keep the front velocity v fixed along the evolution.
The results obtained from sampling over many realizations of disorder are summarized in figure 3. It shows the dependence of the average value of the minimum gap D min , the largest value of mixing matrix denoted by W max as well as the threshold velocity v t as a function of the smoothness α of the inhomogeneous magnetic-field front for a disordered Ising chain. The average is taken over 1000 realizations from which statistics is built to determine the quantiles corresponding to the hardest (95%) and easiest (5%) cases, displayed by dotted lines. All quantities increase monotonically as a function of α in the range of values considered, with the exception of the threshold velocity that levels off for large values of α.
For smooth fronts corresponding to small values of α, a power-law fit yields a Wm ax 1.1 similar to the Ising case with no disorder where  a W . The gap D min , however, disappears faster than polynomially in that limit which results in a monotonic dependence of v t on α (without disorder a Dm in 1 2 ). The maximum value of velocity is obtained for a value of α close to but lower than unity, when the inhomogeneous front extends over a few sites. We also notice that the optimal α is smaller than 1 (i.e. the limit of driving one spin at a time).
For a more quantitative scaling prediction, we consider the distribution of the instantaneous gap ( ) D n f , where the front is traveling inside the chain as in the middle part of figure 2(a). In doing so, we disregard the configurations in which the front is entering or leaving the chain and the gap is large. We denote by ( ) a D P , In order to numerically verify equation (12), for each realization of disorder we sample values of D n f for a couple of hundreds equidistant instances of n f  corresponding to the front traveling inside the chain; see figure 2(a). We collect the values of the gap obtained that way for a couple of thousands realizations and from the histogram we extract the probability distribution ( ) q a D P , as a function of α. Depending on α the statistics is built from >10 6 points from >2000 realizations of disorder. The distributions for different  a 1 collapse onto each other corroborating our scaling prediction; as seen for  a 1 64 in figure 4. In addition, this distribution coincides with the one obtained for homogeneous system at the critical point (g n = g c ), when the scaling variable  figure 4, the distribution differs from the universal one. We expect that this happens due to the effective finite-size effect, where the size of the critical region is so small that non-universal contributions are still relevant. While the analytical strong disorder renormalization group approach of [39,40] cannot be directly used in the presence of inhomogeneous (position correlated) field, we expect that a number of initial decimation RG steps would be necessary to approach universal fixed point trajectory.
The typical value of instantaneous gap scales as ( ) a D~--C exp n 1 3 f in the limit a  0. Here, we are mostly interested in the minimal gap D min , which would be determined by the behavior of the tail of ( ) q a D P , corresponding to small energies. The derivation in [49] suggests that we can expect this tail to be Gaussian, which is indeed consistent with the data in figure 4 (see footnote 4). As the front travels though the chain, it samples the distribution ( ) q a D P , in a continuous way; see figure 2(a). We can estimate the dependence of the minimal gap D min on the system size N by assuming that, to leading order, N instances are drawn from the distribution ( ) q a D P , and determining the probability distribution for the minimal value. With a Gaussian tail this means that any fixed quantile of the minimal (global) gap (e.g. plotted in figure 3(a) for N=512), vanishes slower than a polynomial in N, making this dependence a subleading correction. The leading contribution related to the system size is the time needed for the front to travel though the chain with fixed velocity which is proportional to N (for fixed α).
Summarizing, the above analysis suggests the existence of a finite threshold velocity v t for non-zero α and a maximum at α near unity, when the front extends over few sites. However, given our analysis in terms of the the instantaneous ground state and the first excited state, the values of v t could in principle be overestimated. This is especially true for large α close to unity, when the first excited state is not well separated from the rest of the spectrum. This could be addressed by using adiabatic theorem taking into account all excited states, see e.g. [50,51]. We however take a different approach in the next section, namely, the numerically-exact simulation of the full dynamics.

Optimization of the protocol and residual energy
Given the existence of a finite threshold velocity v t discussed in the previous section, we next explore the possibility of implementing adiabatic dynamics under inhomogeneous driving. In particular, we focus on the minimization of the residual energy Q of the final state.  5 Both in homogeneous and inhomogeneous case we use ( )   D = + 2 1 2 (see [38] and equation (B2)) which is the minimal gap in the sector with given parity. Interestingly, we obtain the collapse even though in the inhomogeneous case  = 0 1 , which is related to part of the chain being in the ferromagnetic phase which could break the symmetry. We note, that the scaling results in [48,49] were obtained for  D = 2 1 .
We note that the total time required for the inhomogeneous protocol to be completed reads where the first term corresponds to the time needed for the middle of the front to travel through the chain and the second term accounts for additional time needed for the magnetic field to reach the final value for all the spins. In the homogeneous limit, the second term in (13) dominates and the total time reduces to | |t g g i f Q (in this limit  ¥ v as a  0). By contrast, in the strongly inhomogeneous limit (when ) the first term dominates and the second one constitutes a small over-head. Here, to compare different protocols we fix the value of T and choose the velocity v according to equation (13) for a given α and N.
A direct analysis of the performance of the inhomogeneous driving scheme is shown in figure 5 that displays the dependence of Q on both v and α for a fixed value of the ramp time T. The slope α in this plot interpolates between the nearly homogeneous transition and the limit of a steep front where only one spin is driven at a time. We observe that for short time scales the homogeneous driving is optimal. However, for longer ramps, when the velocity of the front is small enough, there is a sharp minimum for intermediate values of α where the dynamics is effectively adiabatic. The value of the threshold velocity that marks the appearance of that minimum is consistent with v t , obtained in figure 3 for intermediate and small values of α. Notice, however, that the actual minimum of Q is reached for a slightly smaller value of α than suggested by that figure.
This efficiency in suppressing excitations is shown to be fairly insensitive to the hardness of the problem as quantified by the quantiles considered. Approximately the same landscape is observed for quantiles with = q 50% and = q 99%, where q marks the percentage of realizations which have smaller residual energy. Further, the optimal value of α weakly depends on the quantile. In particular, for = T 10 4 , a » 1 16 is optimal for harder cases with = q 99% and a » 1 32 is optimal for = q 50% 6 . Notice however that in both cases the value of the residual energy is almost that of an adiabatic transition, with the fidelity larger than 0.9999 in both cases. In this limit the actual position of Q might also depend on the additional smoothing of the critical front in equation (3) that is nonanalytic at the point between the piecewise linear and constant sections of the front [52].
Finally, we identify scenarios for the supremacy of the inhomogeneous driving scheme over its homogeneous counterpart in figure 6. In particular, we plot in figure 6(a) the time needed to reach a given quality of the solution, as quantified by the inverse of residual energy density. We compare the two schemes, when the optimal value of α for the inhomogeneous scheme is found from the landscape studies as in figure 5. The homogeneous transition (or a  0) is shown to be optimal for short time scales in figure 6(a). However, for long time scales the residual energy after such a quench is expected to scale only as ( ) Q N T 1 log 2 [37,38], making it unpractical to reach the adiabatic limit. The inhomogeneous driving approaches this limit for nonzero α for long enough time scales such that the velocity of the inhomogeneous front is reduced below the threshold value. This is further confirmed in figure 6(b) where we compare the performance of homogeneous and inhomogeneous protocols with fixed a = 1 32 for a time-scale of = T 10 4 , for different system sizes. As Figure 5. Optimization of the inhomogeneous driving protocol. Landscape plots of residual energy for various time scales and ( ) a v for N=512. We show the = q 50% and = q 99% quantiles (q is the percentage of realizations with smaller residual energy) obtained from simulation of 1000 realizations (500 for = T 10 4 ). While for short time scales the best residual energy is obtained for homogeneous driving, for longer time scales the sharp minimum appears for intermediate values   a 1 16 1 32. The residual energy for optimal smoothness,  a 0.03 at T=10000 timescale, is five orders of magnitude smaller than both the standard annealing strategy with a = 0 and the inhomogeneous scheme based on flipping one spin at a time, shown for  a = 2 2 2.82 (see footnote 5).
the system size is reduced for given T, the velocity of the front is proportionally smaller, which is the reason behind the weak increase of residual energy with growing N in the inhomogeneous case.

Conclusions
In summary, we have analyzed the driving of weakly-disordered spin chains with a time-dependent magnetic field. Under spatially homogeneous driving, the minimization of the residual energy in the final state is essentially constrained by the adiabatic theorem. Long-evolution time are then required to minimize excitations. As an alternative scheme, we have proposed the use of an inhomogeneous magnetic field that sweeps out the system at a well-controlled velocity. In this scenario the spatial modulation of the magnetic field introduces an effective system size that favors adiabatic dynamics. The dependence of the residual energy of the final state on the latter and the sweeping velocity is not monotonic. Upon optimization with respect to these two parameters, we have identified the supremacy of inhomogeneous driving over homogeneous schemes in reducing the residual energy of the final state. In this article we have focused on the case of weak-disorder where the couplings J n are nonzero. We shall discuss the case of strong disorder with possibly vanishing couplings J n in a subsequent article [53].
Our results can prove useful in the design of novel quantum annealing protocols with inhomogeneous control fields on disordered spin systems with the potential to outperform conventional schemes [54] and might be applicable to open systems [55]. Inhomogeneous schedules with controllable local magnetic field can be implemented on the next generations of quantum annealers that are currently under development by D-Wave System and the Google Quantum AI laboratory. Our approach might be applied for finding higher quality of solutions for constrained optimization problems over standard adiabatic quantum computation (with homogeneous fields), as the corresponding embedded problems on annealing hardware Hamiltonians would involve finding low-energy states of disordered spin glasses on low-dimensional lattices. Figure 6. Supremacy of optimal inhomogeneous driving over homogeneous schemes. Comparison of the homogeneous and best inhomogeneous protocol for different ramp times and system sizes. (a) 99% quantile of residual energy. The symbol + denotes the best inhomogeneous strategy, where the optimal α is extracted from landscape plots in figure 5. For long-enough timescales it is advantageous to use larger a » 1 16, while for short ramps homogeneous driving a  0 is better. Circles show the residual energy for a homogeneous quench in the same ramp time. (b) Comparison of homogeneous (red) and inhomogeneous (blue) strategies with fixed a = 1 32 (T = 100 00). Solid lines indicate the median of 1000 realizations for 4 system sizes while the dotted lines mark 1% and 99% quantiles, respectively.
For convenience we introduceˆ †   = H a Ha, where  a is a column vector consisting of operators a n and H is á N N 2 2 matrix. We separate the matrix describing full Hamiltonian = + H H H J g into parts corresponding respectively to the coupling and magnetic field, which are both block diagonal, , and the quasiparticles energies  n are arranged in ascending order. We note that some care is needed for a system with degenerated ground state,  = 0 1 (within numerical precision), in which case one has to take the proper linear combination of eigenvectors of H to eigenvalue 0, to ensure that O 0 is orthogonal.
The parity operatorˆ( with the initial condition ( ) = = O t O 0 0 which corresponds to starting in the ground state of the initial Hamiltonian. We numerically solve this differential equations by employing 4th order time dependent Suzuki-Trotter decomposition [56], which is symplectic and allows to greatly speed up the calculations: we split the Hamiltonian matrix H into parts corresponding to H J and H g , that are block-diagonal in the original basis. This facilitates the propagation at intermediate steps involving terms of the form