Dissipation in adiabatic quantum computers: Lessons from an exactly solvable model

We introduce and study the adiabatic dynamics of free-fermion models subject to a local Lindblad bath and in the presence of a time-dependent Hamiltonian. The merit of these models is that they can be solved exactly, and will help us to study the interplay between non-adiabatic transitions and dissipation in many-body quantum systems. After the adiabatic evolution, we evaluate the excess energy (average value of the Hamiltonian) as a measure of the deviation from reaching the target final ground state. We compute the excess energy in a variety of different situations, where the nature of the bath and the Hamiltonian is modified. We find a robust evidence of the fact that an optimal working time for the quantum annealing protocol emerges as a result of the competition between the non-adiabatic effects and the dissipative processes. We compare these results with matrix-product-operator simulations of an Ising system and show that the phenomenology we found applies also for this more realistic case.


I. INTRODUCTION
The recent experimental advances in the field of quantum technologies have drastically enhanced our capability to control the quantum coherent dynamics of manybody systems in a variety of physical systems, ranging from atomic and molecular optics, to trapped ions, and cavity/circuit quantum electrodynamics. These progresses have made real the possibility to experimentally realise quantum simulators [1] as well as the implementation of the first quantum algorithms [2,3].
Together with the progresses in implementing quantum gates and concatenate them, i.e., by realising standard circuit computation [4], recently adiabatic quantum computation (AQC) [5] and quantum annealing [6] have received a tremendous boost thanks to the experiments performed with D-Wave machines [7][8][9][10]. The strategy underlying adiabatic quantum computation [5,6] is based on the fact that any quantum algorithm can be formulated in terms of identifying the global minimum (ground state) of a given function (Hamiltonian) over a set of many local minima. On the experimental side, quantum effects were seen to survive on eight [11,12], sixteen [13] and even in more than one hundred qubits [14]. Whether these machines hold already at present the so called "quantum supremacy" or not is still under debate [15,16]. However, it is clearly very important to understand the actual mode of operation of these adiabatic computers, in order to understand the limit of their performances and to push it forward.
A key problem in this framework is to understand the role of dissipation and decoherence on adiabatic quantum computers. This question amounts to understanding the key features that control the adiabatic evolution of a many-body open quantum system. Let us first state the general problem. Suppose to be able to follow the quantum dynamics of an appropriate time-dependent Hamiltonian: (1) f (t) being a generic function of time, with f (t in ) = 0 and f (t fin ) = 1. The Hamiltonian H in sets the initial condition as its ground state, while the sought solution to the problem is entailed into the ground state of H fin . If the control time is much larger than the typical inverse gap between the ground state and the first excited state, the system will adiabatically follow its instantaneous ground state |ψ 0 (t) . Reaching the ground state of H fin by adiabatic evolution is the way AQC works [5,6]. As long as the evolution is unitary, the only source of errors is due to excitations generated by non-adiabatic effects in the dynamical evolution. The complexity of the adiabatic algorithms is reflected in the scaling of the minimum gap with the number of qubits. In general, AQC also requires a special form of f (t) to gain a speedup as compared to classical algorithms, see Ref. [17] for a prominent example of an adiabatic Grover search and Ref. [18] for a review. Moreover, optimal controlled ramps can provide additional speedups [19]. An alternative protocol, which is more general than Eq. (1), would incorporate an extra (possibly non-linear) term, such as f (t) [1−f (t)]H E , with H E being properly chosen Hamiltonian. This may have beneficial effects on the minimum gap of the system, and therefore greatly enhance the AQC performance [20].
It is however clear that, especially for long annealing times, another important source of defects is related to incoherent fluctuations induced by finite temperature, or more in general by the unavoidable coupling of the system to some external environment. In this case the quantum state of the system will be mixed, described by a density matrix ρ(t), satisfying a dynamical equation of the form: The first term on the right-hand side describes the coherent unitary time evolution, which is ruled by a manybody time-varying Hamiltonian H(t), according to the quantum annealing protocol (1). The second term on the right-hand side of Eq. (2) accounts for the coupling to the environment, and its form will depend on the nature of noise and dissipation as well as the form of the coupling between the many-body system and the external bath. The dissipator is a completely positive, trace preserving map, that -in general -drives the system into a fixed point, that is a steady state or a steady-state manifold (if the fixed point is not unique). As such, it is inducing a decay in the system towards the steady state. Note that the differential form of Eq. (2)  . This variety of possible answers is reflected in the wide spectrum of cases already considered in the literature. Within the plethora of possible scenarios, it is however important to establish some general trends that may serve as guidelines in going deeper in this formidable problem.
This type of analysis was first performed [21] in the context of the Kibble-Zurek (KZ) mechanism for defect formation [22,23]. The KZ mechanism is related to the fact that, when crossing a gapless critical point of H(t), no matter how slow the variation of the Hamiltonian in time is, the adiabatic theorem is violated and a finite density of defects will be produced. More than thirty years ago, Kibble put forward a scaling argument aimed at predicting the size and the number of such defects [22], while the mechanism yielding the correct scaling was found later by Zurek [23], roughly dividing the dynamics in either adiabatic or impulsive, according to the distance from the critical point. The KZ mechanism has been tested in a variety of quantum toy models at zero temperature, including ordered and disordered systems, as well as for crossing isolated or extended critical regions (see, e.g., Ref. [24] for a review). When the annealing velocity is progressively increased, a crossover behaviour sets in between the KZ scaling and the generation of excitations due to faster quenches, where the dynamics can be described by the underlying classical model [25].
The proliferation of defects due to Landau-Zener transitions is intimately related to the occurrence of errors in the AQC. While, in the unitary case, the number of defects decreases on increasing annealing time, the environ-ment will be dominant for long annealing times. In this regime, one expects a defect formation which is almost independent on the annealing protocol. This picture was confirmed and detailed in Ref. [21], where the scaling in the crossover between the KZ-dominated regime and the environment-dominated regime was also found. In the presence of spatially correlated noise, additional intermediate regimes emerge due to the comparison of the correlation length of the noise and the correlation length of the system [26].
The role of temperature and external noise was further considered in the context of AQC in several papers, showing that in some cases it may be beneficial in reaching the target ground state. Work has been done on the comparison between AQC and classical approaches using thermal hopping [27], on the use of quantum diffusion, showing better performance than closed-system quantum annealing [28] as well as on the crucial role of noise-induced thermalisation in AQC that can outperform simulated annealing [29]. Further, the relation of thermally-assisted tunneling to quantum Monte Carlo has been studied in Ref. [30]. In addition to this, the effect of noise of a thermal environment [31] and decoherence [32] have been studied as well.
It should be kept in mind that the adiabatic dynamics of dissipative many-body system is linked to the understanding of the Landau-Zener problem of a two-level system coupled to an environment, an extensive studied problem in many different areas [33][34][35][36][37][38].
The problem of describing the adiabatic dynamics of a many-body open quantum system is a formidable problem and approximations are necessary. It is however of fundamental importance to have some non-trivial examples where the outcomes of the analysis are not hampered by any approximation. The aim of this work is to present some simple, yet non-trivial, examples where the adiabatic dynamics can be analysed in full details. Most importantly, the phenomenology that we will extract is related to the dynamics of certain spin systems that are very close to the relevant implementations of AQC. This means that the exactly solvable models that we consider here may be used as a very useful benchmark to test important approximations in more complex cases.
The bath we will deal with is Markovian. This means that the dissipative term in Eq. (2) can be written in the Lindblad form where L n are suitable local Lindblad operators that describe the environment (to be defined later) and κ n are the corresponding couplings, which have to be positive for a Markovian Lindblad master equation. The choice of local Lindblad operators that we are going to study does not lead to a thermal state in the steady state. We will dwell in more details on this point later. We will consider quadratic fermionic models whose Lindblad dynamics can be worked out analytically [39][40][41], for different types of local system-bath coupling. The dynamics of this class of models can be studied exactly, and it will help us to clarify several features of the interplay between non-adiabatic effects and incoherent transitions due to the external bath. In order to understand to which extent the results we find can be applied to more realistic cases, later we will consider a spin-1/2 one-dimensional Ising model. In such case, for the incoherent spin decay/pumping, the master equation cannot be mapped into local fermionic operators, therefore we will resort to a numerical study based on a matrixproduct-operator (MPO) representation of the density matrix [42,43]. As we will discuss in more details in the rest of the paper, the overall phenomenology remains unchanged, thus reinforcing the fact that the exactly solvable models introduced here can be very useful benchmarks. We should remind that, in general, the thermodynamic properties of low-dimensional systems can be strongly affected by the dimensionality: for example, thermal fluctuations wash out quantum fluctuations of finite-temperature systems in one dimension, but not anymore in two dimensions. Furthermore, methods that are very powerful in one dimension, as the Jordan-Wigner transformation, are not applicable in higher dimensions. However, in the special case of free fermions, the system's thermodynamics is not expected to change with the dimensionality.
In all the situations we have addressed, we find a robust evidence of the fact that an optimal working time for the quantum annealing protocol emerges as a result of the competition between the non-adiabatic effects and the dissipative processes. There is an optimal time for which the final state is closest to the true final ground state. The scaling of such optimal time (and that of the corresponding generated defects) can be accurately predicted by assuming that the number of defects produced during the time evolution is a sum of the two contributions due to non-adiabaticity and dissipation/decoherence.
For larger working times it may happen that, depending on the type of system-bath coupling, an overshooting point sets in, where the density of the generated defects is larger than that for an infinitely slow annealing, which would adiabatically drive the system through the instantaneous steady state. While this kind of behaviour cannot appear in the unitary scenario, where the defects production is monotonic non-increasing with the annealing speed, in the system-bath scenario this can emerge even for small systems, being eventually related to the spectral structure of the Liouvillian and not necessarily to many-body characteristics.
The paper is organised as follows. In Sec. II we define the Hamiltonian models, the various dissipation schemes, and the annealing problem under investigation. We then study the departure from the instantaneous ground state, in the presence of dissipative processes which may cause decay or dephasing, both for a translationally invariant free-fermion model III, and for an Ising spin chain (Sec. IV). We end with a discussion of our findings and with the concluding remarks in Sec. V. Technical details on the calculations for the quadratic fermionic model are provided in the Appendix.

II. ADIABATIC DYNAMICS WITH LOCAL DISSIPATION: FROM ISING SYSTEMS TO FERMION CHAINS
One of the simplest (and exactly solvable) models exhibiting a quantum phase transition is the onedimensional Ising chain [44]. This is defined by the Hamiltonian where σ α n (α = x, y, z) are the spin-1/2 Pauli operators for the n-th spin of the chain, while J and Γ(t) are respectively the coupling strength between neighbouring spins and the transverse magnetic field. We assume periodic boundary conditions, in such a way as to preserve translational invariance of the model. We will also work in units of = 1 and set the energy scale by fixing J = 1.
It is possible to realise quantum annealing in the Ising model by tuning the parameter Γ(t) in time according to a linear ramping, as for example: where τ is related to the ramping speed. The choice (5) ensures that, during the annealing procedure, the system will encounter a critical point in which the ground-state energy gap closes and defects will start to appear. Duality arguments [45] show that the phase transition occurs at Γ c = 1. The system is driven from a paramagnetic phase, where all the spins are aligned along the field direction z [i.e., the ground state of the initial Hamiltonian H(t in ) ∝ − n σ z n , since Γ(t in ) = +∞], to a doublydegenerate ferromagnetic phase, where the spins are all pointing along the coupling direction x [i.e., the ground state of the final Hamiltonian H(t fin ) = − n σ x n σ x n+1 , since Γ(t fin ) = 0].
Several papers have already addressed the KZ scaling of defects with τ in the paradigmatic Ising model, both for the clean system [46,47] and for the disordered system [48,49]. Here we are going to add the effect of the coupling to an external environment, modelled through a master equation of the form in Eq. (2).
To retain analytic solvability of the full open-system problem [39,40], we will however start from a mapped version of the Ising chain (4) into a free-fermion model. The latter can be achieved by employing a Jordan-Wigner transformation (JWT), which maps the spin operators in terms of spinless fermions: where σ ± n = 1 2 (σ x n ± iσ y n ), while c n (c † n ) denotes the fermionic annihilation (creation) operator on site n, obeying the anticommutation relations {c † m , c n } = δ m,n , {c m , c n } = 0. The resulting Hamiltonian for an Ising chain of length L is quadratic in such operators and reads: where A, B respectively are a symmetric and an antisymmetric L × L matrix whose sole non-zero elements are To enforce periodic boundary conditions, the following matrix elements are non-zero as well: NF denotes the parity of the number of fermions N F = n c † n c n , which commutes with H. The Hamiltonian (7) can be exactly diagonalized using a Bogolibuov transformation [50,51].
In the following we will completely relax the requirement on fermion-parity dependent boundary conditions, and we simply assume that anti-periodic boundary conditions for the fermions are always enforced. This assumption is perfectly justified for a purely coherent evolution, where the fermion parity is conserved and the initial ground state has an even number of fermions. The reason for enforcing the requirement even when considering the open-system adiabatic dynamics is that, as explained below, our Lindblad operators change the fermion parity, and it would be impossible to solve the problem by using a parity-dependent boundary conditions. More in detail, we will study below the adiabatic dynamics of H(t) under the action of three different types of memoryless local environments. We model them in such a way that the Lindbladian (3) is a sum of terms that act uniformly (κ n = κ, ∀n) on each site n of the chain: ii) L (2) n = c n decaying mechanism, iii) L (3) n = c † n c n dephasing mechanism.
Note that while the fermionic dephasing environment (10) can be directly mapped into a dephasing environment of spins, since there is direct local mapping c † n c n → 1 2 (σ z n + 1), the same is not true for the decay and pumping. Indeed, when mapping from a spin operator σ − n (σ + n ) to a fermionic operator c n (c † n ), the Jordan-Wigner transformation (6) includes a non-local operator (string). We should stress that the naming of these Lindblad operators has been chosen with respect to their action on a system with a local, diagonal Hamiltonian, and not with respect to their actual effect on the system that we study. The choice of these specific Lindblad operators is motivated by the possibility to study, in an essentially exact way, the competition between the unitary dynamics and dissipative effects, focusing on features that do not depend qualitatively on the form of the coupling to the environment.
In the next sections we will first address analytically the fermionic model (7) with dissipation provided by Lindblad terms as in (9)-(10), without making any reference to the mapping with spins. The effect of the nonlocal part of the JWT will be only discussed in Sec. IV, where we will study numerically the Ising model (4) with a spin-decay mechanism provided by L (1a) n = σ − n that generalises Eq. (9) to spins. We restate that this choice of local Lindblad operators does not lead to thermalisation in the steady state, as for this purpose non-local terms would be required. However they cover a special interest since in several experimental implementations, such as circuit-QED, cold-atom settings or trapped ions, this kind of local damping is the relevant one.
In order to quantify the loss of adiabaticity during the annealing protocol, originating both from the closure of the Hamiltonian gap and from the dissipative processes, we are going to study the excess energy ε per site, at the end of the annealing. The excess energy at a given time t expresses the difference between the instantaneous energy during the annealing, is the solution of the master equation (2) at time t, and the ground-state energy E 0 (t) = Tr H(t)|ψ 0 (t) ψ 0 (t)| of the instantaneous Hamiltonian system described by H(t): Using the aforementioned Bogoliubov transformation, the second term of Eq. (11) can be computed straightforwardly, while the first term E(t) is non trivial (see the Appendix). For the Ising spin system we will resort to a fully numerical MPO approach. We point out that a related quantity of interest is the density of defects N ≡ 1 2L n 1 − σ x n σ x n+1 , which, in the case of ordered chains and at the end of the annealing, is equivalent to the excess energy ε(0), apart from trivial constants.

III. FREE-FERMIONIC SYSTEM
We first analyse a fermionic system described by the Hamiltonian in Eq. (7), where anti-periodic boundary conditions are imposed. A Fourier transform drastically helps in the diagonalisation of the unitary problem, since the different momentum modes decouple (see App. A). Note that for the sake of simplicity, here we only consider one-dimensional systems, but our analysis of fermions can be easily extended to larger dimensionalities, since a larger dimension will affect calculations only by changing the Brillouin zone.
Let us concentrate on the case in which each lattice site is coupled to some external bath through a pumping mechanism as in Eq. (8). The master equation (2) during the annealing protocol can be easily integrated via a straightforward generalisation of the time-dependent Bogoliubov method already employed by Dziarmaga [47], as detailed in App. B. The crucial point resides in the fact that, as for the Hamiltonian, the dissipative part of the Lindbladian with L Excess energy ε(τ ) Final excess energy as a function of the annealing time, for the free-fermion model (7) coupled to an environment which induces a pumping mechanism, as in Eq. (8): The various data sets denote different values of the dissipative coupling κ, as listed in the legend. Here we simulated the annealing protocol of Eq. (5) for chains of L = 10 3 sites. Black squares denote data for κ = 0, which obey a power-law behaviour for τ > 1 with the KZ scaling exponent γ = 0.5 (dashed line).
modes at different momenta, once a Fourier transform has been employed. As a consequence, the density matrix at time t factorizes into different contributions for the various modes: The relevant Hilbert space for each positive momentum k has dimension 4, and thus the Liouvillian dynamics can be easily followed inside it. We recall that, for the unitary Schrödinger dynamics, a further decomposition into independent 2 × 2 problems was possible, due to the additional conservation of the fermionic parity (which is now violated by the dissipative decaying terms). The excess energy per site ε during the annealing protocol is thus obtained via a numerical integration of the linearized Liouville equations for each k mode (B2). For numerical convenience, we restricted the initial point of the annealing procedure (5) to t in = −5τ , and checked that the results are not appreciably affected by this choice [49]. We studied systems up to L = 10 3 sites and annealing times up to τ = 10 3 ; a fourth-order Runge-Kutta integration procedure with time step dt = 10 −2 has been employed. Fig. 1 shows the behaviour of the excess energy at the end of the annealing protocol, ε(0), for various values of the dissipation strength κ, as a function of the annealing time τ . In the absence of dissipation (κ = 0), we recover the Kibble-Zurek (KZ) scaling [46,47] ε(τ ) ∼ 1/τ γ with γ = 1/2, which can be obtained by the knowledge of the Ising critical exponents associated to the phase transition at Γ c = 1 across which the system is driven. A finite dissipation κ > 0 induces a competition between the KZ mechanism of defect generation due to the crossing of a gapless point (which is progressively reduced, with increasing annealing time τ ), and the production of defects generated by the incoherent driving itself. Such competition clearly emerges in Fig. 1 as a non-monotonic behaviour, which generates an optimal working point for the annealing procedure in the presence of dissipation.
Let us now have a closer look at the non-monotonicity, and focus on the optimal (minimal) value ε opt reached by the excess energy, and on the corresponding annealing time τ opt . Figure 2 displays how such quantities depend on κ. Our numerical data nicely agree with a power-law behaviour over more than two decades of κ values, such that ε opt ∼ κ 1/3 and τ opt ∼ κ −2/3 . Below we show that this behaviour can be easily predicted by assuming that the KZ production of defects is totally independent of that generated by the dissipation. The above mentioned competition is thus explained in terms of an incoherent summation of the two (independent) contributions.

A. Scaling of the optimal point
We start from the observation that, after the annealing procedure, the final state of the closed system can be easily written as a Bogoliubov state where excitations are provided by pairs of quasiparticles with equal and opposite momenta [47]: Here |0 indicates the Bogoliubov vacuum corresponding to the final ground state of H(0), α k and β k are complex amplitudes, while the momentum k can take L/2 positive values from 0 to π (see App. A for details).
In the dissipative case, we will not only have those doubly excited states |1 k , 1 −k = γ † k γ † −k |0 , but also singly excited states such as |1 k = γ † k |0 and |1 −k = γ † −k |0 , which represent further sources of defects. Indeed, by using the Bogoliubov transformation, we can rotate the Master equation in this frame. This allows us to write down the dynamical equation for 1 k |ρ k |1 k . We find with for the specific choice of L n = c † n . In the adiabatic regime where the KZ scaling argument holds and for small dissipation, the density of defects is much smaller than 1, so that 0|ρ k |0 can be approximated by its initial value 1. Note that, since the density of defects N is written as N = k γ † k γ k in the Bogoliubov basis, excitations of the form |1 k only contribute to the positive values of k, while excitations due to coherent dynamics |1 k 1 −k contribute to both, k and −k. Following this, the incoherent part of density of defects can be estimated according to where the last equality has been obtained after a change of variables from t to Γ(t) = −t/τ , and observing that the summation over k > 0 after the integral over Γ yields a constant factor L/2. Assuming now that the mechanisms of defect generation due to KZ and due to dissipation are unrelated [52], we have: From this expression for the total density of defects, the optimal annealing time minimizing the defects production can be thus estimated by the condition ∂ τ N (τ )| τopt = 0. A direct calculation gives with a corresponding density of defects The predictions given by these equations are in nice agreement with our numerical data shown in Fig. 2, keeping in mind that ε = 2N .
To further highlight the role of the dissipation during the annealing procedure, we have also analysed the excess energy at the end of the annealing, after subtracting the corresponding excess energy in the absence of dissipation: Note that, in order to properly define the quantity ∆, we have manifested in Eq. (21) the κ-dependence of ε. After rescaling such quantity as ∆(τ ) → ∆(τ )/κ, we observe a fairly good data collapse with τ , as plotted in Fig. 3. In addition, our data obey a linear scaling as a function of the annealing time except for deviations induced by bigger values of κ (rather than by longer annealing times τ ) in the regime where the excess energy is nearly saturated to its maximal value (see also Fig. 1). We have checked that the behaviour of ε(τ ) − ε(+∞) towards saturation decays with a power law as ∼ τ −1 , which is in accordance with Ref. [53]. The observations made above point toward a substantial independence of the role played by the dissipation, with respect to the KZ mechanism. The incoherent coupling to the external bath acts uniformly and irrespective of the adiabaticity condition ruled by the ground-state energy gap.

B. Interplay between pumping and decay
Here we study the interplay between pumping and decaying mechanism and the question whether the steady Excess energy ε(τ ) Fig. 1, but for a free-fermion model coupled to an environment which induces a decay mechanism, as in Eq. (9): L (2) n = cn. We observe the same initial trend as for the pumping mechanism, however for longer annealing times we observe an overshooting before the saturation sets in. state of a system subject to both mechanisms is thermal or not. For this, first we focus on the annealing protocol in the presence of a uniform incoherent decay mechanism only, induced by the Lindblad operators L (2) n = c n . The behaviour of the final excess energy ε(τ ) as a function of the annealing time is shown in Fig. 4 for different values of the dissipation strength κ. As one can see from the figure, at relatively small annealing times the trend is qualitatively analogous to that obtained for the incoherent pumping (see Fig. 1). The non-monotonic behaviour of ε(τ ) reveals the presence of an optimal working point, where the number of defects is minimal. However for larger times τ we also recognize the appearance of an overshooting point, where the energy defects become larger than those reached for an infinitely slow annealing. Here as well, we have checked that the behaviour of ε(τ ) − ε(+∞), after such overshooting point, decays with a power law as ∼ τ −1 , and again a linear scaling with κ [53].

FIG. 4. Same plot as in
To better highlight the overshooting behaviour, let us recall that, contrary to the incoherent pumping mechanism, the incoherent decay will drastically affect the completely filled ground state of the initial Hamiltonian at Γ(t in ) = +∞, since it would tend to empty the system and thereby increasing the energy in the system. Consequently in the limit τ → ∞, where we can assume the system always to be in the instantaneous steady state, its energy will be E(t) = 2Γ(t) > 0, so it will approach its final value E(t fin ) = 0 from above. Since for 1 τ < ∞ we know that the dynamics approximately follows this open adiabatic dynamics, it is reasonable to expect that its instantaneous energy will follow a similar trend, in particular it will approach its final value from above as well, and the corrections due to finite τ will result in the observed overshooting. In Fig. 5 (upper panel) we analysed the minimum excess energy that is reached at the optimal working point, and the corresponding annealing time. Their behaviour with κ again follows a power law which is similar to the pumping case, as discussed in Sec. III A. Note that the argument leading to the scaling predictions for a pumping environment holds as well for a decaying environment, only the function in Eq. (16) changes. However, this does not influence the scaling behaviour discussed here, but only the pre-factors. We also stress that, in the decay case, the integral involved in this calculation strongly depends on the value −5τ used to replace the initial value of the field by −∞, which is reasonable since we have seen that this environment creates defects already long before the quantum critical region is reached. As a consequence, the scaling behaviour of ε opt and the corresponding τ behaves in accordance with Eqs. (19)- (20).
In the lower panel of Fig. 5 we have repeated a similar analysis for the maximum excess energy at the overshooting point and the corresponding annealing time, as a function of the dissipation strength. We observe that such annealing time τ max scales linearly with κ, while the change of the maximum excess energy ε max is relatively small, since it varies by less than 10% over almost two orders of magnitude.
For a better understanding of the overshooting, in Fig. 6 we show the instantaneous excess energies for different annealing times during the protocol. For very FIG. 6. Instantaneous excess energy as a function of the external field Γ(t), for various annealing times τ and fixed dissipation strength κ = 0.1 of the decaying environment. One observes that for long annealing times (τ = 1000) the system follows the instantaneous steady state, which has an energy 2Γ(t), and at the end it saturates toward ε(∞) = 1. Intermediate times show the same trend but are not following the open adiabatic dynamics as closely, thus resulting in a higher final excess energy ε(τ ) > ε(∞). For very short annealing times (τ = 1), the influence of the dissipation is much smaller and the final excess energy is smaller than ε(∞).
small annealing times we see that the instantaneous steady-state energy is far away from the actual dynamics and no overshooting takes place. For long annealing times, the excess energy increases hugely in the beginning and then follows the (open) adiabatic dynamics, while the behaviour is similar for intermediate annealing times, but not as drastic. As a consequence, there is an intermediate regime where the annealing time τ is big enough such that an overshooting can take place, and the final excess energy ε(τ ) will be bigger than in the infinite-time limit ε(∞).
To underline the difference between the two kinds of dissipation (pumping/decay), in Fig. 7 we plotted the instantaneous excess energy for the same parameters (τ = 10 3 , κ = 10 −1 ), but different type of dissipation. We observe that, as stated above, in the pumping case the excess energy mostly increases in the last fifth of the protocol, which is close to the quantum critical point. The decaying scenario shows a completely different behaviour, rather following adiabatically the instantaneous steady state of the system. Now, we turn our attention to the interplay between pumping and decay and how the overshooting observed for pure decay is influenced. For this we study the final excess energy as a function of the ratio between pumping and decaying, η = κ pump /κ decay . In Fig. 8 we show the results for values of η ranging from 0 (no pumping), where we observe the biggest overshooting, to 1, where External field Γ(t)  the overshooting is completely disappeared. Note that the optimal working point does not change much with varying η since, for the given parameters, the contribution by the pumping mechanism in this regime are far smaller than the one by the decay. The scaling of the maximum value of the overshooting as a function of η is shown in Fig. 9. An explanation for the diminishing overshooting can be given when looking at the dependence of the instantaneous steady-state energy during the proto- col: if η is smaller than 1, it decreases from an initial positive value to 0 linearly, such that an overshooting is possible. For η = 1, the instantaneous steady-state energy is constant equal to 0 such that an overshooting due to adiabatic dynamics is prevented. For η > 1, the instantaneous steady-state energy approaches 0 linearly from below, again preventing an overshooting. Finally, we comment on the issue of thermalisation: A single qubit subjected to both incoherent pumping (L (1) = c † ) and decay (L (2) = c) processes would relax to a thermal state whose inverse temperature β is related to the ratio of the strengths of the two Lindblad operators. Since in our case the translational invariant quadratic Hamiltonian H(t) factorizes into many Hamiltonians (each one describing a mode of pseudomomentum k) whose Hilbert spaces are, each of them, essentially two-dimensional, this might raise the question if our system shows thermalisation as well. Indeed, the steady state of each of these modes can be approximated by a thermal state with very high fidelity (> 98%) for the complete range of physical relevant coupling strengths. However, the corresponding inverse temperature β (k B = 1) of each mode depends on k, and therefore the complete steady state is not well approximated by a thermal state of a single parameter β.

C. Dephasing
Up to now all the discussion was based on a systembath coupling scheme which induces a decay/pumping mechanism. There is however a complementary effect of decoherence, where the dissipation can generate pure dephasing. This can be easily obtained through diagonal Lindblad terms L Excess energy ε(τ )  the onsite fermionic number operator), as in Eq. (10). As detailed in App. C, despite the translational invariance, in such case the solution to the master equation (2) cannot be trivially written in a tensor structure as that in Eq. (12). As a matter of fact, the Lindbladian D[ρ] now transforms into a non-local object, where the different momentum modes are now coupled together. Therefore it is more suitable to solve a close set of 4L differential linear equations for the relevant two-point correlators [40], see Eq. (C15). By employing a fourth-order Runge-Kutta integration procedure of those equations, with time step dt = 10 −2 , we were able to reach annealing times up to τ = 10 3 .
The main results of our analysis are summarized in Fig. 10, where we plot (upper panel) the excess energy ε(τ ) at the end of the annealing, as a function of the annealing time τ . Comparing these data with those of Fig. 1, we immediately recognize a qualitatively analogous trend as for the pumping mechanism. In particular, the non-monotonic behaviour again reveals a competing effect between the KZ mechanism and the incoherent dephasing. Quantitative differences are barely visible on the scale of the two figures. We observed a slight worsening of the annealing protocol, for the same value of τ , the excess energy being slightly larger than that of the previous case. As we did previously, we also analysed the excess-energy difference ∆(τ ) rescaled by κ (bottom panel). Its scaling with τ is completely analogous to that in Fig. 3, with data growing linearly with the annealing time, and eventually deviating for sufficiently large values of κ and τ . Finally we recall that the argument of Sec. III A for determining the scaling of the optimal working point for the annealing protocol as a function of κ holds also in this case. Indeed the corresponding data (with the same power-laws), shown in Fig. 11, are closely similar to those of Fig. 2).
Summarizing the results of our analysis on the quantum annealing in a translationally invariant free-fermion model interacting with a local environment, the emerging scenario for the different types of dissipation is the following. For all the three incoherent mechanisms we observe a competition which leads to the onset of an optimal working point for the annealing procedure at a given τ opt rate. On the other side, for larger values of τ an overshooting point appears only in the presence of a decay mechanism, due to the fact that the instantaneous energy approaches the steady-state value ε(τ ) = 1 from below (while the opposite happens for the pumping and for the dephasing). Finally, we analysed how the final excess energy approaches the τ → ∞ limit: while for pumping and decay we observed a behaviour |ε(τ ) − ε(∞)| ∼ τ −1 , for dephasing we found |ε(τ ) − ε(∞)| ∼ exp(−τ ).

IV. ISING CHAIN
Let us now go back to the spin-1/2 language and discuss the effects of the coupling to an external bath on the quantum annealing of the Ising chain, Eq. (4). We first notice that dephasing can be induced by a Lindblad term L (3a) n = σ z n , which is readily mapped into the local fermionic operator (2c † n c n − 1), through the JWT of Eq. (6). In such case, one would thus recover the dephasing mechanism for free fermions (we refer to Sec. III C for details). On the other hand, incoherent pumping/decay would be induced by L (1a) n = σ + n and L (2a) n = σ − n , respectively; in that case, when mapping into fermions, the appearance of the JW string operator forbids an analytic treatment as the one discussed previously. Let us thus concentrate on the latter scenario.
We employ a numerical method based on an efficient approximation of the many-body density matrix in terms of a MPO [42,43]. We expect this to be valid whenever the amount of correlations in the system is sufficiently small to satisfy an area-law scaling for the bipartite entanglement in the operator space. The time evolution is performed by means of the time-evolving block decimation (TEBD) algorithm, after a Trotter decomposition of the Liouvillian superoperator in the right-hand side of Eq. (2). In our simulations of the annealing protocol (1) for the Ising model we considered systems up to L = 20 sites, using MPOs with a bond link m ≈ 250 and adopting a typical Trotter step dt = 10 −2 . We adopted the same time dependence of the field Γ(t) as in Eq. (5), where for practical convenience we started from t in = −3τ , and verified that (on the scales of the figures shown below) the results are not affected by this choice. As detailed below, we found an emerging physical scenario which is consistent to that previously discussed in Sec. III, already for small sizes L 10.
Our numerical results for the annealing in the presence of incoherent decay, showing the final excess energy as a function of τ and for various dissipation strengths κ, are summarized in Fig. 12. Despite the KZ power-law scaling cannot be seen for limited system sizes (not even in the absence of dissipation), the non-monotonicity of the various curves for κ = 0 clearly emerges as a result of the open-system dynamics. We ascribe this behaviour to the emerging picture described in Sec. III, where we discussed much longer systems of free fermions. Indeed, in Fig. 13 we repeated the same analysis for the scaling of the optimal working time τ opt and of the corresponding optimal excess energy ε opt with the dissipation strength, finding a similar power-law behaviour. The exponents do agree within 20% of relative difference. We point out that we were not able to fully resolve the overshooting  behaviour in this case, since it would require longer annealing times. However this is already visible in Fig. 12, for the curve corresponding to κ = 0.1. Moreover, we also checked that the scaling with τ of the final excess-energy difference ∆(τ ) is again linear for sufficiently small values of κ and τ , as for the fermionic model.

V. CONCLUSIONS
We presented an extensive study of the adiabatic dynamics of free-fermion models, being driven across their quantum critical point, within an open-system approach using local Lindblad operators. Using the excess energy, we quantified the deviations from adiabatic dynamics of the ground state, showing a competition between the unitary dynamics following a KZ mechanism and incoherent defect generation due to dissipation. While being local, the studied environment covers a wide range of possible sources of dissipation -varying from decay and pumping, separately or simultaneously, to dephasing-and at the same time showing a consistent behaviour for all of them: The competition between the two processes leading to an optimal working point. This can be modelled by the Ansatz of independent processes, which brings to a scaling behaviour that predicts the observed optimal working point in a fairly accurate way. For larger annealing times, we highlighted the possibility to observe an overshooting point, where defects become larger than those reached for an infinitely slow annealing. This effect is intrinsically due to the coupling with an external bath, which drives the system toward the steady state according to the Liouvillian dynamics of the master equation.
Furthermore, we studied the one dimensional Ising chain, which is closely related to the free-fermion models, by means of a matrix-product-operator technique, where we found the same behaviour for small system sizes as well, suggesting a generic nature of the observed phenomena. Within the framework of free-fermion models, a generalization to higher dimensions is straightforward.
The annealing procedure of Eq. (5) in this context has been already studied in Ref. [47]. The approach consists in employing a Fourier transform of the type where the operators c ( †) k satisfy canonical anticommutation relations for fermions as well, and the index k takes values (assuming L to be even, without loss of generality) k = ±1 π L , ±3 π L , . . . , ±(L − 1) π L . The resulting Hamiltonian in Fourier space takes the form (A3) Since H conserves the fermionic parity, the global Hilbert space H can be written as a direct sum over different k subspaces: H = ⊕ k>0 H k , and indeed H = k>0 H k as in Eq. (A3). Each subspace at fixed k > 0 is built from the two states {|0 , |1 k , 1 −k } or {|1 k , |1 −k }, depending on the parity number (even or odd, respectively). The ground state is found in the even parity sector, as one can see from the diagonalization of H k .
The Hamiltonian H k at a given fixed time t, as extrapolated from Eq. (A3), can be readily diagonalized by means of a Bogoliubov transformation [50,51] so that the ground state is annihilated by all the quasiparticle operators γ k . The system dynamics once the parameter Γ(t) is varied can thus be found by employing the time-dependent Bogoliubov method [54], which makes the Ansatz that the instantaneous system wave function |ψ(t) is annihilated by a set of quasiparticle operatorsγ k(H) , in the Heisenberg representation, which are defined though the transformation c k(H) = u k (t)γ k(H) + v * −k (t)γ † −k(H) . This Ansatz satisfies the Heisenberg equation with the constraintγ k(H) |ψ(t) = 0, provided the coefficients u k (t) and v k (t) obey the time-dependent Bogoliubov-de Gennes equations i d dt u k = −2u k [Γ(t) + cos k] + 2v k sin k, i d dt v k = +2v k [Γ(t) + cos k] + 2u k sin k. (A6) These equations can be integrated starting from the initial condition Γ(−∞) = +∞, and thus mapping them into a Landau-Zener problem [47].
where P l = 1 − δ l,1 (l = 1, . . . L) are the L components of the vector P . The first term on the right-hand side can be manipulated as j A m,j F l+m−j = A m,m F l +A m,m+1 F l−1 +A m,m−1 F l+1 = A l,l F l + A l,l−1 F l−1 + A l,l+1 F l+1 where in the second line we used the fact that A T = A; moreover, due to translational invariance of the model, it is possible to shift both indices of A together. Proceeding in an analogous way, for the other terms of (C10) we find j B m,j K l+m−j = − B · K l , (C12) where the 4L × 4L matrix M is given by and the vector P has been defined above, while the L components of the vector Q are given by Q l = 1 + δ l,1 (with l = 1, . . . L).
Eventually, the instantaneous energy E(t) ≡ H(t) can be calculated from the quadratic Hamiltonian (A1): (F n,n+1 −G n,n+1 −K n,n+1 +I n,n+1 )−2Γ F n,n where we used anti-commutation relations to express all terms such that the resulting index of the translational invariant correlators is non-negative. To match boundary conditions and parity considerations, we use L odd and are therefore in the negative parity sector. The initial conditions, for a given Hamiltonian H(t) (that is a certain value of Γ) can be immediately found by a Bogoliubov transformation that generalizes Eq. (A4) to non-homogeneous quadratic systems [48]: The transformation satisfies the properties γ µ γ † ν 0 = δ µ,ν and γ µ γ ν 0 = γ † µ γ † ν 0 = γ † µ γ ν 0 = 0, where . . . 0 indicates the expectation value over the ground state of H(t in ), that is with Γ = +∞. This yields From these equations, exploiting the translational invariance of the system, we can choose the initial conditions of the system (C15) by selecting the first column of each of those four matrices.