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

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–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–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–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–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.

Zoom In Zoom Out Reset image size

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.

We consider a chain of N spins described by the random Ising Hamiltonian

$$\begin{eqnarray}&&\hat{H}=-{ \mathcal J }\displaystyle \sum _{n=1}^{N-1}{J}_{n}{\sigma }_{n}^{x}{\sigma }_{n+1}^{x}-{ \mathcal J }\displaystyle \sum _{n=1}^{N}{g}_{n}{\sigma }_{n}^{z},\end{eqnarray} \tag{ 1 }$$

with quenched disorder (constant in time) in the nearest-neighbors couplings J_n. We set ${ \mathcal J }={\hslash }=1$ in subsequent calculations, that is equivalent to use ${\hslash }/{ \mathcal J }$ as a unit of time. Realizations of disorder are drawn from a uniform distribution

$$\begin{eqnarray}P({J}_{n})=\left\{\begin{array}{cc}1 & \mathrm{for}\,{J}_{n}\in (1/2,3/2),\\ 0 & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 2 }$$

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 $\mathrm{log}({g}_{c})=\overline{\mathrm{log}(| {J}_{n}| )}$ for uniform g_n = g_c ($\simeq \,0.9558\ldots$), and belongs to the universality class of the infinite-randomess 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}_{n}({t}_{i})={g}_{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}_{n}({t}_{f})={g}_{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}_{n}\geqslant 0.5$ 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–7]. The mechanism exploits the divergence of the equilibrium relaxation time $\tau [\varepsilon ]={\tau }_{0}/| \varepsilon {| }^{z\nu }$, known as critical slowing down, as a function of the dimensionless distance to the critical point $\varepsilon =({g}_{c}-g)/{g}_{c}$. In the proximity of g_c the modulation of the magnetic field can be linearized in the form $g(t)={g}_{c}(1-t/{\tau }_{Q})$, so that $\varepsilon =t/{\tau }_{Q}$, where ${\tau }_{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 $\hat{t}={({\tau }_{0}{\tau }_{Q}^{z\nu })}^{\tfrac{1}{1+z\nu }}$, and satisfies $\tau (\hat{t})=| \varepsilon /\dot{\varepsilon }|$. KZM estimates the size of the domains in the broken symmetry phase using the equilibrium value of the correlation length $\xi [\varepsilon ]={\xi }_{0}/| \varepsilon {| }^{\nu }$. At the freeze-out time, it scales as a power-law of the quench rate, e.g., $\hat{\xi }=\xi [\varepsilon (\hat{t})]={\xi }_{0}{({\tau }_{Q}/{\tau }_{0})}^{\tfrac{\nu }{1+z\nu }}$. In the Ising model without disorder, $\nu =z=1$, and the KZM prediction for the density of excitations reads $d\sim \xi {(\hat{t})}^{-1}\sim \tfrac{1}{\sqrt{{\tau }_{Q}}}$ [41–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., $Q=\langle \hat{H}{\rangle }_{\tau }-{E}_{\mathrm{gs}}$. For the Ising model, the value of Q after the quench scales in the same way as the density of excitations.

The presence of disorder, that drives the universality class of the critical point towards the infinite-randomnes fixed point, severely modifies the dependence of the relaxation time on the distance from the critical point. It is found that $\tau [\varepsilon ]\simeq {\tau }_{0}/| \varepsilon {| }^{1/| \varepsilon | }$, where the critical exponent $z=1/2| \varepsilon | +O(1)$ effectively diverges as the system approaches the critical point ($z\to \infty$ as $\varepsilon \to 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 $d\sim 1/{\mathrm{ln}}^{2}{\tau }_{Q}$ in a slow transition (for large ${\tau }_{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

$$\begin{eqnarray}{g}_{n}(t)=\left\{\begin{array}{ll}{g}_{i}, & n-{vt}\gt \tfrac{{g}_{i}-{g}_{f}}{2\alpha },\\ \tfrac{{g}_{i}+{g}_{f}}{2}+\alpha (n-{vt}), & | n-{vt}| \leqslant \tfrac{{g}_{i}-{g}_{f}}{2\alpha },\\ {g}_{f} & n-{vt}\lt \tfrac{{g}_{f}-{g}_{i}}{2\alpha },\end{array}\right.\end{eqnarray} \tag{ 3 }$$

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}_{f}={vt}$ 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 $\alpha \to 0$ and $v\to \infty$ while keeping $\alpha v={\tau }_{Q}^{-1}$ fixed, and (ii) driving of one spin at a time, when $\alpha \approx 1$ [47].

In absence of disorder, the KZM can be extended to account for an inhomogeneous scenario [27–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\gg {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\ll {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 $\alpha \ll 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–33]

$$\begin{eqnarray}&&{\hat{\xi }}_{i}\sim {\alpha }^{-\tfrac{\nu }{\nu +1}}.\end{eqnarray} \tag{ 4 }$$

When the gap at the critical point vanishes polynomially with the system size this leads to the opening of instantaneous gap that scales as ${\hat{{\rm{\Delta }}}}_{i}\sim {\hat{\xi }}_{i}^{-z}\sim {\alpha }^{\tfrac{z\nu }{\nu +1}}$. By combining the characteristic time and length scales in the problem one can then estimate the threshold velocity as

$$\begin{eqnarray}&&{v}_{t}\sim {\alpha }^{\tfrac{(z-1)\nu }{\nu +1}}.\end{eqnarray} \tag{ 5 }$$

In the Ising model without disorder z = 1 and v_t is a constant independent of α ($\alpha \ll 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\lt {v}_{t}$ and the effectively homogeneous regime for $v\gt {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\to \infty$ 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.

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,

$$\begin{eqnarray}&&\phantom{\rule{}{1.0ex}}| \langle 0,t| \displaystyle \frac{{\rm{d}}\hat{H}}{{\rm{d}}{t}}| 1,t\rangle \phantom{\rule{}{1.0ex}}| =\phantom{\rule{}{1.0ex}}| \langle 0,t| \displaystyle \frac{{\rm{d}}\hat{H}}{{{\rm{d}}{n}}_{f}}| 1,t\rangle \phantom{\rule{}{1.0ex}}| \displaystyle \frac{{{\rm{d}}{n}}_{f}}{{\rm{d}}{t}}={\rm{\Omega }}({n}_{f})v,\end{eqnarray} \tag{ 6 }$$

where we have defined ${\rm{\Omega }}({n}_{f})=| \langle 0,t| \tfrac{{\rm{d}}\hat{H}}{{{\rm{d}}{n}}_{f}}| 1,t\rangle |$. The instantaneous energy gap ${{\rm{\Delta }}}_{{n}_{f}}$ can be parameterized both as a function of the front position n_f  and the time of evolution t,

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{n}_{f}}={\rm{\Delta }}({n}_{f})={\rm{\Delta }}(t)={E}_{1}(t)-{E}_{0}(t),\end{eqnarray} \tag{ 7 }$$

where ${E}_{\mathrm{0,1}}(t)$ are the energies of the instantaneous ground state $| 0,t\rangle$ and the first excited state $| 1,t\rangle$ of Hamiltonian (1), respectively. Since the parity operator $\hat{P}={\prod }_{i=1}^{N}{\sigma }_{n}^{z}$ 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

$$\begin{eqnarray}&&\displaystyle \frac{| \langle {\rm{d}}\hat{H}/{\rm{d}}{t}\rangle {| }_{\mathrm{0,1}}}{{\rm{\Delta }}{(t)}^{2}}=\displaystyle \frac{v{\rm{\Omega }}({n}_{f})}{{\rm{\Delta }}{({n}_{f})}^{2}}\ll 1,\end{eqnarray} \tag{ 8 }$$

which allow us to estimate the value of threshold sweep velocity v_t below which the dynamics of the quench is effectively adiabatic as

$$\begin{eqnarray}&&v\ll {v}_{{n}_{f}}={v}_{t}({n}_{f})=\displaystyle \frac{{\rm{\Delta }}{({n}_{f})}^{2}}{4{\rm{\Omega }}({n}_{f})}.\end{eqnarray} \tag{ 9 }$$

The factor of 4 in the definition of v_t above is introduced so that ${v}_{t}({n}_{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

$$\begin{eqnarray}&&{\rm{\Omega }}({n}_{f})=| \langle 0,{n}_{f}| \displaystyle \sum _{n}{g}_{n}^{\prime }({n}_{f}){\sigma }_{n}^{z}| 1,{n}_{f}\rangle | ,\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}{g}_{n}^{\prime }({n}_{f})=\left\{\begin{array}{ll}0, & | n-{n}_{f}| \gt \tfrac{{g}_{i}-{g}_{f}}{2\alpha },\\ -\alpha , & | n-{n}_{f}| \leqslant \tfrac{{g}_{i}-{g}_{f}}{2\alpha }.\end{array}\right.\end{eqnarray} \tag{ 11 }$$

This analysis is restricted to quasi-adiabatic dynamics governed by adiabatic following of the ground-state and $| \,0\rangle \,\mbox{--}$ $| 1\rangle$ transitions. As a result, it is expected to fail when transitions to higher excited states are dominant, e.g., for large α ($\approx \,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 $\alpha =1/32$. While the gap fluctuates as the front travels through the chain, it remains finite. We define ${{\rm{\Delta }}}_{\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 $\alpha \gt 0$. In that case, the critical front can be thought of as probing different local realizations of disorder where the effective size of the system ${\hat{\xi }}_{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 ${{\rm{\Delta }}}_{\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 ($\alpha =0$) the typical energy gap at the critical point is expected to vanish as $\sim \exp (-C\sqrt{N})$ [38, 39, 48, 49], where C is a nonuniversal constant.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 2(b) we show the mixing term ${\rm{\Omega }}({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}_{t}({n}_{f})$ below which the evolution is expected to stay adiabatic. We define v_t as the minimum of ${v}_{t}({n}_{f})$ for a given realization; see figure 2(c), where its value is of order unity for $\alpha =1/32$. Since ${v}_{t}({n}_{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 ${{\rm{\Delta }}}_{\min }$, the largest value of mixing matrix denoted by ${{\rm{\Omega }}}_{\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 ${{\rm{\Omega }}}_{\max }\sim {\alpha }^{1.1}$ similar to the Ising case with no disorder where ${\rm{\Omega }}\simeq \alpha$. The gap ${{\rm{\Delta }}}_{\min }$, however, disappears faster than polynomially in that limit which results in a monotonic dependence of v_t on α (without disorder ${{\rm{\Delta }}}_{\min }\sim {\alpha }^{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 ${\rm{\Delta }}({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 $P({{\rm{\Delta }}}_{{n}_{f}},\alpha )$ the distribution of the instantaneous gap for different α, normalized according to ${\int }_{0}^{\infty }{\rm{d}}{{\rm{\Delta }}}_{{n}_{f}}P({{\rm{\Delta }}}_{{n}_{f}},\alpha )=1$. For homogenous driving ($\alpha =0$) of a finite system at the critical point the equivalent distribution over realizations of disorder is universal if one considers the scaling variable $\theta ={N}^{-1/2}\mathrm{ln}{\rm{\Delta }}$ [38, 39, 48, 49]. In the spirit of the inhomogeneous KZM [30], for an inhomogeneous system the characteristic length scale which governs the behavior of the system is expected to scale as ${\hat{\xi }}_{i}\sim {\alpha }^{-\tfrac{\nu }{1+\nu }}\sim {\alpha }^{-2/3}$, see equation (4). This suggests that the relevant scaling variable for small finite α is

$$\begin{eqnarray}&&{\theta }_{{\rm{\Delta }}}={\hat{\xi }}_{i}^{-1/2}\mathrm{ln}{{\rm{\Delta }}}_{{n}_{f}}={\alpha }^{1/3}\mathrm{ln}{{\rm{\Delta }}}_{{n}_{f}},\end{eqnarray} \tag{ 12 }$$

where ${\hat{\xi }}_{i}$ plays the role of an effective size of the critical system. For the random Ising model, the critical exponent $\nu =2$ [39] characterizes the equilibrium value of $\xi [\varepsilon ]$, that describes mean correlations dominated by rare pairs of strongly correlated spins and should be relevant for the low energy part of the spectrum.

In order to numerically verify equation (12), for each realization of disorder we sample values of ${{\rm{\Delta }}}_{{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 $P({\theta }_{{\rm{\Delta }}},\alpha )$ as a function of α. Depending on α the statistics is built from $\gt {10}^{6}$ points from $\gt 2000$ realizations of disorder. The distributions for different $\alpha \ll 1$ collapse onto each other corroborating our scaling prediction; as seen for $\alpha \leqslant 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 $\mathrm{log}({\rm{\Delta }})/\sqrt{0.46N}$ is used⁵ . The factor of $\simeq \,0.46$, which we found by collapsing the distributions in figure 4, can be interpreted as the prefactor in the scaling ${\hat{\xi }}_{i}\sim {\alpha }^{-2/3}$, that sets the effective size of the critical region.

Zoom In Zoom Out Reset image size

On the other hand we observe that for larger values of α, such as e.g. $\alpha ={2}^{-4}$ in 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 ${{\rm{\Delta }}}_{{n}_{f}}\sim \exp (-C{\alpha }^{-1/3})$ in the limit $\alpha \to 0$. Here, we are mostly interested in the minimal gap ${{\rm{\Delta }}}_{\min }$, which would be determined by the behavior of the tail of $P({\theta }_{{\rm{\Delta }}},\alpha )$ 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 $P({\theta }_{{\rm{\Delta }}},\alpha )$ in a continuous way; see figure 2(a). We can estimate the dependence of the minimal gap ${{\rm{\Delta }}}_{\min }$ on the system size N by assuming that, to leading order, N instances are drawn from the distribution $P({\theta }_{{\rm{\Delta }}},\alpha )$ 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.

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.

We note that the total time required for the inhomogeneous protocol to be completed reads

$$\begin{eqnarray}&&T=\displaystyle \frac{N}{v}+\displaystyle \frac{| {g}_{i}-{g}_{f}| }{\alpha v},\end{eqnarray} \tag{ 13 }$$

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 $| {g}_{i}-{g}_{f}| {\tau }_{Q}$ (in this limit $v\to \infty$ as $\alpha \to 0$). By contrast, in the strongly inhomogeneous limit (when $N\gg | {g}_{i}-{g}_{f}| /\alpha$) 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.

Zoom In Zoom Out Reset image size

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}$, $\alpha \approx 1/16$ is optimal for harder cases with $q=99 \%$ and $\alpha \approx 1/32$ is optimal for $q=50 \%$ ⁶ . 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 $\alpha \to 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\sim 1/\mathrm{log}{(T)}^{2}$ [37, 38], making it unpractical to reach the adiabatic limit. The inhomogeneous driving approaches this limit for non-zero α 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 $\alpha =1/32$ for a time-scale of $T={10}^{4}$, for different system sizes. As 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.

Zoom In Zoom Out Reset image size

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.

We acknowledge funding support from UMass Boston under project P20150000029279 (AdC) and Narodowe Centrum Nauki (NCN, National Science Center) under Project No. 2013/09/B/ST3/01603 (MMR), and Quantum Artificial Intelligence Laboratory at Google (MM).

In this section, we show the results of the analysis based on adiabatic theorem, see section 3, applied to the case without disorder, J_n = 1. In figures A1 (a)–(c) we show, respectively, the instantaneous gap, mixing term ${\rm{\Omega }}({n}_{f})$ and, the estimation on local value of threshold velocity (from the combination of the two using equation (9)), for the slope of the front $\alpha =1/32$ and N = 512. They can be directly compared with the disordered case presented in figure 2.

Zoom In Zoom Out Reset image size

We define the threshold velocity based on the value in the bulk (i.e. when the front is inside the chain)—we neglect here a small peak appearing in figure A1(c) when the front is entering the chain. The scaling of the gap in the bulk can be derived analytically, ${\rm{\Delta }}\simeq \sqrt{8\alpha }$ ($\alpha \ll 1$) [34]. For small α we fit ${\rm{\Omega }}\simeq \alpha$. Accordingly, as can be seen in figure A1(d), the vanishing of the gap with decreasing α is compensated by the vanishing of the mixing term that results in the saturation of v_t to a constant value for $\alpha \ll 1$. This can be contrasted with the disordered case—figure 3(c)—where the threshold velocity is monotonically decreasing with (small) α. Such a dependence is expected for models with finite critical exponent $z\gt 1$, see equation (5) [35].

The threshold velocity can be found analytically [34] as v_t = 2 and equals the largest velocity of quasiparticles at the critical point of the clean Ising model. We fix the factor of 4 in the denominator of equation (9) to match the correct value of the threshold velocity in the limit of $\alpha \ll 1$. In that limit v_t gives a sharp boundary between $v\gt {v}_{t}$ where the system gets excited when the front is traveling through the chain, and $v\lt {v}_{t}$, where there are no excitations appearing in the bulk.

Figure A1(d) also indicates that for $\alpha \approx 1$ the threshold velocity is becoming smaller as the scaling prediction, based on assumption that the front is smooth enough, $\alpha \ll 1$, does not longer hold. This is as well the limit where the first excited state is not well separated from the rest of the spectrum and the estimation of threshold velocity as presented in section 3 is expected to break down.

The Hamiltonian in equation (1) is solved in a standard way, exploiting its mapping onto a system of free fermions via the Jordan–Wigner transformation ${\sigma }_{n}^{z}=1-2{c}_{n}^{\dagger }{c}_{n};$ ${\sigma }_{n}^{x}+{\rm{i}}{\sigma }_{n}^{y}=2{c}_{n}{\prod }_{m\lt n}(1-2{c}_{m}^{\dagger }{c}_{m})$, where c_n are fermionic annihilation operators. It is convenient to introduce Majorana fermions ${a}_{2n-1}={c}_{n}+{c}_{n}^{\dagger }$, ${a}_{2n}={\rm{i}}({c}_{n}-{c}_{n}^{\dagger })$, which are Hermitian and satisfy anticommutation relations $\{{a}_{m},{a}_{n}\}=2{\delta }_{m,n}$. In this base the Hamiltonian reads:

$$\begin{eqnarray}\hat{H} & = & -\displaystyle \sum _{n=1}^{N-1}\left(\displaystyle \frac{{{\rm{i}}{J}}_{n}}{2}{a}_{2n}{a}_{2n+1}+{\rm{h.c.}}\right)\\ & & -\displaystyle \sum _{n=1}^{N}\left(\displaystyle \frac{{{\rm{i}}{g}}_{n}}{2}{a}_{2n-1}{a}_{2n}+{\rm{h.c.}}\right)\\ & = & {\hat{H}}_{J}+{\hat{H}}_{g}.\end{eqnarray} \tag{ B1 }$$

For convenience we introduce $\hat{{H}}={\vec{a}}^{\dagger }{H}\vec{a}$, where $\vec{a}$ is a column vector consisting of operators a_n and ${H}$ is a $2N\times 2N$ matrix. We separate the matrix describing full Hamiltonian ${H}={{H}}_{J}+{{H}}_{g}$ into parts corresponding respectively to the coupling and magnetic field, which are both block diagonal,

$$\begin{eqnarray*}{{H}}_{g} & = & \underset{n=1}{\overset{N}{\displaystyle \bigoplus }}\left(\begin{array}{cc}0 & -{{\rm{i}}{g}}_{n}/2\\ {{\rm{i}}{g}}_{n}/2 & 0\end{array}\right),\\ {{H}}_{J} & = & 0\oplus \underset{n=1}{\overset{N-1}{\displaystyle \bigoplus }}\left(\begin{array}{cc}0 & -{{\rm{i}}{J}}_{n}/2\\ {{\rm{i}}{J}}_{n}/2 & 0\end{array}\right)\oplus 0.\end{eqnarray*}$$

The static properties of (instantaneous) system are found by numerical diagonalisation of matrix ${H}$. This amount to employing Bogoliubov transformation to a new base of Majorana fermions $\vec{a}={{O}}_{0}\vec{b}$, where orthogonal matrix ${{O}}_{0}$ brings the Hamiltonian into canonical form $\hat{H}=-{\rm{i}}{\sum }_{n=1}^{N}{\epsilon }_{n}{b}_{2n-1}{b}_{2n},$ that is ${{O}}_{0}^{T}{{HO}}_{0}={\bigoplus }_{n=1}^{N}\left(\begin{array}{cc}0 & -\,{\rm{i}}{\epsilon }_{n}/2\\ {\rm{i}}{\epsilon }_{n}/2 & 0\end{array}\right)$, with ${\epsilon }_{n}\geqslant 0$. The ground state of $\hat{H}$ is a vacuum state annihilated by all annihilation operators in that base $({b}_{2n-1}-{{\rm{i}}{b}}_{2n})| 0,{n}_{f}\rangle =0$, and the quasiparticles energies ${\epsilon }_{n}$ are arranged in ascending order. We note that some care is needed for a system with degenerated ground state, ${\epsilon }_{1}=0$ (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 $\hat{P}={\prod }_{i=1}^{N}{\sigma }_{n}^{z}={\prod }_{i=1}^{N}({{\rm{i}}{a}}_{2n-1}{a}_{2n})$ commutes with the Hamiltonian $[\hat{H},\hat{P}]=0$, and the relevant instantaneous gap is calculated as energy of two-quasiparticles excitation,

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{n}_{f}}=2({\epsilon }_{1}+{\epsilon }_{2}).\end{eqnarray} \tag{ B2 }$$

In figure 2 we follow the instantaneous ground state in the given parity subspace and we make sure that $\langle 0,{n}_{f}| \hat{P}| 0,{n}_{f}\rangle$ is fixed—and in our case equal 1. In case the true (numerical) ground state is a state with parity −1, we fix it be exciting one-quasiparticle which corresponds to ${\epsilon }_{1}\to -{\epsilon }_{1}$, ${b}_{2}\to -{b}_{2}$ and ${[{{O}}_{0}]}_{n,2}\to -{[{{O}}_{0}]}_{n,2}$. The mixing term ${\rm{\Omega }}({n}_{f})$ in equation (10) is calculated from $\langle 0,{n}_{f}| {\sigma }_{n}^{z}| 1,{n}_{f}\rangle =\langle 0,{n}_{f}| {{\rm{i}}{a}}_{2n-1}{a}_{2n}{d}_{1}^{\dagger }{d}_{2}^{\dagger }| 0,{n}_{f}\rangle$ using the Wick's theorem, where ${d}_{n}^{\dagger }=({b}_{2n-1}+{{\rm{i}}{b}}_{2n})/2$ creates the quasiparticle with energy ${\epsilon }_{n}$.

The time evolution in section 4 is simulated in the Heisenberg picture $\tfrac{\partial }{\partial t}{a}_{n}(t)={\rm{i}}[\hat{H}(t),{a}_{n}(t)]$. For a free fermionic problem like in equation (B1) time-dependent operators can be expanded in the base of original Majorana fermions, $\vec{a}(t)={O}(t)\vec{b}$. This leads to the time-dependent Bogoliubov equations

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{O}(t)=-4{\rm{i}}{H}(t){O}(t),\end{eqnarray} \tag{ B3 }$$

with the initial condition ${O}(t=0)={{O}}_{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 $\exp (-{\rm{i}}{\rm{d}}{t}{{\rm{H}}}_{J(g)})$ that can be efficiently calculated without diagonalizing the full Hamiltonian at each time step. Finally, the energy of the final state (here g_f = 0) is found as $-{\sum }_{n=1}^{N-1}{J}_{n}\langle {\sigma }_{n}^{x}{\sigma }_{n+1}^{x}{\rangle }_{t={t}_{{\rm{final}}}}=-{\sum }_{n=1}^{N-1}{J}_{n}{\rm{i}}\langle {a}_{2n}({t}_{{\rm{final}}}){a}_{2n+1}({t}_{{\rm{final}}})\rangle$.