Engineering non-equilibrium quantum phase transitions via causally gapped Hamiltonians

The inhomogeneity slope and its horizontal and vertical velocities of inhomogeneity can be characterized by a set of hyper-parameters $\{\alpha ,{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}_{{fix}},t)$ as:

$$\begin{eqnarray}&&{v}^{v}(n,t)=-\partial g(n,t)/\partial t,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{v}_{k}^{h}(n,t)=\partial n({g}_{{fix}},t)/\partial t,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\alpha (n,t)=\partial g(n,t)/\partial n.\end{eqnarray} \tag{ 9 }$$

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 IQA with respect to possible temporal and/or spatial inhomogeneities.

In this class, we consider a general form of independent or separated space and time quantum fluctuations to drive the annealing dynamics

$$\begin{eqnarray}&&g(n,t)=\displaystyle \sum _{l=0}{h}_{l}(n){\lambda }_{l}(t).\end{eqnarray} \tag{ 10 }$$

An example of this class will be a periodic spatial inhomogeneity (standing wave) combined with spatially-independent time-local inhomogeneity as:

$$\begin{eqnarray*}&&g(n,t)={h}_{0}{\lambda }_{0}(t)+{\lambda }_{1}(t)\displaystyle \sum _{k=-\infty }^{\infty }{a}_{k}{{\rm{e}}}^{{\rm{i}}\pi {kn}},\end{eqnarray*}$$

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.

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 [{\theta }_{k}(| | n-{n}_{k}| | -{v}_{k}(n)t)]:$

$$\begin{eqnarray}&&g(n,t)={g}_{c}\left\{1+\displaystyle \sum _{k=1}^{M}{\omega }_{k}\tanh [{\theta }_{k}(| | n-{n}_{k}| | -{v}_{k}(n)t)]\right\},\end{eqnarray} \tag{ 11 }$$

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 $\tanh [{\theta }_{k}(n-{v}_{k}(n)t)]\,\simeq$ ${\theta }_{k}(n-{v}_{k}(n)t)$, then for each cluster we get:

$$\begin{eqnarray}&&{g}_{k}(n,t)={g}_{c}\{1+{\theta }_{k}(n-{v}_{k}(n)t)]\}.\end{eqnarray} \tag{ 12 }$$

In the first example of this type, we consider an inhomogeneity with constant v_k for each cluster of the form ${g}_{k}(n,t)={g}_{c}\{1+\tanh [{\theta }_{k}(n-{v}_{k}t)]\}$ which yields $n={\tanh }^{-1}({g}_{k-{\rm{fix}}}/{g}_{c}-1)/{\theta }_{k}+{v}_{k}t.$

$$\begin{eqnarray*}\begin{array}{rcl}{v}^{v}(n,t) & = & {g}_{c}[1-{\tanh }^{2}[{\theta }_{k}(n-{v}_{k}t)]]{\theta }_{k}{v}_{k}={\alpha }_{k}(n,t){v}_{k},\\ {v}_{k}^{h}(n,t) & = & {v}_{k},\\ {\alpha }_{k}(n,t) & = & {g}_{c}[1-{\tanh }^{2}[{\theta }_{k}(n-{v}_{k}t)]]{\theta }_{k}.\end{array}\end{eqnarray*}$$

In the next example of this type, we consider linear approximation of activation function near critical point which yields ${g}_{k}(n,t)={g}_{c}\{1+{\theta }_{k}(n-{v}_{k}t)]\}$ and $n=({g}_{k-{\rm{fix}}}/{g}_{c}-1)/{\theta }_{k}+{v}_{k}t.$ Thus we have: ${v}^{v}={g}_{c}{\theta }_{k}{v}_{k}={\alpha }_{k}{v}_{k}$, ${v}_{k}^{h}={v}_{k}$, and ${\alpha }_{k}={g}_{c}{\theta }_{k}$.

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)=\mathrm{const}.$, and linear velocity ${\lambda }_{0}(t)=(1-t/{\tau }_{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/{\tau }_{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 $\theta \to 0$ and $v(n)\to \infty$, while $\theta 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/{\tau }_{Q}$, ${v}_{k}^{h}(n,t)=\infty$, and α (n, t) = 0.

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 equation (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 [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 figure 2, in which different snap-shots of time-dependent transverse fields are plotted along the chain. In figure 2(a), we illustrate the trivial/homogeneous schedule. In figure 2(b) we explore the effects of periodic critical fronts by constructing the schedule

$$\begin{eqnarray}&&g(n,t)={g}_{i}(1-t/T)+a\cos ({kn})\sin (\pi t/T),\end{eqnarray} \tag{ 13 }$$

where T is the total evolution time. Finally, in figure 2(c) we illustrate an example of multiple-critical-fronts strategy.

Zoom In Zoom Out Reset image size

In the latter, in order to decide the position of the cluster, we employ a simple preprocessing. It is based on observation, that for strictly 1D geometry, the instantaneous local gap may be set by a single, very weak link. We discuss it in more detailed in section 3. Such weak link sets the local time-scale needed for adiabatic transition dividing the chain into two weakly interacting parts. We want to place the borders between the clusters at such links, as (i) they would require the longest time to align according to the weak coupling and (ii) the energy penalty for placing the defects there is the smallest. To that end, we look for the largest cluster—starting at the end of the chain or at the end of the previous cluster—where ${\min }_{{{\rm{cluster}}}_{k}}({J}_{n,n+1})\cdot \kappa \gt {v}_{k}$. Here, v_k would be the velocity of the front for this (kth) cluster, and κ is a parameter fixing the exact value of the threshold (in practice κ ≈ 2). This condition would allow for adiabatic transition as if the energy gap were set by single links only. If the condition is not satisfy the considered candidate for cluster is cut at its weakest link, creating a new smaller cluster where we check the condition again. The procedure is repeated until full chain is divided into clusters. As a result, 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}_{k}^{v}=(| {g}_{f}-{g}_{i}| +\alpha {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, 2... L (counting from the beginning of the cluster) the magnetic field is constructed as:

$$\begin{eqnarray}\begin{array}{c}g(n,t)=\left\{\begin{array}{ll}{g}_{i}, & \left|n-\displaystyle \frac{L}{2}\right|-{tv}\,\gt \displaystyle \frac{{g}_{i}-{g}_{f}}{2\alpha },\\ \displaystyle \frac{{g}_{i}+{g}_{f}}{2}+\alpha \left(\left|n-\displaystyle \frac{L}{2}\right|-{tv}\right), & | | n-\displaystyle \frac{L}{2}| -{tv}| \leqslant \displaystyle \frac{{g}_{i}-{g}_{f}}{2\alpha },\\ {g}_{f} & \left|n-\displaystyle \frac{L}{2}\right|-{tv}\,\lt \displaystyle \frac{{g}_{f}-{g}_{i}}{2\alpha }.\end{array}\right.\end{array}\end{eqnarray} \tag{ 14 }$$

The probability distribution of topological defects for these protocols are presented in figure 3. The defect density is evaluated for a strongly disordered instance of a spin chain consist of up to about 1000 qubits with ${J}_{n,n+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.

Zoom In Zoom Out Reset image size

The scaling of annealing or quenching time as a function of inverse residual energy is plotted in figure 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 ${\epsilon }_{Q}\sim {(\mathrm{log}T)}^{-\gamma }$. We obtain γ ≈ 3.8 ± 0.4. It can be compared with γ ≈ 3.4 which has been obtained by from Caneva et 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.

Zoom In Zoom Out Reset image size

On the other hand, in type II protocol we observe ${\epsilon }_{Q}\sim {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 time-scales (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 $\sim 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 for ${J}_{n,n+1}$ drawn from a uniform distribution in [−1, 1].

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 as $\hat{v}\sim \hat{\xi }/\hat{t}$. Here, $\hat{\xi }$ and $\hat{t}$ are the effective length scale and time scale that the system experience according to KZM given by the relations (2) and (3), respectively. This yields $\hat{v}={\tau }_{Q}^{\nu (1-z)/(1+z\nu )}$. Next we consider the control parameter ε(n, t) to be position dependent. We linearize the relation at the critical front for fixed position n_fixed as

$$\begin{eqnarray}&&\varepsilon (n,t)=\alpha ({n}_{{\rm{fixed}}}-{vt})=-\alpha {vt}+\mathrm{const}.\end{eqnarray} \tag{ 15 }$$

This gives the local annealing rate ${\tau }_{Q}(n)=1/(\alpha v)$. Causality implies that if $v\gg \hat{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 for $\hat{v}$, which leads to

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

Alternatively, we could look at the instantaneous Hamiltonian resulting from inhomogeneous front in equation (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 as

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

It follows from a variant of KZM argument, so-called KZM in space [46–49], where one asks about characteristic distance from n_fixed 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–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 equation (17) allows to estimate the gap of the instantaneous Hamiltonian as

$$\begin{eqnarray}&&{\hat{{\rm{\Delta }}}}_{i}\sim {\alpha }^{\tfrac{z\nu }{\nu +1}}.\end{eqnarray} \tag{ 18 }$$

By combining those characteristic scales we obtain a threshold velocity v_t which is again given by equation (16). The meaning of this threshold velocity is however different here. Namely, we can expect that if the velocity v in equation (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 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, in figure 5 we illustrate how the defects can be created when two clusters are merging together. This highlights the interplay of the causal effects and different boundary conditions on the defect suppression.

Zoom In Zoom Out Reset image size

As we have discussed in section 1, the disordered systems have a completely different critical phenomena as $z\nu \to \infty$. In this case one has to modify equation (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 equation (17). We elaborate on this in the next section in the context of multiple critical fronts driving strategy.

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 derive 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 site-dependent transverse field which is the largest energy scale within a local neighborhood of our system; that is ${g}_{i}\gt {J}_{{ij}}\forall 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 as ${\tilde{J}}_{{jk}}=\tfrac{{J}_{{ij}}{J}_{{ik}}}{{g}_{i}}$, if J_jk = 0 or otherwise ${\tilde{J}}_{{jk}}=\max \left({J}_{{jk}},\tfrac{{J}_{{ij}}{J}_{{ik}}}{{g}_{i}}\right)$. 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}}\geqslant {J}_{{ik}}\forall k$ and ${J}_{{ij}}\geqslant {J}_{{lj}}\forall 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 with 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 ${\tilde{g}}_{i}=\tfrac{{g}_{i}{g}_{j}}{{J}_{{ij}}}$.

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 figure 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 causal dependence emerging in the multiple critical fronts strategy, see the SDRG visualizations of our protocols in figure 6.

Zoom In Zoom Out Reset image size

Next, we discuss an upper-bound for a global threshold velocity such that the multiple critical fronts 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 figure 6. When the fronts are sufficiently separated, the RG flow makes the corresponding clusters exponentially decoupled. Consequently, when looking at the possible transitions, we can consider each front independently. The transition matrix elements for each cluster k:

$$\begin{eqnarray}&&\left|\langle 0,t| \displaystyle \frac{{\rm{d}}\hat{H}}{{{\rm{d}}{n}}_{k}^{f}}| 1,t\rangle \right|\displaystyle \frac{{{\rm{d}}{n}}_{k}^{f}}{{\rm{d}}{t}}={{\rm{\Omega }}}_{k}({n}_{k}^{f}){v}_{k},\end{eqnarray} \tag{ 19 }$$

where n_k^f encodes the (time dependent) position of the front and ${v}_{k}=\tfrac{{{\rm{d}}{n}}_{k}^{f}}{{\rm{d}}{t}}$ is its vertical velocity. ${{\rm{\Omega }}}_{k}({n}_{k}^{f})$ defines local instantaneous ground state to first excited state transition matrix elements for that front. 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 causal 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, ${{\rm{\Delta }}}_{k}({n}_{k}^{f});$ that is $\tfrac{{{\rm{\Omega }}}_{k}({n}_{k}^{f})}{{{\rm{\Delta }}}_{k}^{2}({n}_{k}^{f})}$. Then the low-energy quasi-adiabatic dynamics is expected for

$$\begin{eqnarray}&&v\ll {v}_{t}=\mathop{\min }\limits_{k}\displaystyle \frac{{{\rm{\Delta }}}_{k}^{2}({n}_{k}^{f})}{{{\rm{\Omega }}}_{k}({n}_{k}^{f})}.\end{eqnarray} \tag{ 20 }$$

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 equation (20) with the threshold velocity which was derived in the previous section. There, it was calculated as $\hat{v}={\hat{\xi }}_{i}/{\hat{t}}_{i}$, where ${\hat{\xi }}_{i}$ and ${\hat{t}}_{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. that the gap appearing in equation (20) is squared comparing to other expression), we note that ${{\rm{\Omega }}}_{k}({n}_{k}^{f})$ and ${{\rm{\Delta }}}_{k}({n}_{k}^{f})$ are not independent. We expect that ${{\rm{\Omega }}}_{k}({n}_{k}^{f})/{{\rm{\Delta }}}_{k}({n}_{k}^{f})\sim {\xi }_{i}^{-1}$, i.e. it is directly proportional to the size of the effective critical region. We show this in figure 7 and discuss further below. Importantly, this allows us to focus on the scaling of the gap.

Zoom In Zoom Out Reset image size

The effective size of the critical region (for a given front) is given by equation (17). In our case ν = 2, which gives ${\xi }_{i}\sim {\alpha }^{-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, ${{\rm{\Delta }}}_{k}({n}_{k}^{f})\sim {{\rm{e}}}^{-\mathrm{const}\sqrt{{\xi }_{{i}}}}\sim {{\rm{e}}}^{-\mathrm{const}\cdot {\alpha }^{-1/3}}$ [37].

To be more precise, we consider P(Δ, α) as a distribution of the minimal relevant instantaneous gap Δ, 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={\alpha }^{1/3}\mathrm{log}{\rm{\Delta }}$, with the expected $P({\rm{\Delta }},\alpha ){\rm{d}}{\rm{\Delta }}=\tilde{P}(x){\rm{d}}{x}$, and universal distribution $\tilde{P}(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({\rm{\Delta }},N){\rm{d}}{\rm{\Delta }}=\tilde{P}(x){\rm{d}}{x}$ with $x={N}^{-1/2}\mathrm{log}{\rm{\Delta }}$ [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 equation (17).

We verify those scaling predictions in figure 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 (which would be related with one of the fronts), which we distinguish by finding the minimal energy eigenstate with the corresponding transition matrix element Ω 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 Δ 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 figure 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\in [-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 $\sim {{\rm{e}}}^{-\mathrm{const}\cdot {\alpha }^{-1/3}}$, which can be attributed to many-body effects, is larger.

The analysis above regards local, instantaneous gap related with the critical fronts traveling though the chain. To quantify the difficulty of the problem, we consider the minimal gap encounter during such quench. To that end we focus on the distribution $P({\rm{\Delta }},\alpha ){\rm{d}}{\rm{\Delta }}=\tilde{P}(x){\rm{d}}{x}$. Importantly $\tilde{P}(x)$, as calculated for the homogeneous critical case in [56, 58], has a Gaussian tail for $| x| \gg 1$,⁴ i.e.

$$\begin{eqnarray}&&\tilde{P}(x)\sim {{\rm{e}}}^{-{{ax}}^{2}}.\end{eqnarray} \tag{ 21 }$$

Figure 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⁵ . Now, let the probability that Δ is smaller than some ${{\rm{\Delta }}}_{\min }^{q}$ be

$$\begin{eqnarray}&&P({\rm{\Delta }}\lt {{\rm{\Delta }}}_{\min }^{q},\alpha )={\int }_{-\infty }^{{x}_{q}}\tilde{P}(x){\rm{d}}{x}=\epsilon ,\end{eqnarray} \tag{ 22 }$$

where ${x}_{q}=\mathrm{log}({{\rm{\Delta }}}_{\min }^{q}){\alpha }^{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

$$\begin{eqnarray}&&P({{\rm{\Delta }}}_{\min }\gt {{\rm{\Delta }}}_{\min }^{q})={(1-\epsilon )}^{N}=q,\end{eqnarray} \tag{ 23 }$$

where Δ_min is the minimum from the sample of N instantaneous gaps. This equation defines ${{\rm{\Delta }}}_{\min }^{q}$ as a q-quantile for the global minimal gap. Now, we obtain ${{\rm{\Delta }}}_{\min }^{q}$ from the above equation. This gives $\tfrac{1}{N}\mathrm{log}q=\mathrm{log}(1-\varepsilon )\simeq -\epsilon \sim \tfrac{1}{2{{ax}}_{q}}{{\rm{e}}}^{-{{ax}}_{q}^{2}}$, which is obtained by expanding the error function for the Gaussian tail in equation (21), to the leading order in small x_q. Solving this equation in the leading order we obtain

$$\begin{eqnarray}&&{\alpha }^{1/3}\mathrm{log}{{\rm{\Delta }}}_{\min }^{q}={x}_{q}\approx -\sqrt{\displaystyle \frac{1}{2a}\mathrm{log}\displaystyle \frac{N}{2a\mathrm{log}{q}^{-1}}}.\end{eqnarray} \tag{ 24 }$$

This suggests that if we fix a quantile q, then in the asymptotic limit

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\min }^{q}\sim \exp (-\mathrm{const}\cdot {\alpha }^{-1/3}\cdot \sqrt{\mathrm{log}N}),\end{eqnarray} \tag{ 25 }$$

i.e. it is vanishing slower than any polynomial with increasing N. This has to be compared with homogeneous gap scaling as: ${{\rm{\Delta }}}_{\min }^{q}\sim \exp (-\mathrm{const}\cdot \sqrt{N})$, which vanishes as stretched exponential with N. We note here that similar analysis in case when $\tilde{P}(x)$ would have exponential tail for large negative x would give polynomial dependence for the minimal gap on N.

In order to fully analyze the tail of the gap distribution we also consider the situation that the minimum gap is enforced by very small local link rather than many-body 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 quasi-adiabatically. To that end, let us assume that the links are drawn from uniform distribution ${J}_{i}\in [-1,1]$. Probability that a single link is weaker then some is equal , or $P(| {J}_{i}| \gt \epsilon )=1-\epsilon$. Let us consider that q is the probability that all the links are stronger than this is $q={(1-\epsilon )}^{N}$. For small this yields

$$\begin{eqnarray*}&&\epsilon \simeq -\displaystyle \frac{\mathrm{log}q}{N}.\end{eqnarray*}$$

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 figure 7 for larger values of α.

Finally, in figure 8 we illustrate the relation between the transition matrix elements Ω and the gap Δ. To that end, for simplicity we consider protocol with single front and, as in figure 7, sample the values of minimal relevant gap and corresponding Ω 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={\alpha }^{1/3}\mathrm{log}({\rm{\Delta }})$ as above, and $y=\mathrm{log}({\alpha }^{-2/3}{\rm{\Omega }}/{\rm{\Delta }})$ to reflect the expected relation ${\rm{\Omega }}/{\rm{\Delta }}\sim {\xi }_{i}^{-1}\sim {\alpha }^{2/3}$. We plot the obtained distributions of P(x, y) in figure 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 Ω/Δ—i.e. where it is most relevant—coincide with the local minima of energy gap. Apart from those point, Ω/Δ is quickly approaching zero, which is reflected by elongated shape of the distribution P(x, y) in the direction of small (irrelevant) y.

Zoom In Zoom Out Reset image size