Engineering non-equilibrium quantum phase transitions via causally gapped Hamiltonians

We introduce a phenomenological theory for many-body control of critical phenomena by engineering causally-induced gaps for quantum Hamiltonian systems. The core mechanisms are controlling information flow within and/or between clusters that are created near a quantum critical point. To this end, we construct inhomogeneous quantum phase transitions via designing spatio-temporal quantum fluctuations. We show how non-equilibrium evolution of disordered quantum systems can create new effective correlation length scales and effective dynamical critical exponents. In particular, we construct a class of causally-induced non-adiabatic quantum annealing transitions for strongly disordered quantum Ising chains leading to exponential suppression of topological defects beyond standard Kibble-Zurek predictions. Using exact numerical techniques for 1D quantum Hamiltonian systems, we demonstrate that our approach exponentially outperform adiabatic quantum computing. Using Strong-Disorder Renormalization Group (SDRG), we demonstrate the universality of inhomogeneous quantum critical dynamics and exhibit the causal zones reconstructions during SDRG flow. We derive a scaling relation for minimal causal gaps showing they narrow more slowly than any polynomial with increasing size of system, in contrast to stretched exponential scaling in standard adiabatic evolution. Furthermore, we demonstrate similar scaling behaviour for random cluster-Ising Hamiltonians with higher order interactions.

Controlling non-equilibrium dynamics of quantum many-body systems is of one the main challenges in condensed matter physics and quantum control. Such complex quantum systems have very rich parameter space and unusual dynamical properties that makes them very hard to simulate and control as they are driven through critical regions 1 . The main difficulties arise from the fact that these systems generally contain high degree of disorders and effectively low dimensions such that they are not prone to exact analytical treatment or meanfield approximations. In principle their dynamics can be mapped to the dynamics of spin-glass systems that are driven/quenched by external control fields and could experience various first and second-order quantum phase transitions and Griffiths singularities 2,3 . Quantum dynamics of such complex systems, except trivial cases, would be out-of-equilibrium when they are quenched in any finite time. However, such rich dynamical properties could lead to novel computational resources 4-7 provided that we obtain sufficient degree of control over their dynamics.
Adiabatic quantum computation (AQC) has been developed as a particular paradigm that utilize the continuous-time dynamics of driven many-body quantum systems for solving optimization tasks 4,5 . In this model, the solution of a hard combinatorial optimization problem is encoded in the ground state of an interacting many-body system which can be prepared adiabatically from an initially trivial ground state, provided that time evolution is much longer that the inverse of minimum gap square 5 . One of the major challenges to AQC, that has been largely ignored in the quantum computing literature, is that for many realistic problems the analog quantum annealer will inevitably contain a significant amount of quenched disorder smearing the corresponding quantum phase transitions for pure systems. Thus, the required time-scale for satisfying adiabatic limit could grow as a stretch exponential due to Griffiths singularity 3 , even in the absence of any first order phase transitions. The Griffiths effects have pronounced consequences for finitedimensional quantum systems, much stronger than in the classical counterparts. In fact, near-term quantum processors are best examples of low-dimensional quantum systems due to the inherent locality of physical interactions and geometrical constraints on the degree of connectivity 8 . After the embedding of a computational problem into quantum annealers, or their digital simulations 9 , they will inevitably react to quantum fluctuations inhomogeneously at the physical level. Consequently, near-term quantum processors will typically experience locally inhomogeneous and smeared first and second order phase transitions, even if we drive them with an external field which is homogeneous in space. In particular, annealing schedules exhibit multiple vanishing gaps between ground state and first excited state, see Fig. 1(a), leading to exponentially long annealing time-scales. In practice, we always have a finite annealing time-scale that would inherently violate the adiabaticity condition, even for finite-size systems, leading to emergence of domain walls or topological defects that emerge at a relatively wide effective quantum critical region. This is in sharp contrast to a single, well-defined quantum critical point for pure system, where their density of defects can be estimated via Kibble-Zurek mechanism (KZM) in the thermodynamics limit [10][11][12][13][14] . As of today, there is no known way to guarantee the quality of solutions, given finite space-time physical resources, and there is no constructive or algorithmic way to improve performance for such analog quantum information processors within a given accuracy. These issues have lead us to the following fundamental questions: Is it possible to engineer quantum phase transition in disordered systems by inhomogeneous control fields to enforce spatiallyinduced gaps between low energy sector and higher energy states (see Fig. 1(b))?
Here, we present a general approach for controlling quantum critical dynamics. We introduce different classes of spatial and/or temporal inhomogeneous protocols to drive strongly disordered quantum spin chains across a quantum phase transition and minimize the residual energy of the final state. This is achieved by creating governing Hamiltonian with multiple critical fronts that can synchronize the local phase transitions in space and time. In each local region, the number of spins that simultaneously experience the critical dynamics is controlled by the length scale and shape of the inhomogeneity in which the magnetic field is modulated. Causality is introduced as the main control strategy to spatially coordinate symmetry breaking events among neighboring regions by finding the appropriate degree of inhomogeneities and the speed of critical fronts to reduce the number of topological defects. We explore the conditions for an optimal suppression of domain walls and show that we can beat the standard homogeneous KZM prediction for the density of the topological defects for strongly disordered transverse Ising problem in 1D. Moreover, we show that these phenomena can similarly be observed for systems with k-local physical interactions. We demonstrate that inhomogeneous driving can be exponentially faster for such systems than conventional (homogenous) schemes such as adiabatic quantum annealing. Furthermore, we show that the universality of quantum critical phenomena holds for inhomogeneous quantum critical dynamics even in the presence of strong disorder.
The outline of this paper is as follows: In section I, we review causal origins of topological defects in the context of KZM for pure and disorder systems. In section II, we describe two general classes of inhomogeneous quantum annealing (IQA), type I and type II, for reconstructing phase transitions and present numerical results for strong-disorder 1D transverse Ising model. In this context, we show that AQC can be understood as a trivial form of either type I or type II IQA. In section III, we first provide a phenomenological theory of Exact numerical simulations of instantaneous eigenenergies of a random quantum Ising model as a function of evolution time for two distinct algorithms: (a) Homogeneous or standard adiabatic quantum evolution, (b) Inhomogeneous non-adiabatic quantum annealing. In (a) the instantaneous ground state energy is represented by red lines and it is set at 0. Yellow marks the instantaneous excited energy states that can be reached if evolution time-scale becomes comparable to the inverse of instantaneous gap. In contrast blue lines shows unaccessible excited eigenenergies due to vanishing instantaneous Hamiltonian state-to-state transition moments (see section III). In panel (b) it can be observed that the complexity of problem can be substantially changed via inducing certain casually-induced local gaps separating low-energy sector (red lines) from instantaneously accessible excited states (yellow lines). the emerging local gaps and its connection to threshold velocities for critical fronts. We then derive an expression for distribution of local gaps as a function of inhomogeneity slope with a logarithmic correction on the system size. We demonstrate universality of critical fronts shapes via strong-disorder renormalization group (SDRG) techniques. We also discuss how the shape of inhomogeneity is related to its penetration depth into disordered phase. A generalization of our work for k-local Hamiltonian system is presented in section IV. A detailed treatment of our work as a generalization of KZM and discussions on lower-and upper-bounds for the shape of critical fronts is provided in a separate manuscript 15 . The generalization to spin-glass systems will be presented in another subsequent work 16,17 .

I. CAUSAL ORIGIN OF TOPOLOGICAL DEFECTS
We start by reviewing the Kibble-Zurek mechanism (KZM) for pure systems (in absence of any disorders) which has been developed as the phenomenological theory to describe the breakdown of adiabaticity in critical systems [10][11][12][13][14] . The theory provides a rough estimate for the density of topological defects that arise when a quantum or classical many-body system is driven through a continuous critical point at a finite rate. The key observation is that in the vicinity of a quantum critical point a system at the thermodynamical limit effectively stops following the adiabatic evolution for any finite quench rate -no matter how slow it is driven. This results in emergence of universal KZ length scale which depends on the quench rate and manifests itself, among others, in the density of topological excitations.
The time dependent evolution of system can be expressed by a Hamiltonian as: is a control parameter with value g c at the critical point, and H p is the Hamiltonian of interest or the "problem Hamiltonian". Near a critical point the characteristic energy scale of the system behaves as ∼ 1/|ε| zν at the thermodynamical limit. The system experiences a divergence of the equilibrium relaxation time, τ = τ 0 /|ε| zν , as well as a divergence of the equilibrium correlation length, ξ = ξ 0 /|ε| ν , where ε = (g c −g)/g c is the dimensionless distance to the critical point. The ν and zν are the critical exponents that characterize the universality class of the phase transition. The derivation below assumes that the exponents are well defined, i.e. they do no dependent on ε and describe pure power-law dependence, and that there are no other relevant long-distance scales in the problem.
The speed of information, or the speed of second sound, is on the order of the ratio of critical length-scale to the critical time-scale v s ∼ ξ/τ = (ξ 0 /τ 0 )|ε| ν(z−1) . (1) A causal separation near a critical point for any pair of spins could emerge if their relative distance is much larger than length scale that the information can propagate with the corresponding second sound velocity, v s , for given finite quenching time interval. Consequently, choices of broken symmetry for spins belong to two different casual zones are not necessarily related. This is the origin of topological defects formation. The Lieb-Robinson bound 18 , which characterizes the maximum speed of information in quantum many-body systems with local interactions, provides an upper-bound for v s . We note that v s can achieve its Lieb-Robinson upperbound when z = 1, such as the prototypical 1D transverse Ising model. Within the vicinity of g c the quenched external field can be linearized in the form g(t) = g c (1 − t/τ Q ), such that ε(t) = t/τ Q , where τ Q is the quench rate and the critical point is crossed at t = 0. The parameter regime close to the critical point in which the system is not able to adiabatically adjust to the slowly changing external field, and effectively, to zero-th order approximation not responding, is called frozen or impulse regime. The freezing occurs at a particular time scalet in which the relax-ation time τ (t) becomes approximately equal to quench rate ε/ε. Thus, by setting τ (t) = |ε(t)/ε(t)| we arrive at This equation gives the KZ time-scale relevant to describe the universal behavior of the system slowly quenched though the critical point. The corresponding length-scale is a power-law of the quench rate as well This length scale can be used to estimate the size of the domains in the broken symmetry phase. Consequently, the density of defects is expected to vanish as d ∼ ξ(t) −D , where D is the dimensionality of system and we assume that the defects are sufficiently robust and do not relax quickly during the subsequent evolution. This is the key predication of KZM. For example, in the well-studied case of 1D Ising model in absence of any disorder we have ν = z = 1. The KZM prediction for the density of excitations reads d ∼ ξ(t) −1 ∼ τ −1/2 Q in that case [19][20][21] , which can indeed be verified analytically 20 . The above argument was later developed into full dynamical scaling hypothesis, which allows to obtain similar power-laws for other observables of interest [22][23][24][25][26] .
Understanding casual effects in disordered systems near a critical point and any attempt for estimation of density of defects requires careful analysis and not much is known outside of specific cases. Experimentally, quenches from the superfluid to the Bose glass were reported 27 , with the resulting residual energies vanishing very slowly with the increasing quench rate. Full theoretical understanding is still missing in this case. Theoretical investigations are mostly limited to the class of systems with the critical point in the universality class of so-called infinite-disorder fixed point. Here, we are interested in systems belonging to this class. We first consider the prototypical example of a random transverse Ising Hamiltonian for a chain of N spins,Ĥ = − N n=1 g(n)σ x n − N −1 n=1 J n,n+1 σ z n σ z n+1 , with quenched (fixed) disorder in the nearest-neighbors couplings J n,n+1 . In this article we assume that they are drawn from the flat distribution over interval [−1, 1]. The unit of time is set by = 1. Using strong-disorder renormalisation group (SDRG) techniques, the equilibrium properties of this model were first evaluated by Daniel Fisher 28 . For a homogenous or uniform transverse field in the model, the distribution of disorders induces a critical point that can be evaluated by relation g c = exp(log(|J n,n+1 |)). For uniform distribution of J n,n+1 ∈ [−1, 1] this yields a critical value of g c = e −1 0.367879. It should be pointed out that the critical point for similar systems in two dimensions 29 and in presence of dissipation 30 are also known to belong to this universality class. We use numerical SDRG to demonstrate universality of our non-equilibrium protocols in the section III. We also generalize our results to Ising model with certain k-local interactions in the section IV.
The presence of disorder, changes the universality class of the critical point of the Ising model from ν = z = 1 to ν = 2 and z → ∞, and thus quantitatively and qualitatively modifies the dependence of correlation length and density of defects on the quench time-scale. Most importantly, using SDRG techniques, it was evaluated that as the system approaches the critical point the gap of random Ising model scales as [ε] |ε| 1/|ε| [28], and consequently the critical exponent zν = 1/|ε| + O(1) diverges as ε → 0. For that reason the KZM derivation described earlier has to be modified to take this into account 31,32 . The characteristic time-scalet follows from the condition |ε(t)/ε(t)| = 1/(τ Q |ε(t)|) ≈ κ|ε(t)| 1/|ε(t)| , where κ is a constant factor on the order of one. The above relation can be solved in the limit of infinitely long annealing time, ln(τ Q ) 1, yielding 31 The density of defects is then suppressed logarithmically with quench time d ∼ 1/ ln 2 τ Q , which is quadratically faster than simulated annealing, where defects scale as d ∼ 1/ ln τ Q 33,34 . The existence of these logarithmic scaling laws implies that one has to run exponentially long annealing times to reduce the residual energy of the final state. However, as we will show in the next section one can recover a polynomial scaling by driving the system with a spatially inhomogeneous transverse field.

II. CAUSAL CONTROL OF TOPOLOGICAL DEFECTS WITH MULTIPLE CRITICAL FRONTS
From carefully studying defect formation under homogeneous drive fields, one can see how a new way of suppressing or controlling topological defects can emerge by being aware of causal separation of subsystems due to the extremely small values of velocity for information propagation near a critical point according to Eq. (1). In other words, one can try causal synchronization of the local phase transitions by inhomogeneous driving fields, as far as the critical front do not move faster than a threshold velocity corresponding to the speed of information, see [35][36][37][38] for a quantum case and 39-44 for classical counterpart. Note that this is fundamentally different than the standard annealing paradigm which is guided by the inverse of a global gap of a quantum Hamiltonian system which provides an upper-bound for relaxation time scales according to the adiabatic theorem. In other words, adiabaticity provides a sufficient condition for annealing time and it is not necessary to get low-energy states or even the ground state of disordered Ising systems.
Here we provide a phenomenological description of causally-induced non-equilibrium quantum phase transi-tions. Specifically, we develop an algorithmic quantum annealing approach to create a causal sequence of locally gapped Hamiltonians. We note that for strongly disordered systems in low dimensions there is a quantum Griffiths region that is spread in the disordered and ordered phases, i.e. on both sides of a critical point 28 . Within the Griffiths region the system undergoes effective local phase transitions that are space-time separated in nature even if the control fields are homogeneous. The key observation is that one can create situations in which the choices of symmetry-breaking events in a local neighborhood that have already experienced phase transitions earlier could influence the symmetry breaking events elsewhere, provided that the control fields have certain inhomogeneous spatiotemporal structures. These symmetry breaking events are perceived by the rest of the system, which is still in a disordered phase, as effective boundary conditions influencing their local fields.
In order to develop an algorithmic quantum annealer, here we construct a general class of inhomogeneous quantum annealing schedules. They are a function of a fixed total quench time or annealing time T ∼ τ Q , proportional to the annealing rate of the homogeneous quench τ Q introduced in the previous Section. The performance of the algorithms are evaluated by computing the precision Q of approximating the ground state. Here we mostly focus on the random instances of strongly-disorder spin chains, nevertheless our construction is general and can be applied to higher-dimensional systems 17 . There are two main reasons for such a choice. First, for 1D case we can simulate their dynamics exactly using a mapping to free-fermionic system, as e.g. in Ref. 31, 32, and 37. Also the critical behavior of such systems when driven via homogeneous external fields have been studied extensively, thus the new non-equilibrium physics of such systems when driven inhomogeneously can be better benchmarked and appreciated. The overall Hamiltonian for a system of N spins under a inhomogeneous driver field can be written as: and the quality of an output state is characterized by a normalized residual energy as: is the quantum state of the system at the final annealing time. |ψ gs is the ground state of the classical time-independent Hamiltonian, or the problem Hamilto- We note that for pure systems, where J nm = J, the normalized residual energy can be related to Kibble-Zurek correlation lengthξ by Q ∼ |J|ξ −D whereξ −D is the density of topological defects and D is the dimension of system. Here we assume that the inhomogeneous drive field is a transverse field that can be locally modulated for ev-ery individual spin. This Hamiltonian can be realized with the near-term quantum annealing technologies currently being developed at the D-Wave Quantum Computing Systems and Google Quantum AI Lab.
As we describe in the next section, for any given instance of disorders {J nm } as we drive the system toward the quantum critical point, the system responds to quantum fluctuations within M distinct "clusters", which are related to the emergence of rare local regions within the Griffiths phase. As we will show, the number and locations of clusters can be estimated via a simple preprocessing step that grows linearly with the size of the chain for 1D system. The generalization to higher dimensional system is presented in Ref. 17 .
In each cluster we drive the many-body system by a transverse Ising Hamiltonian with some local structure. Thus, we drive these M clusters simultaneously into some space-time separated inhomogeneous transitions, where each h l is the time-independent global magnetic field which has a spatial structure and each λ l (t) is spatially uniform but it can generates nonlinear dependence acterize various spatiotemporal dependencies of traveling quantum critical fronts, where ||n − n k || denotes a distance measure of node n from some center node n k per cluster where we trigger the quantum fluctuations. The center of these spatiotemporal inhomogeneities can be shifted linearly in time by v k (n, t)t with a spatiotemporal dependence for each k cluster. However, for simplicity, for rest of this work we concentrate on a constant critical front motion for each cluster; that is v k (n, t) = v k .
In the following section, we define the shape and two different kind of velocities for critical fronts and provide two simple examples of type I and II annealing, namely periodic inhomogeneous annealing, and mutliple-criticalfronts inhomogeneity.

A. Shape and velocity of critical fronts
The inhomogeneity slope and its horizontal and vertical velocities of inhomogeneity can be characterized by a set of hyper-parameters α, v h , v v corresponding, respectively, to local slope of the instantaneous field in space and it's spatial (horizontal) and temporal (vertical) velocities, that are defined by derivatives of g(n, t) and n(g f ix , t) as: v v (n, t) = −∂g(n, t)/∂t, α(n, t) = ∂g(n, t)/∂n. Thus, we can derive closed form expressions over these hyperparameters for different annealing schedules. To appreciate the generality of the shape of g(n, t), we consider two concrete and qualitatively distinct classes for inhomogeneous quantum annealing with respect to possible temporal and/or spatial inhomogeneities.

B. Type I inhomogeneous quantum annealing: space and time separated inhomogeneity
In this class, we consider a general form of independent or separated space and time quantum fluctuations to drive the annealing dynamics An example of this class will be a periodic spatial inhomogeneity (standing wave) combined with spatially-independent time-local inhomogeneity as: where each term in the spatial contribution in the second term corresponds to an estimated cluster size. We provide a simple illustration of these periodic spatial inhomogeneities in the next section. We note that KZM -in the context of pure systems -was also extended to quenches that are homogeneous in space, but nonlinear (inhomogeneous) in time 25,45 . Such inhomogeneity adjusts the quench rate to the distance from the critical.
Consequently it allows to reduce the number of generated defects.
C. Type II inhomogeneous quantum annealing: Spatiotemporal inhomogeneities In this class, we build a sufficiently general example by creating a mutliple-critical-fronts annealing schedule in M clusters where critical fronts in each cluster are moving with the speed v k (n) and each are governed by a separate activation function tanh where g c is the critical value of transverse field. We note that there is no particular significance for our choice of activation function here. As an important special class of the above driver field, we linearize the activation function in each cluster near quantum critical point, that is , then for each cluster we get: In the first example of this type, we consider an inhomogeneity with constant v k for each cluster of the In the next example of this type, we consider linear approximation of activation function near critical point

D. Standard AQC: Absence of any inhomogeneity
We note that the standard or homogeneous quantum annealing schedule which has been extensively studied and numerically benchmarked for almost two decades can be considered as the extreme limit of either type-I or type-II of IQA. In the former case we have one spatially uniform transverse field, g 0 (n) = const., and linear velocity λ 0 (t) = (1 − t/τ Q ), where τ Q is the overall annealing time-scale. Thus we have the familiar form of homogeneous transverse field, which is linear in time, as: g(n, t) = g 0 (1 − t/τ Q ). In order to see AQC as a extreme limit of type-II IQA, we must note that the homogeneous transfer field can be considered as a single critical front with a trivial flat shape with infinite velocity; that is θ → 0 and v(n) → ∞, while θv(n), which is basically the vertical velocity, will be equivalent to the inverse of annealing time and thus will be finite. The hyperparameters for homogenous annealing respectively become v v (n, t) = 1/τ Q , v h k (n, t) = ∞, and α(n, t) = 0.

E. Exponential suppression of defects
To illustrate the power of multi-front critical control, we numerically investigate two concrete forms of type I and II inhomogeneous annealing as described above and compare their performance against standard QA. All the simulations in this section are done using the Jordan-Wigner transformation that maps the Hamiltonian in Eq. (5) onto the system of free fermions where it can be solved numerically in a standard way. For details of these techniques, we refer the readers to the Appendix B of Ref. 37. For our numerical analysis here the cluster formation that we invoke is simple and has linear scaling with the system size for Ising chains. In the next section we use SDRG to examine construction and scaling of causal gaps.
Examples of the type I and II annealing schedules for one random instance of Ising chain are given in Fig. 2, in which different snap-shots of time-dependent transverse fields are plotted along the chain. In Fig. 2(a), we illustrate the trivial/homogeneous schedule. In Fig. 2(b) we explore the effects of periodic critical fronts by constructing the schedule where T is the total evolution time. Finally, in Fig. 2(c) we illustrate an example of multi-critical-point strategy.
In the latter, in order to decide the position of the cluster, as a pre-processing, we consider growing the clusters one spin at a time from one end of the chain. We keep the position of the weakest link in the cluster. After adding one spin, we check if the new J n is the new weakest one. We find the largest cluster satisfying min cluster k (J n ) · κ > v k , where v k is the velocity of the front for this cluster, and κ is a parameter fixing the exact value of the threshold. For a total fixed annealing time, all velocity v k are cluster dependent; allowing optimization of the computational parameters over the available time. For a cluster of size L k sites the (vertical) velocity is given by v v k = (|g f − g i | + αL k /2)/T for a given fixed total annealing time T . Each cluster is driven separately, and the inhomogeneous front is brushing from the middle of each cluster to both ends simultaneously. For each cluster spanning spins n + 1 . . . n + L the magnetic field is constructed as: The probability distribution of topological defects for these protocols are presented in Fig. 3. The defect density is evaluated for a strongly disordered instance of a spin chain consist of up to about 1000 qubits with J i,i+1 sampled randomly from [−1, 1]. Here, the inhomogeneous annealing can be regarded as a many-body quantum control strategy which can significantly reduce the number of topological defects by synchronizing the symmetry breaking events and can brush the reminder of defects into the weakest J ij where they act as defect sinks. In other words, not only do these non-adiabatic paths suppress the emergence of domain walls, but also minimize the energy cost per defect by several order of magnitudes. The scaling of annealing or quenching time as a function of inverse residual energy is plotted in Fig. 4 for 1000 random instances of 1D Ising chains ranging from 256 spins to 1024 spins. It can be observed that the annealing time is improved exponentially over standard AQC scheme for typical, 50% quantile, as well as harder instances, 99% quantile (each corresponds to the residual energy in which x% instances have smaller values). In the homogeneous protocol we have Q ∼ (log T ) −γ . We ob-tain γ ≈ 3.8±0.4. It can be compared with γ = 3.4 which has been obtained by from Caneva at. al. 32 from smaller values of quench times T and for slightly different protocol, where the magnetic field was also disordered. What is important here is that γ is larger than 2, i.e. the value of exponent governing the scaling of defect density. This results from defects being more likely to appear on the links with smallest |J n,n+1 |. For uniform distribution in [−1, 1], γ = 4 would correspond to defects appearing only at the weakest links.
On the other hand, in type II protocol we observe Q ∼ T −1.03 (fit for N = 1024, T > 100 and 50% quantile). Multi-front protocol have been constructed by fixed α = 1/8, where optimality of such choice is shown in panel (c) for T = 1000. Similar plots for other timescales (not shown) suggest that in this system this value is optimal independently if T 1. Panel (d) shows the mean number of clusters in type II annealing for different time-scales. For large times it should scale as ∼ N/T 0.5 which results from the preprocessing procedure. This follows from the expected size of clusters which can be solved in given T , given that J n,n+1 is drawn from uniform distribution in [−1, 1].

III. SCALING RELATIONS FOR CAUSAL GAPS
Here, we introduce the notation of causally gapped Hamiltonians that are created by time-dependent multiple critical fronts introduced in the previous section. In particular, we drive a scaling relation for the distributions of minimum casual gaps as the function of system size and inhomogeneity slope. The core idea is that there is an effective threshold velocity v k t that determines the suppression of topological defects formation within each cluster in disordered systems. If the velocity of the front in each cluster is much larger than this threshold speed, v k v k t for all k, then the effect of the spatial variations of the control field become irrelevant in the sense that we will recover the standard critical dynamics created by homogeneous driving, which can be understood by the standard KZM. However, when we drive each critical front such that v k v k t , the length scale and shape of the critical front becomes highly relevant allowing to suppress creation of the topological defects in each cluster. The shape of the critical front determines the number of spins that simultaneously experience criticality creating an effective (finite-size) energy gap. Otherwise, the homogeneous system would be gapless at the critical point in the thermodynamical limit.

A. Casual gap for pure systems
In the absence of disorder and with sufficiently smooth critical front, that is α 1, one can invoke a variant of KZM to estimate when the inhomogeneity of the driving front is relevant. This question can be regarded from two perspectives 35,36 . Firstly, starting from the limit of homogeneous driving, we note that the relevant speed of information at the critical front can be expressed asv ∼ξ/t. Here,ξ andt are the effective length scale and time scale that the system experience according to KZM given by the relations (2) and (3), respectively. This yieldsv = τ ν(1−z)/(1+zν) Q . Next we consider the control parameter ε(n, t) to be position dependent. We linearize the relation at the critical front for fixed position n f ixed as ε(n, t) = α(n f ixed − vt) = −αvt + const. (15) This gives the local annealing rate τ Q (n) = 1/(αv). Causality implies that if v v, then the choice of symmetry breaking which happened earlier along the chain cannot influence what is happening later and we recover the independent defect formation assumed in the standard KZM. The self-consistency condition allows to express such threshold velocity as a function of our main control parameter α. This is obtained by inserting the above annealing rate into the expression forv, which leads tov Alternatively, we could look at the instantaneous Hamiltonian resulting from inhomogeneous front in Eq. (15). We focus here on the instantaneous ground state of such system, which is interpolating -in space -between order and disorder phases. They are spatially separated by an effective critical regime which size, called a penetration depth, can be estimated aŝ It follows from a variant of KZM argument, so called KZM in space [46][47][48][49] , where one asks about characteristic distance from n f ixed up to which the system is able to locally adjust to ε(n) changing in space as if it were locally homogeneous with local correlation length determined by ε(n). Apart from the Ising model [46][47][48][49] , the interplay between inhomogeneous external field and criticality was studied in the context of spin−1 Bose-Einstein 50 , the XY model 51 , and the XXZ model 52 .
The finite size of the effective critical region in Eq. (17) allows to estimate the gap of the instantaneous Hamiltonian as∆ By combining those characteristic scales we obtain a threshold velocity v t which is again given by Eq. (16). The meaning of this threshold velocity is however different here. Namely, we can expect that if the velocity v in Eq. (15) is v v t , then the system would be able to follow its instantaneous ground state.
It should be noted that whenever z = 1 the threshold velocity becomes constant and we get a sharp suppression The critical fronts are moving with vertical velocity below the threshold and as a result, following the quench, we obtain a system without defects. In (b) we consider the setup with two clusters and two fronts merging in the middle of the chain. Again, the fronts velocity is below the threshold and there are no defects inside clusters. However, as the clusters were independent, defect can be created as they merge, reflecting the possibility that two independent clusters might break the symmetry differently. The probability of having such defect is equal 1/2 as can indeed be seen from the numerics. In (c) we consider the system with periodic boundary conditions. Effectively, in this setup we have one cluster which is selfmerging at the end of the quench. Importantly, such process does not lead to the creation of a defect -in contrast to the situation in panel (b).
of topological defects whenever we drive the system with critical front that is slightly below v t . The value of v t becomes equal to 2 for a transverse Ising chain when all couplings are equal to one 35 . This suppression of defects and causal synchronization, however, could be affected when we drive the system with a multi-front strategy. For a simple pure system driven by two critical fronts, both moving with velocity smaller than v t , we illustrate how the defects can be created when two clusters are merging together in Fig. 5, which highlights interplay of the causal effects and different boundary conditions on the defect suppression As we have discussed in section I, the disordered systems have a completely different critical phenomena as zν → ∞. In this case one has to modify Eq. (18) accordingly 37 by taking into account that at the critical point the gap is expected to vanish as a stretched exponential with the system size -effectively given in our case by Eq. (17). We elaborate on this in the next section in the context of multiple critical fronts driving strategy.

B. Universality of causal gaps via SDRG
In this section we use a combination of analytical and numerical SDRG techniques to show that the distribution of causal gaps are universal irrespective of the shape of inhomogeneities for strongly-disorder Ising systems. We also drive a scaling relation for the dependence of the minimal causally-induced gap on the actual system size and slope of inhomogeneity. We first present implementations of SDRG for disordered spin chains under various schedules.
The core concept of all RG techniques is to re-express the parameters which define a problem by coarse-graining some microscopic degrees of freedoms. In each step of the RG flow we arrive at effective Hamiltonian terms that have fewer and much simpler parameters acting on a lower energy and larger (macroscopic) length scale, such that certain physical or computational aspect of interest in the original problem remain unchanged. The procedure can be recursively repeated until the Hamiltonian is no longer changing which indicates that we have arrived at the fixed point of the RG flow. In SDRG, that is specifically developed for disordered systems and has been generalized to higer-dimensional systems 53 , the largest energy scales is systematically removed via two different operations: site decimation and bond decimation 28,54 , see 55 for a review.
A site decimation occurs whenever we have a sitedependent transverse field which is the largest energy scale within a local neighborhood of our system; that is g i > J ij ∀j. Site decimation means that we basically lock the spin i to the direction of its local transverse field. Such spin will be effectively decouples from the rest of the system. Emerging new couplings are generated between all neighbors of the decimated spin i. These effective couplings can be computed within second order perturbation theory asJ jk = Jij J ik gi , if J jk = 0 or otherwisẽ J jk = max(J jk , Jij J ik gi ). A bond decimation is performed in similar fashion. A bond decimation occurs whenever we have two sites i and j interacting via J ij that is the largest energy scale within a local neighborhood of our system; i.e. J ij ≥ g i and J ij ≥ g j ; and also J ij ≥ J ik ∀k  The RG flow starts with physical spins (as tree leaves represented by dots on a grey scale) and ends at the RG fixed point (roots) via other macrospins (color dots) in descending energy scales. The black physical spins are those that are already in the ordered or symmetry breaking phase while the grey means that they are still partially or fully in symmetric phase. The snapshot at g = 0.8 highlight how radically different these two strategies are: casually independent cluster formation in HQA and naturally growing of a central cluster nucleus in IQA. The other snapshots at g = 0.4, 0.37, 0.3 are chosen to be close to the critical point. The shaded blue area are casual zones that occurs for 0.2 SDRG energy cutoff for first two panels showing qualitatively and quantitively different cluster formations for IQA and HQA. In the last two panel the shaded grey area correspond to causal zones for the fixed point of SDRG. The key observation here is that HQA has several causally separated holes in the shaded area where certain physical spins cannot communicate their choice of symmetry breaking with other spins. These are the places that topological defects are highly likely. Nevertheless, within each causal zone forming at low-energy scale there is also some chance of defects for macrospins that are causally separated in higher energy scales (indicated by different colors within each causal zone). and J ij ≥ J lj ∀l. Bond decimation simply means that we lock the two sites i and j into a macrospin by projecting the combined pair into their local ground states. No additional bonds or coupling between any spins are generated in this case. All the spins that were previously coupled to at least one of the sites are now interacting to the combined cluster. For the spins that were coupled to both i and j we invoke a maximum selection rule. Effective transverse field at the emerging macrospins becomes g i = gigj Jij . We apply the above RG procedure to our time dependent Hamiltonian by considering each time as a snapshot for different instances of static spin chains, see Fig. 6. The SDRG simulation confirms our assumption of causally independent clusters in the homogeneous strategy, in the low-energy or long wavelength limit. In contrast, we observe inter-cluster casual dependence emerging in the multiple critical front strategy, see the SDRG visualizations of our protocols in Fig. 6.
Next, we discuss an upperbound for a global threshold velocity such that the multiple critical front strategy can lead to suppression of excitation between the low energy manifold and the rest of excited states. If we have full parallel annealing in all clusters simultaneously, we essentially interrupt the causality of symmetry breaking events between different clusters as illustrated in Fig 6. Consequently, we can consider each front independently. The transition matrix elements for each cluster k: where Ω k (n f ) defines local ground state to first excited state transition matrix elements for each cluster k and v k = dn k f dt denoting the cluster dependent vertical critical front velocity. For simplicity of analysis in following sections we choose a constant and uniform velocity for all the critical fronts in various clusters i.e. v k = v. Due to casual independence of clusters in multiple criticality, the adiabatic condition for low-energy states (approximate solutions) can be characterized by maximum over of all possible local transition matrix over its local gap, ∆ k (n f ); that is k (n f ) . Then the low-energy quasi- Universal scaling of the gap distributions with respect to the shape of inhomogeneity. In order to achieve smooth shapes of transverse fields at the border of various clusters, we consider a type-I inhomogeneity (standing critical wave) with 4 clusters in a system of N = 1024. By comparison to Fisher and Young SDRG simulations for universality of quantum Ising models 56 , we can conclude that our inhomogeneous protocol has changed the effective length scale of system from N to O(α −2/3 ). Alternatively, we can consider that multiple critical front strategy has changed the critical dynamical exponent of full system from z → ∞ to z → 0; that is the system always stay close to a critical point during almost the entire evolution, however it never experience a true quantum phase transition on the whole system (shortcutting adiabatic evolution).

adiabatic dynamics is expected for
In the following we are going to discuss the universality of the threshold velocity, reflected in its dependence on the shape of critical front characterized by the slope α. Firstly, it is worth drawing the connection between the condition in Eq. (20) with the threshold velocity which was derived in the previous section. There, it was calculated asv =ξ i /t i , whereξ i andt i were the characteristic length scale (the penetration depth) and time scales (given by the inverse of the gap) related with the slope of inhomogeneity α. Similarly, the adiabatic condition is sometimes formulated as, v ∆ · Γ, where ∆ is the energy gap, and Γ estimates the width of the region (in the driving parameter space) for which the gap is close to its minimal value, see e.g. 57 -in analogy to the avoided level crossing and the Landau-Zener problem. To resolve this seeming inconsistency (i.e. the gap appearing in Eq. (20) is squared comparing to other expression), we note that Ω k (n k ) and ∆ k (n k ) are not independent. We expect that Ω k (n k )/∆ k (n k ) ∼ ξ −1 i , i.e. it is directly proportional to the size of the effective critical region. We show this in Fig. 7 and discuss further below. Importantly, this allows us to focus on the scaling of the gap.
The effective size of the critical region (for a given front) is given by Eq. (17). In our case ν = 2, which gives ξ i ∼ α −2/3 . The fronts are independent if their respective distances, or the sizes of the clusters, are ξ i . The typical gap related with such critical front is then expected to scale as a stretched exponential in ξ i , ∆ k (n k ) ∼ e −const √ ξi ∼ e −const·α −1/3 37 .
To be more precise, we consider P (∆ n f , α) as a distribution of the instantaneous gaps ∆ n f , as each front is traveling within a cluster for a fixed α. We argue that we can observe the universality of local gap probability distributions by introducing a rescaled log-gap x = α 1/3 log ∆, with the expected P (∆, α)d∆ =P (x)dx, and universal distributionP (x). In contrast, the universal scaling of the gap in the finite size homogeneous system at criticality is known to be described by universal distribution P (∆, N )d∆ =P (x)dx with x = N −1/2 log ∆ 56,58 . Our ansatz is then a directly consequence of an assumption that the effective size of the critical region in our protocol is no longer given by N but instead is characterized by the penetration depth in Eq. (17).
We verify those scaling predictions in Fig. 7. We consider the system with two clusters, i.e. 4 critical fronts to highlights the independence of the front and that the scaling prediction naturally carry on to the case of multiple fronts. We calculate the minimal relevant gap (related with one of the fronts), which we distinguish by finding the minimal energy eigenstate with Ω n f above some threshold. We calculate the gap via numerical SDRG which ends at a system of few spins which is subsequently exactly diagonalized. While the SDRG procedure which we use is introducing some errors, we check that the results can be essentially reproduced by the numerically exact solutions based on free-fermionic picture for transverse-field Ising models. Using SDRG allows us to highlight that the fronts are independent as they become effectively decoupled during SDRG flow. We compute the distribution of ∆ n f by sampling, for each disorder instance, from different equidistant front positions as they are moving though the chain. We disregard the beginning and the end of the protocol when the system is far from criticality and focus on the relevant intermediate part where we have independent critical fronts. In Fig. 7 we collect the statistics using 5000 disorder instances. Indeed, we can see that the peaks of the rescaled distribution collapse validating the scaling anzats. The tail corresponding to small energies, where the distributions do not properly collapse, is non-universal and results from occasional very weak links (as J ∈ [−1, 1]). They give rise to gaps of similar order (occurring again when the magnetic field acting on the neighboring spins is again almost switched off), which are characteristic for strictly 1D system. We elaborate more on this later in this section. The presence of such gaps is especially pronounced for larger values of α when the typical gap ∼ e −const·α −1/3 , which can be attributed to many-body effects, is larger.
The analysis above regards local, instantaneous gap re- Fig. 7, especially in the limit of small α when gaps attributed to weak links are less relevant, indicates that the above holds also in our case 60 . Now, let the probability that ∆ is smaller than some ∆ q min be where x q = log(∆ q min )α 1/3 . In order to find the probability distribution for minimal gap we assume that we are sampling N times from P (∆, α) (or more precisely proportional to N times in the limit of large N ). That way, the probability where ∆ min is the minimum from the sample of N instantaneous gaps. This equation defines ∆ q min as a q-quantile for the global minimal gap. Now, we obtain ∆ q min from the above equation. This gives 1 N log q = log(1 − ε) − ∼ 1 2axq e −ax 2 q , which is obtained by expanding the error function for the Gaussian tail in Eq. (21), to the leading order in small x q . Solving this equation in the leading order we obtain This suggests that if we fix a quantile q, then in the asymptotic limit i.e. it is vanishing slower than any polynomial with increasing N . This has to be compared with homogeneous gap scaling as: ∆ q min ∼ exp(−const · √ N ), which vanishes as stretched exponential with N . We note here that similar analysis in case whenP (x) would have exponential tail for large negative x would give polynomial dependence for the minimal gap on N .
In order to fully analyse the tail of the gap distribution we also consider the situation that the minimum gap is enforced by very small local link rather than manybody gap of critical system of effective size ξ i . This issue, which is essentially an artifact of 1D systems, can be largely avoided by a multi-front strategy where such links are placed in-between clusters which are driven quasiadiabatically. To that end, let's assume that the links are drawn from uniform distribution J i ∈ [−1, 1]. Probability that a single link is weaker then some is equal , or P (|J i | > ) = 1 − . Let's consider that q is the probability that all the links are stronger than this is q = (1 − ) N . For small this yields − log q N .
The minimal gap related with such a weak link is of similar size, and such effects become relevant in 1D geometry. Finally, we should note that if we consider the distribution of the logarithm of the gap x, the uniform distribution and related weak links directly translate to the exponential tail of P (x) -mentioned in the previous paragraph -which is indeed still visible in Fig. 7 for larger values of α. Finally, in Fig. 8 we illustrate the relation between the transition matrix elements Ω n f and the gap ∆ n f . To that end, for simplicity we consider protocol with single front and, as in Fig. 7, sample the values of minimal relevant gap and corresponding Ω n f as the front is traveling though the chain. The results are collected as probability distribution, where the statistics is collected from 5000 instances. In order to illustrate the universal behavior we employ rescaled variables. For the gap x = α 1/3 log(∆ n f ) as above, and y = log(α −2/3 Ω n f /∆ n f ) to reflect the expected relation Ω n f /∆ n f ∼ ξ −1 i ∼ α 2/3 . We plot the obtained distributions of P (x, y) in Fig. 8 for several different values of α. We observe that they are roughly similar in agreement with our prediction. We should note that local maxima of Ω n f /∆ n f -i.e. where it is most relevant -coincide with the local minima of energy gap. Apart from those point Ω n f /∆ n f is quickly approaching zero, which is reflected by elongated shape of the distribution P (x, y) in the direction of small (irrelevant) y.

IV. RANDOM CLUSTER-ISING HAMILTONIAN
The numerical results in the previous sections were confined to 1D geometry where, in the final ground state, the neighboring spins are aligned according to corresponding J n,n+1 . In this section we show that the general approach discussed in this article does not hinge on the possibility to align nearest-neighbor interacting spins in strictly 1D Ising geometry. To that end we consider random cluster-Ising Hamiltonian as follows: where the first two terms contains the problem Hamiltonian and the last term is the (inhomogeneous) driving term given by external transverse field. The quench dynamics generated by this Hamiltonian can be simulated analogously as for the Ising chain from the previous section as -using Jordan-Wigner transformation -it can be mapped onto a free fermionic system. First, lets consider the pure and homogeneous model with J n = J > 0, K n = K > 0 and g n = g > 0. Depending on the relative strength of those terms the model has 3 distinct phases, see e.g. [61][62][63][64][65] . It is convenient to set K+J = 1. When the magnetic field is dominating, g > 1, the system is in paramagnetic phase analogously to the Ising model. We are going to initialize the evolution in this phase. In the opposite limit, g < 1, the system is in ferromagnetic phase for J + g > K and symmetryprotected topological order phase for K > J + g. For J = g = 0 it reduces to the cluster Hamiltonian 66 . The phases are separated by critical points with z = 1 or z = 2.
We consider random couplings J n and K n , which makes the final target state far from being trivial for a 256 spin system. We present the results of the quench in Fig. 9. We consider two different disorder distributions, where in Fig. 9(a) J n are dominating; that is J n ∈ [−0.75, 0.75] and K n ∈ [−0.25, 0.25], and in Fig. 9(b) all J n = 0 and K n ∈ [−1, 1]. For homogeneous driving the residual energy is vanishing logarithmically with the quench time in both cases, which is similar behavior as for the random Ising model. It can be observed that there is crossover of the performance for sufficiently long times and the homogeneous protocol  9. Scaling of the total annealing time as a function of inverse residual energy for homogeneous control (blue lines) versus inhomogeneous protocol: Here we use two fronts spreading from the center of the chain for N = 256. The Ising Hamiltonian in Eq. (26) has generally three different paramagnetic, ferromagnetic and topological insulator phases. In (a) the Hamiltonian contains 2-and 3local random terms versus in (b) only 3-local random terms are present. It can be observed that in both cases non-zero α significantly outperform the homogeneous or α = 0 schedule, with an exponential advantage of optimal inhomogeneous protocol α = 1/8 when we have an interplay of two-and threebody interactions.
is considerably outperformed by an inhomogeneous driving fields with the optimal slope. The advantage for the case of random J n could be exponential, see Fig. 9(a). In this case, the spatial inhomogeneity allows the system to reach the quality of solution (small residual energies) which are practically unattainable within the homogeneous approach. Here, we use a version with single cluster and two critical fronts spreading from the center of the chain. It should be noted that different optimal shape (α) in those two cases are different. This optimal value of α can be found numerically for given distribution of disorders as a simple preprocessing, or hyper-parameter characterization, similar to the spirit of finding the optimal annealing time or number of sweeps for simulated annealing or quantum Monte Carlo solvers.

V. CONCLUSION AND FUTURE WORKS
We have presented a model for engineering quantum phase transitions in disordered systems by manipulating information flow among clusters that are formed within a quantum critical region. We have shown that space-time inhomogeneities in the control fields could lead to reconstruction of casual zones (light cones), such that symmetry breaking events can be synchronized suppressing the density of topological defects and/or redistributing their spatial arrangements. We have used exact diagonalization techniques for 1D systems to show an exponential speedup of non-adiabatic inhomogeneous quantum annealing over standard adiabatic quantum computing, even in the presence of higher order interactions. By application of renormalization group techniques we have demonstrated that the effective causal gaps exhibit universality with respect to the shape of inhomogeneity. We have derived a scaling relation showing such effective gaps have sub-polynomial scaling with the system size, in contrast to stretch exponential for homogeneous control strategies. In a subsequent work 17 , we will provide a detailed theoretical discussion of our work as a generalization of KZM for disordered systems including various bounds for the shape of critical fronts and threshold velocities under different assumptions. We will also discuss how our approach can be applied to lowdimensional spin-glass problem Hamiltonians 16,17 . During the preparation of this manuscript an exponential speedup for inhomogeneous quantum annealing of p-spin model was reported showing ferromagnetic first-order phase transitions can be smeared with inhomogeneous control strategies 67 .