Multiple-scale integro-differential perturbation method for generic non-Markovian environments

Non-Markovianity may significantly speed up quantum dynamics when the system interacts strongly with an infinite large reservoir, of which the coupling spectrum should be fine-tuned. The potential benefits are evident in many dynamics schemes, especially the continuous-time quantum walk. Difficulty exists, however, in producing closed-form solutions with controllable accuracy against the complexity of memory kernels. Here, we introduce a new multiple-scale perturbation method that works on integro-differential equations for general study of memory effects in dynamical systems. We propose an open-system model in which a continuous-time quantum walk is enclosed in a non-Markovian reservoir, that naturally corresponds to an error correction algorithm scheme. By applying the multiple-scale method we show how emergence of different time scales is related to transition of system dynamics into the non-Markovian regime. We find that up to two long-term modes and two short-term modes exist in regular networks, limited by their intrinsic symmetries. In addition to the effective approximation by our perturbation method on general forms of reservoirs, the speed-up of quantum walks assisted by non-Markovianity is also confirmed, revealing the advantage of reservoir engineering in designing time-sensitive quantum algorithms.


I. INTRODUCTION
Many mathematical techniques have been developed to identify intrinsic time scales in dynamical systems. Multiple-scale perturbation [1], for example, is a highlydeveloped approximation method that solves complex dynamics in a perturbative way, by introducing trial variables as different time scales which are often of physical importance themselves. It is used in many fields, especially in such quantum optics topics [2] as spontaneous radiation processes [3][4][5] and nonlinear solitons [6,7].
A summary of its applications in quantum optical problems can be found in Ref. [8]. Also originated from quantum optics is the study of open quantum systems [9], in which the dynamics of a quantum system is enriched when the system is open to an infinitely large reservoir from which microscopic interaction occurs. In particular, when the interaction spectrum is fine-tuned the open system may go through a quantum phase transition [10] and exhibit such non-Markovian features as bidirectional exchange of energy and coherence between the system and the reservoir, accompanied by singularity in time-dependent system variables [11]. Assisted by non-Markovianity, unexpected atypical dissipation and decoherence behaviors [12] can arise in open systems, e.g., the sudden death of entanglement [13,14]. It was recently found that non-Markovianity changes quantum speed limits [15,16] and can significantly speed up quantum dynamics [17]. * xm@bu.edu † hes@bu.edu The benefits of properly utilizing and engineering non-Markovian reservoirs [18,19] for designing time-sensitive quantum algorithms and protocols are obvious.
In this paper, we examine one specific quantum dynamics scheme, the continuous-time quantum walk, of which the concept was first introduced in Ref. [20]. As the name suggests, a walker moves continuously in time, navigating among different sites. Unlike a classical random walker, the propagation of the quantum walker is coherent, i.e., besides the randomness inherited from quantum-mechanical probability amplitudes the coherence between different sites also governs the system dynamics [21]. The study of continuous-time quantum walks is mathematically based on network theory [22][23][24] and is closely related to other fields in quantum information theory, e.g., universal quantum computation [25], quantum algorithms [26][27][28], and perfect state transfer [29]. A continuous-time quantum walk can be experimentally implemented, usually by such quantum optical systems as waveguides [30] or Rydberg atoms [31].
One of the breakthroughs in quantum information theory is the finding of a quadratical speed-up (∼ O( √ N )) in a quantum search algorithm, namely, Grover's algorithm [32]. Grover's algorithm has proven to be equivalent to a continuous-time quantum walk on a complete network [26]. In addition, a speed-up by continuous-time quantum walks on a star network has also been found [27], and this has encouraged further study of continuous-time quantum walks on regular networks. Here, a network is said to be regular if it has a repeating pattern in its network topology. This should not be confused with k-regular graphs which are defined differently. To expand the theoretical structure, research on open-system quantum walks has also been conducted [33,34] in order to deal with noisy environments [35,36], but exact solutions or model-based descriptions are rare, especially in the non-Markovian regime where the system dynamics is inevitably governed by strong memory effects. On the other hand, we expect that many interesting features granted by non-Markovianity should also exist in the quantum walk scheme, among which the most useful is the steep decrease of quantum speed limit [17].
To construct a practical methodology and verify the non-Markovianity-assisted speed-up, we propose a generalized version of the multiple-scale perturbation method that now can work on integro-differential equations. Strong memory effects can be studied by perturbatively expanding the memory kernel. The advantage of using a perturbation method is being able to derive closed-form approximate solutions where accuracy and complexity are controllable. We further apply the perturbation method to a continuous-time quantum walk enclosed in a non-Markovian reservoir. Such a model naturally follows an error correction algorithm scheme, with the reservoir a collection of independent "error" sites. We find that two time scales of different physical importance emerge when the studied system moves into the non-Markovian regime. With the accuracy of our perturbation method being guaranteed, we investigate how quantum-walk dynamics is affected by the coupling strength between system and reservoir, as well as by the intrinsic network topology. The expected speed-up is confirmed by looking into the four different eigenfrequencies hidden in the non-Markovian dynamics, inherited from symmetries of regular networks.
It is unclear yet whether this speed-up can be utilized as a quantum resource, given that non-Markovianity can be detrimental for certain delicate quantum tasks [37]. We expect that non-Markovianity should be the most useful for tasks that are speed-focused, e.g., quantum simulation [38]. Our method is suitable to do a more general study on this issue.
The rest of this paper is organized as follows. Section II introduces the multiple-scale perturbation method and its generalization to integro-differential equations, which we expect to have broader applications than quantum walks. Section III introduces the concept of continuoustime quantum walks and the quantum algorithm picture behind the interaction with non-Markovian reservoirs. Section IV presents a test of the performance of our multiple-scale approximation method for general reservoirs and examines the relation between time scales and non-Markovianity. Section V describes continuoustime quantum walks on some regular networks, including complete networks, star networks, rings, and square lattices, and reveals the existence of up to four system modes in terms of the two time scales and how they connect to network topology.
Remaining questions in this crossover study of non-Markovian memory effect and continuous-time quantum walks include the possibility that there are more than two independently important time scales in the non-Markovian model. We hope to understand these time scales better in a systematic manner in the future. We could also apply our method to quantum walks on complex networks [39] where the statistics of disorder could take unpredictable new forms.

II. MULTIPLE-SCALE INTEGRO-DIFFERENTIAL PERTURBATION METHOD
Often in dynamical systems there is no exact solution y(x) to the system dynamics. It is difficult to acquire accurate and reliable approximations. The regular perturbation method [1] is a powerful approach in finding an approximation of the unknown dynamics. By introducing a small dimensionless perturbation parameter α, any y(x) can be expressed in terms of a power series in α and can be well approximated in a closed form by finite leading terms. In practice, however, given an intricate system, the complex behavior of y(x) often invalidates the perturbation approximation. In the power series, each term follows the same approximate over-simplified dynamics, causing the complexity of y(x) to be limited. Thus the complex behavior is unavoidably lost when the infinite series is broken down.
The multiple-scale perturbation method [1] has then been introduced to overcome this difficulty. The basic idea is to add more degrees of freedom in terms of new independent variables into the system. Different characteristics are captured by different variables, an attempt to retrieve the complexity of y(x) in each single perturbation term. The procedure of the multiple-scale perturbation method goes as follows: first, two (or more) different and artificial scales are chosen as functions of x and α; the scales are considered as new independent variables, allowing the ordinary differential equations to be converted into partial differential equations and expanded into perturbation series.
The additional degree of freedom from the extra variables is eventually constrained by requiring that higher perturbation terms diverge no more quickly than lower perturbation terms [1].
The classical Duffing equation is a neat example appearing in the study of anharmonic oscillators, which cannot be solved by a regular perturbation method due to its small nonlinear term [8]. In fact, the secular term in the first-order perturbation correction y (1) (x) is proportional to x sin (x), which is unbounded and incorrect [8]. Now, by choosing two independent scales, u = x and v = αx, one has, locally, in the neighborhood of x. A term-by-term perturbation expansion yields up to the first two orders. Thus one has y (0) (u, v) = h(v) cos u+k(v) sin u. h(v) and k(v) are determined next by requiring −4 y (0) 3 − 2 ∂ 2 ∂u∂v y (0) = 0 so that y (1) does not diverge. With initial conditions y(0) = 1 and y (0) = 0 the final result reads which is not divergent and exhibits the first-order correction to the frequency [8].
Based on the same thought, we generalize the method to an integro-differential equation which is of a general form Here, F(d/dx, y(x)) represents an arbitrary differential term(s) and G(x − x ) is a convolution kernel. The difficulty arises when an integral term is added, that the convolution is not a local operation and thus the artificial scales cannot be considered independent. The integral term has to be dealt with indirectly. If the kernel is holomorphic near x = x , then G(x−x ) can be expanded as a series of (x − x ) living in its neighborhood without any singularity. p and q are understood as integers, w.l.o.g. The expansion further suggests that locality can be regained from Eq. (1) if the perturbation procedure is fine-tuned. Here, the trick is to bring Eq. (1) into higher and higher differential orders, meanwhile trying to cancel out integral terms or at least make them comparably smaller than differential terms. The iterative procedure of the multiple-scale perturbation method thus contains three steps for each perturbation order: 1. Let d/dx act on both sides of the integrodifferential equation to solve.
Recall that (d/dx) n+1 dx (x − x ) n y(x ) = n!y(x), a guarantee of being locally holomorphic near x.
2. Introduce artificial scales and replace d/dx by partial differential operators in terms of the new variables.
3. Extract the lowest order terms. They are the corresponding perturbation correction to be calculated. The rest terms are saved for next iteration.
During the iteration, more and more integral terms are differentiated. With a successful choice of scales and α, all perturbation corrections should only consist of local terms, of which the solutions are easily carried out by the usual multiple-scale approach [1].
Note that perturbation methods enable us to understand the functional importance of different dynamic terms and parameters. The solutions generated by perturbation methods are also in a closed form that is more suitable for practical purposes, in part because the precision of numerical calculations in closed-form functions are more controllable than in complex functions. The study of integro-differential equations has traditionally involved the use of integral transforms, but few solutions have closed-form expressions. We will see that the regular perturbation method also cannot be used if we are to expand the solutions perturbatively in the transformed domain. The different behaviors of intrinsic scales in the system simply cannot be captured, the perturbation corrections remaining unbounded or erratic. Hence the use of the multiple-scale perturbation method is essential to meet the exact needs.

III. CONTINUOUS-TIME QUANTUM WALKS
The simplest definition of a continuous-time quantum walk (CTQW) involves a connected network G(V, E) which comprises a set V of nodes and a set E of links. The weighted adjacency matrix A of G is an Hermitian The walker is a N -dimensional complex vector c(t) which follows the dynamicsċ Sometimes, instead of an adjacency matrix, a Laplacian matrix is preferred, yet for both the dynamics is equivalent on regular networks [40]. The first systematic study of Eq. (3) from a statistical physics perspective dates back to the Anderson localization model [41], where the diffusion of a single electron that scatters in-between N sites, if described well by short-range interaction and tight-binding approximation, can be reduced to simpler dynamics in the form of Eq. (3). The diffusion behavior in terms of the electron state |ψ(t) = i c i (t) |i is determined by the amount of "impurity" in the Hamiltonian H. Between any two sites i and j the interaction i| H |j is simply A ij .
Another perspective-which is more general-is to consider a spin system with Hamiltonian where the occupation number i σ i + σ i − in the Fock space is manifestly conserved, and view the continuous-time quantum walk as a spin diffusion in the one-exciton Hilbert subspace H sub = C N . The operators σ ± are Pauli operators. We note that statistical studies of such a spin system beyond one exciton have undergone difficulties and eventually led to the theory of many-body localization [42], which is however not our focus here.  [9]. The whole system can also be viewed as the union of two networks: an independent network G and a complete bipartite network

A. Quantum walk in a reservoir
The solution of Eq. (3) is relatively simple and has been studied well even in the N → ∞ limit [43]. However, by introducing a large reservoir that is described by a free Hamiltonian with infinite modes, H E = ∞ k=1 ε k a † k a k , the system dynamics will be dramatically changed if further undergoing an interaction term, where w i g k is the coupling between spin i and mode k, and a k and a † k are annihilation and creation operators of either bosons or fermions. The subspace H sub = C N ⊕ C ∞ is still independent though, conserved by i σ i where it does not matter whether a k and a † k follow commutation or anti-commutation relations. The quantum walker state can be written as is described by open system dynamics [9]. After some calculations we derive, in the matrix formulation, or simplyċ + iAc = −G * c where * is the convolution operation. We have and G(t) = w (w * ) G(t), the memory kernel G being The Fourier transform of G(t) is called the spectral density J(ω) of the reservoir [9]. When J(ω) is close to a constant, the spectrum is white-noise like, yielding a memory-less kernel G(t) ∼ δ(t). The open system is Markovian and is only subjected to exponential decay. However, predicted by the open quantum system theory, a "colored" spectrum will cause non-Markovian effect, and in the extreme case J(ω) ∼ δ(ω) the system undertakes to-and-fro oscillation without sign of dissipation. In general, G(t) does not have to be a single form. If each element G ij (t) has independent spectrum density, then the system behavior will be even more complex.
Thus different forms of J(ω) suggest completely different behaviors [44]. We therefore assume that there is a universal dimensionless factor α which compares the bandwidth of J(ω) to its peak. For multimodal distributions α is interpreted as the average for different modes. In the weak-coupling regime, α 1, the Markovian open system is well described by the Born-Markov approximation [12], a regular perturbation approach with the perturbation parameter being α −1 . In the strong-coupling regime, α 1, the non-Markovian system behaves more interestingly, but the regular perturbation method no longer works [9]. Note that this model can be easily realized in cavity QED [45], where H G describes a system of two-level dipoles, H E describes a quantized radiation bath with a k and a † k being electromagnetic fields, and w i g k represents the strength of coupling between dipole i and cavity mode k. The quality factor of cavity Q is important [2] and can be related to the perturbation parameter by α ∼ Q −1/2 . It is of practical interest as well, that all results in this paper can actually be tested through cavity-QED experiments.
Back to network theory, as shown in Fig. 1, the open system and reservoir can be considered together as the union of G(V, E) and a complete bipartite network here is bipartite, i.e., it has two disjoint and independent set of nodes, V and V E , and has links of which each only connects one node from V to one node from V E . The symbol ∀ denotes that the bipartite network actually contains every link from V to V E and thus is complete.
It is a conceptually useful point of view, to regard all nodes that E consists of as "error" sites. By walking into an "error" site, the walker realizes its mistake and attempts to go back to the possible "target" sites in G, but never to enter another "error" site. Therefore, for any quantum search algorithm in terms of CTQW, the dynamics here naturally follows an error correction scheme that functions on a solution space G plus an "error" space E. Within this scheme, it is clear that the error correction is carried out by non-Markovian feedback controls, the efficiency of which is dependent not only on the "error" rate w but also on the spectrum, J(ω). This leads us to better understanding of error correction algorithm design.

B. Perturbation expansion
The last step is to apply the multiple-scale integrodifferential perturbation method to Eq. (4).
We introduce three new dimensionless quantities,t = γt, Here γ ∼ i,k w i g k has a dimension of inverse time. We choose two independent time scales, a primary scale T and an auxiliary scale τ , written as (6) and also,T = γT andτ = γτ . The coefficients {A n } and {B n } are to be determined. We assume thatG(t) is holomorphic neart, derived from Eq. (4), with c(t) = m α m c (m) (t) and I the identity matrix. Note thatT ,τ , andt appear simultaneously in Eq. (7), a result of mixing local partial derivatives and nonlocal integrals together.
Following the perturbation procedure, we let d/dt act on both sides of Eq. (7) and extract the terms of the lowest perturbation order, The solution of Eq. (8) is derived by letting d/dt act again on the rest terms and then extracting the new lowest-order ones. Note that c (1) terms should not dominate c (0) terms, so all c (0) must cancel out with each other, which further implies sτ ), which are to be fixed by initial conditions. Γ Higher-order perturbation corrections are subjected to the same procedure, where {A n } and {B n } are also going to be fixed, an example of which is given in the next section.

IV. NON-MARKOVIAN RESERVOIRS
In this section, we study N = 1 systems and investigate different non-Markovian reservoirs in order to understand how non-Markovianity changes system behavior and invalidates the regular perturbation method.

A. Example: Lorentzian reservoir
A Lorentzian reservoir is one of the few types that can be solved exactly in closed form, which is the reason why it is used here as a preliminary example. The spectral density is where λ is the spectral width, γ is the coupling strength, and ∆ is the off-resonance frequency. The corresponding memory kernel is

Comparison of perturbation methods
We only consider the simplest configuration that ∆ = 0 so the system is on resonance. The solution of Eq. (4) is simply where D = 2γλ − λ 2 , given initial conditions c(0) = 1 andċ(0) = 0 [9]. The parameter D distinguishes the two different coupling regimes: an imaginary D corresponds to the weak-coupling regime where γ < λ/2; a real D corresponds to the strong-coupling regime where γ > λ/2. It thus makes sense to choose α = λ/γ as the universal perturbation parameter. Figure 2 shows approximate solutions of the state population ρ = |c(t)| 2 by the multiple-scale (MS) method compared with the regular perturbation approach that works on the s domain of LT (see Appendix A for details). In Fig. 2, the perturbation solutions of positive orders, LT2 and LT4, as a direct attempt to capture the non-Markovian behavior in the strong-coupling regime, fail to remain bounded. The secular terms in LT2 and LT4 all diverge. In fact, as α → 0, the leading term α p in the memory kernel, Eq. (2), introduces a small boundary layer att ∼ α p , where c(t) may vary so fast that any finite order perturbation term cannot penetrate through the layer and is forced to diverge [1]. LT(-0) and LT(-2), on the other hand, are approximations derived from inversely perturbing the dynamics, i.e., regular perturbation around α −1 . LT(-0) shows an exponential decay and should have done well in the weak-coupling regime. When α −1 1, though, the perturbation series does not converge for the apparent reason.
Instead, MS0 and MS1 are in good agreement with the exact solution. The two closed-form approximations are given by c (0) 2 and c (0) + αc (1) 2 , where i.e., ω and the series expansion ofG(t −t ) yields G 0 = 1/2, The first order perturbation c (1) (T ,τ ) follows the exact form of the zeroth order, and ω (1) ± , while the perturbation correction onT andτ is given by A 1 /A 0 = −1/4 and B 1 /B 0 = 0, respectively. Details of calculation are given in Appendix B.

Non-Markovianity
The rebound of population shown in Fig. 2 is one of the features of non-Markovianity. In general, a quantum Markovian evolution is defined by a set of tracereserving linear maps {E(t, t ), t ≥ t }, where E(t, t ) is the propagator for the open system and is required to be completely positive under the composition law [11] This requirement further gives rise to the Gorini-Kossakowski-Susarshan-Lindblad theorem [11]: a quantum evolution is Markovian if and only if D(t) ≥ 0 for all t, where D(t) is named the dissipator in the density matrix formulation [11]. Given the series expansion [Eq. (6)], D(t) has the following form,    Another feature of non-Markovianity related to time scales is the steep decrease of quantum speed limit in the strong-coupling regime [17]. It is found that the evolution time between two orthogonal pure or mixed states may not be unique if the evolution is non-Markovian [17]. In Eq. (12), there are infinite numbers of singularities along the time axis, set by c(t) = 0. They correspond to all possibilities of the evolution time between the two eigenstates of H G . The first singularity is approximately att MS0 = π(2γλ) −1/2 , ort MS1 = (γλ/2) −1/2 (1 − λ/4γ) −1 arccos λ/ (2γ + λ), derived by MS0 and MS1, respectively. Their relative errors are only to the order of α and α 3 , compared to the exact resultt = (2/D)[π − arctan(D/λ)] which is known as the minimal evolution time [17] between the two orthogonal eigenstates.

B. General reservoirs
It is known that reservoir engineering helps producing squeezed states [19] or generating non-Markovianity [18]. A general and well-behaved approximation method is potentially useful for this practical purpose. Here, we apply the MS perturbation method on more general reservoirs and investigate how the goodness of approximation is related to the form of memory kernel.
As shown in Figs. 3(a)-3(c), we choose three different simple but nontrivial memory kernels, where there are still two parameters, γ and λ, by which α ∼ (λ/γ) 1/2 is to be fully constructed. Their forms are similar to the Lorentzian reservoir, but their series expansions in terms of {G n } are clearly different. We define the characteristic times, T c = γ −1 G To be more specific, in Figs. 3(a) and 3(d), the memory kernel is of the form of a Gauss error function, G(t − t ) = (π 3/2 γλ/ √ 2)erfc(λ |t − t | / √ 2). As expected, the approximate solutions MS0 and MS1 are good enough, because erfc(x) ∼ O(x −1 exp(−x 2 )) descends fast enough when x → ∞, making G(t − t ) more dominated by local behavior. In Figs. 3(b) and 3(e), G(t − t ) = γλ/(1 + λ |t − t |). A descending speed of O(x −1 ) makes G(t − t ) less dominated by local behavior and thus implies the approximation to be less accurate. Note that O(x −1 ) also induces logarithmic singularity in the corresponding spectral density. In Figs. 3(c) and 3(f), it is a special case that we have τ c → ∞ implied by the Gaussian memory kernel, G(t − t ) = √ 2πγλ exp(−λ 2 |t − t | 2 /2). The auxiliary scale τ collapses to an infinitesimal point and thus becomes trivial. Given that there is only one attractor at c(∞) = 0 [44], we believe another time scale that corresponds to decay behavior (which, even though existing, must be comparably small) should be hidden in the system with a complex nonlinear manner. This raises the question of how to identify all time scales in a general reservoir, which is theoretically important yet a nontrivial task.

V. REGULAR NETWORKS
In this section, we investigate how network topology determines the propagation of an open-system walker. We study regular networks, e.g., complete networks, star networks, rings, and square lattices, where repeating patterns exist in their network topologies. All results given hereafter are calculated by MS0, i.e., only to the zeroth order. We omit the superscript, ω (0) ≡ ω, and let A 0 = B 0 = 1, w.l.o.g., for simplification. Only the two leading terms G 0 and G 1 from {G n } are used, and we set them to be G 0 = 1/2 and G 1 = −1/2 which coincide with those of the Lorentzian reservoir. The reservoir itself, however, does not have to be of the exact same type.

A. Binary quantum walk
The simplest quantum walk is a binary system, i.e., N = 2. The walker simply chooses between "Yes" and "No", the two nodes of G. Without any reservoir, the binary system is a simple harmonic oscillator, with a frequency of |ω ± 0 | = α −1 A 12 ∼ T −1 . Here, we assume that there are two eigenfrequencies, ω + 0 and ω − 0 , as we will see how the degeneracy is broken by introducing a reservoir. Figure 4 shows how the state populations ρ 1 = |c 1 | 2 and ρ 2 = |c 2 | 2 change over time and how they are compared with the numerical solution in a Lorentzian reservoir. Besides the goodness of approximation, we note that the oscillation frequency between ρ 1 and ρ 2 is indeed significantly increased, as (8.34/γ) −1 π ≈ 0.377γ 0.08γ, a consequence of the non-Markovianityassisted steep decrease of quantum speed limit [17]. By solving Eq. (8), we actually find three eigenfrequencies, ω + 0 ≈ −2.32γ, ω − 0 ≈ 2.17γ, and |ω ± ∞ | ≈ 0.151γ, which satisfy and ω − 0 are larger than α −1 A 12 ≈ 0.253γ, while |ω ± ∞ | is smaller. We use the notations ω ± 0 and ω ± ∞ for the sake of consistency with the following context where we will see that at most four eigenfrequencies exist, regardless of the coupling ratio of individual nodes.
One way to implement a binary quantum walk in cavity QED is to produce a dipole-dipole interaction [46] between two atoms. Oscillation between the two upper levels of the atoms is regarded as the quantum walk between two nodes. Our perturbation method provides a convenient way to study exchange of states and its speed limit in relevant dipole-dipole interaction experiments.

B. Complete network
When N > 2, network topology starts to have an effect on the system dynamics. We let G to be a complete network, so that Here, κ is the link weight which has a dimension of γ. Figure 5 shows the magnitudes of ω ± 0 and ω ± ∞ [by Eq. (8)] and their corresponding decay rates Γ ± 0 and Γ ± ∞ [by Eq. (9)] as functions of κ. It is worth noting that while Γ + 0 and Γ − 0 are bounded, Γ ± ∞ diverge when κ → 0, which is unphysical. To derive a physical solution, we are forced to reassign Γ ± ∞ to another universal and trivial solution of Eq. (9) which we ignored first, namely, Γ ± ∞ ≡ 0. The deeper implication of abandoning the diverging Γ ± ∞ is, as first revealed in Fig. 3(f), that there should be another global and complex time scale(s) which is responsible for long-term decay. Therefore, the approximate solution of the quantum walker state reads    In fact, by numerical solution, we find that those terms related to the long-term frequencies ω ± ∞ indeed decay much slower, the magnitudes of which can be approximated as constants C ± ∞ in the short run. As κ → ∞, Figs. 5(a)-5(c) together show that ω + 0 is asymptotically close to (N − 1)κ/α (dashed line), which is the intrinsic mode frequency of a complete network. At the same time, Γ + 0 → 0 guarantees that the intrinsic mode frequency is secular in the κ → ∞ limit, and thus we recover a CTQW as a closed system.
When κ → 0, ω + 0 and ω − 0 approach a constant, N 1/2 G 1/2 0 , which is the oscillation frequency induced by the reservoir. Since ω ± 0 are directly relevant to non-Markovianity, we deduce that the existence of short-term frequencies is responsible for improvement of quantum speed limit in the non-Markovian regime. Such improvement, however, is only temporary and is controlled by the short-term decay rates Γ ± 0 → −(2G 0 ) −1 G 1 , because in the long run only long-term frequencies |ω ± ∞ | persist which are smaller than (N − 1)κ/α, indicating a slowdown of dynamics.
We also see that the degeneracy between ω + ∞ and ω − ∞ , as well as Γ + ∞ and Γ − ∞ , is broken by introducing small fluctuation on w, as shown in Figs. 5(c) and 5(f). We expected more eigenfrequencies to be split from the eigenspectrum of Eq. (8), but only see four frequencies at most, no matter how large N is. The hidden degeneracies are thus maintained by the topological symmetry of G. Note that although w is randomly realized in Figs. 5(c) and 5(f), different realizations of w indeed yield similar results, since the ensemble average of different realizations actually converges to the one-shot observation when N becomes large, which is guaranteed by the assumed ergodicity of random matrix theory.
Finally, we note that in spite of its speed-up effect on quantum walks, open-system dynamics usually forbids perfect state transfer [29] in the studied system. Hence, introducing a reservoir brings both advantages and disadvantages, between which the trade-off needs to be taken care of by delicate reservoir engineering.

C. Star network
Instead, we let G to be a star network, i.e., All results are presented in Fig. 6 which shares a lot of similarity with Fig. 5. Note however that the difference    between a star network and a complete network is fundamental, as there is always proximate degeneracy between ω + 0 and ω − 0 , as well as Γ + 0 and Γ − 0 (Fig. 6). In particular, comparing Figs. 6(a) and 6(d) with Figs. 5(a) and 5(d) demonstrates how different network topologies affect the emergence of a new eigenvalue in their corresponding spectrums when N = 3. Also, as shown in Figs. 6(a)-6(c), the intrinsic mode frequency of a star network is √ N − 1κ/α (dashed line), which is smaller than that of a complete network. |ω ± ∞ | ≤ √ N − 1κ/α still holds true. Nevertheless, we see that both Γ ± 0 → 0 as κ → ∞, revealing that the short-term decay rates are lower when κ is not too small-a potential advantage to prolong the speed-up of quantum walks.

D. Ring
If G is a ring, we have In Fig. 7, it is reconfirmed that at most four eigenfrequencies exist in G. The unphysical decay rates Γ ± ∞ still diverge as expected. It is also shown that the degeneracy between ω + ∞ and ω − ∞ is barely broken even with fluctuation on w turned on. This is because the diameter of a ring is ∼ N/2, much larger than the diameter of a complete network or a star network which is only 1 or 2. Thus the fluctuation on w only has a local impact which is averaged at long distance. It is also interesting to see that ω ± 0 and ω ± ∞ all approach 2κ/α asymptotically as κ → ∞. There are no other eigenfrequencies.

E. Square lattice
When G is a √ N × √ N square lattice, the behaviors of eigenfrequencies and decay rates of G (Fig. 8) are similar to those of a ring. Comparing Fig. 8(a) with Fig. 7(a), we see that instead ω ± 0 and ω ± ∞ approach 4κ/α asymptotically as κ → ∞. Since a ring is a 1-D system and a square lattice is 2-D, we conclude that for any finite dimensional system, the spectrum of Eq. (8) should be simpler than infinite dimensional networks.

VI. CONCLUSION
To summarize, we present a multiple-scale perturbation method that works on integro-differential equations in the form of Eq. (1), which can be used to unravel the functional importance hidden in the memory kernel and its related complex dynamics.
The multiplescale perturbation method helps to find closed-form approximate solutions that cannot be derived from regular integral transforms, providing a controllable precision that meets practical needs.
In particular, we study its application to a continuoustime quantum walk on some network G enclosed by a general non-Markovian reservoir E. Such a composite system can be regarded as a quantum walk between possible "target" sites (G) and "error" sites (E) if viewed as a quantum error correction algorithm. We propose two physically-important time scales, a primary time scale T and an auxiliary time scale τ , both existing in the strong-coupling regime where G and E are strongly coupled. Compared to the failure of ordinary perturbations supported by Laplace transform, the multiple-scale method shows sufficient accuracy which should be determined by how fast the memory kernel converges locally. The emergence of a new time scale, as the coupling goes stronger, is closely related to the emergence of non-Markovianity. Next, we investigate the eigenfrequencies and their corresponding decay rates of quantum walks on different regular networks. The speed-up of non-Markovian dynamics is confirmed, which however does not exist in the long run, when the two short-term fast frequencies become negligible and only the two long-term slow frequencies persist, which are smaller than the intrinsic mode frequency (the frequency when no reservoir exists). In addition, the behaviors of quantum walks on rings and square lattices are rather simpler than those on complete networks and star networks, because of the limit from dimensionality.