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

Recent experimental advances in the field of quantum technologies have drastically enhanced our capability to control the quantum coherent dynamics of many-body systems in a variety of physical systems, ranging from atomic and molecular optics, to trapped ions and cavity/circuit quantum electrodynamics. This progress has made real the possibility of experimentally realising quantum simulators [1] as well as implementing the first quantum algorithms [2, 3].

Together with the progress in implementing quantum gates and concatenating 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–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 in eight [11, 12], sixteen [13] and even in more than one hundred qubits [14]. Whether these machines already hold the so-called 'quantum supremacy' or not at present 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 performance and to push them 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 that to be able to follow the quantum dynamics of an appropriate time-dependent Hamiltonian:

$$\begin{eqnarray}&&H(t)=[1-f(t)]\,{H}_{\mathrm{in}}+f(t)\,{H}_{\mathrm{fin}},\quad t\in [{t}_{\mathrm{in}},{t}_{\mathrm{fin}}],\end{eqnarray} \tag{ 1 }$$

f(t) being a generic function of time, with $f({t}_{\mathrm{in}})=0$ and $f({t}_{\mathrm{fin}})=1$. The Hamiltonian ${H}_{\mathrm{in}}$ sets the initial condition as its ground state, while the sought solution to the problem is entailed in the ground state of ${H}_{\mathrm{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 $| {\psi }_{0}(t)\rangle$. Reaching the ground state of ${H}_{\mathrm{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 nonadiabatic 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 speed-up as compared to the classical algorithms, see [17] for a prominent example of an adiabatic Grover search and [18] for a review. Moreover, optimally controlled ramps can provide additional speed-ups [19]. An alternative protocol, which is more general than equation (1), would incorporate an extra (possibly nonlinear) term, such as $f(t)[1-f(t)]{H}_{{\rm{E}}}$, with H_E being a 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 the incoherent fluctuations induced by finite temperature, or more in general by the unavoidable coupling of a system to some external environment. In this case, the quantum state of the system will be mixed, described by a density matrix $\rho (t)$, satisfying a dynamical equation of the form:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=-\displaystyle \frac{{\rm{i}}}{{\hslash }}[H(t),\rho ]+{\mathbb{D}}[\rho ].\end{eqnarray} \tag{ 2 }$$

The first term on the right-hand side describes the coherent unitary time evolution, which is ruled by a many-body time-varying Hamiltonian H(t), according to the quantum annealing protocol (1). The second term on the right-hand side of equation (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, and, in general, it drives the system to a fixed point, which is a steady state or a steady-state manifold (if the fixed point is not unique). As such, it induces a decay in the system towards a steady state. Note that the differential form of equation (2) generally requires a Markovian bath. Understanding the effect of dissipation on AQC amounts to quantifying, in some way, how the state of the system deviates from the final ground state, because of the presence of the extra term ${\mathbb{D}}[\rho ]$ in equation (2).

Does the presence of the environment facilitate the reaching of the final ground state or is it detrimental to the AQC? It is clear that there cannot be a unique answer to this question: the deviations from the unitary case may depend strongly on the form of ${\mathbb{D}}[\rho ]$ in relation to the type of evolution imposed by H(t). 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 into 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 is in time, 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 number of such defects [22], while the mechanism yielding the correct scaling was found later by Zurek [23], roughly dividing the dynamics into 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. [24] for a review). When the annealing velocity is progressively increased, 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 upon increasing the annealing time, and the environment will be dominant for long annealing times. In this regime, one expects the defect formation to be almost independent of the annealing protocol. This picture was confirmed and detailed in [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 a 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 the AQC in several papers, showing that in some cases it may be beneficial to reach the target ground state. Work has been done on comparing the AQC and classical approaches using thermal hopping [27], on the use of quantum diffusion, showing better performance than closed-system quantum annealing [28], and on the crucial role of noise-induced thermalisation in the AQC, which can outperform simulated annealing [29]. In addition, the relation of thermally assisted tunnelling to the quantum Monte Carlo has been studied in [30], and the effect of the noise of a thermal environment [31] as well as decoherence [32] have been studied as well.

It should be kept in mind that the adiabatic dynamics of a dissipative many-body system is linked to the understanding of the Landau–Zener problem of a two-level system coupled to an environment, which is an extensively studied problem in many different areas [33–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 nontrivial examples, where the outcomes of the analysis are not hampered by any approximation. The aim of this work is to present some simple, yet nontrivial examples where the adiabatic dynamics can be analysed in full detail. 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 the 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 equation (2) can be written in the Lindblad form

$$\begin{eqnarray}&&{\mathbb{D}}[\rho ]=\displaystyle \sum _{n}{\kappa }_{n}\left({L}_{n}\rho {L}_{n}^{\dagger }-\tfrac{1}{2}\{\rho ,{L}_{n}^{\dagger }{L}_{n}\}\right),\end{eqnarray} \tag{ 3 }$$

where L_n are the suitable local Lindblad operators that describe the environment (to be defined later) and ${\kappa }_{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 on this point in more detail later.

We will consider quadratic fermionic models, whose Lindblad dynamics can be worked out analytically [39–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 nonadiabatic effects and incoherent transitions due to the external bath. In order to understand to what 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 a case, for incoherent spin decay/pumping, the master equation cannot be mapped into local fermionic operators, and therefore we will resort to a numerical study based on a matrix product operator (MPO) representation of the density matrix [42, 43]. As we will discuss in more detail 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 remember that, in general, the thermodynamic properties of low-dimensional systems can be strongly affected by the dimensionality: for example, thermal fluctuations wash out the quantum fluctuations of finite-temperature systems in one dimension, but not anymore in two dimensions. Furthermore, methods that are very powerful in one dimension, such 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 robust evidence of the fact that an optimal working time for the quantum annealing protocol emerges as a result of the competition between the nonadiabatic 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 an 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 two contributions due to nonadiabaticity and dissipation/decoherence.

For longer 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 defect production is monotonically nonincreasing with the annealing speed, in the system-bath scenario this can emerge even for small systems, eventually being related to the spectral structure of the Liouvillian and not necessarily to the many-body characteristics.

The paper is organised as follows. In section 2 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 3, and for an Ising spin chain (section 4). We end with a discussion of our findings and with the concluding remarks in section 5. Technical details on the calculations for the quadratic fermionic model are provided in the appendix.

One of the simplest (and exactly solvable) models exhibiting a quantum phase transition is the one-dimensional Ising chain [44]. This is defined by the Hamiltonian

$$\begin{eqnarray}&&H(t)=-J\displaystyle \sum _{n}{\sigma }_{n}^{x}{\sigma }_{n+1}^{x}-{\rm{\Gamma }}(t)\displaystyle \sum _{n}{\sigma }_{n}^{z},\end{eqnarray} \tag{ 4 }$$

where ${\sigma }_{n}^{\alpha }(\alpha =x,y,z)$ are the spin-1/2 Pauli operators for the nth spin of the chain, while J and ${\rm{\Gamma }}(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 the translational invariance of the model. We will also work in units of ${\hslash }=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 ${\rm{\Gamma }}(t)$ in time according to a linear ramping, as for example:

$$\begin{eqnarray}&&{\rm{\Gamma }}(t)=-t/\tau ,\quad \mathrm{for}\ \ t\in (-\infty ,0],\end{eqnarray} \tag{ 5 }$$

where τ is related to the ramping speed. The choice of (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 ${{\rm{\Gamma }}}_{{\rm{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}_{\mathrm{in}})\propto -{\sum }_{n}{\sigma }_{n}^{z}$, since ${\rm{\Gamma }}({t}_{\mathrm{in}})=+\infty$), to a doubly degenerate ferromagnetic phase, where the spins are all pointing along the coupling direction x (i.e. the ground state of the final Hamiltonian $H({t}_{\mathrm{fin}})=-{\sum }_{n}{\sigma }_{n}^{x}{\sigma }_{n+1}^{x}$, since ${\rm{\Gamma }}({t}_{\mathrm{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 equation (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:

$$\begin{eqnarray}&&{\sigma }_{n}^{-}=\exp \left({\rm{i}}\pi \displaystyle \sum _{m\lt n}{c}_{m}^{\dagger }{c}_{m}\right){c}_{n},\end{eqnarray} \tag{ 6 }$$

where ${\sigma }_{n}^{\pm }=\tfrac{1}{2}({\sigma }_{n}^{x}\pm i{\sigma }_{n}^{y})$, while c_n (${c}_{n}^{\dagger }$) denotes the fermionic annihilation (creation) operator on site n, obeying the anticommutation relations $\{{c}_{m}^{\dagger },{c}_{n}\}={\delta }_{m,n}$, $\{{c}_{m},{c}_{n}\}=0$. The resulting Hamiltonian for an Ising chain of length L is quadratic in such operators and reads:

$$\begin{eqnarray}&&H(t)=\displaystyle \sum _{m,n}\left[{c}_{m}^{\dagger }\,{A}_{m,n}(t)\,{c}_{n}+\tfrac{1}{2}({c}_{m}^{\dagger }\,{B}_{m,n}\,{c}_{n}^{\dagger }+{\rm{H}}.{\rm{c}}.)\right],\end{eqnarray} \tag{ 7 }$$

where $A,B$ are respectively a symmetric and an antisymmetric L × L matrix whose sole nonzero elements are ${A}_{n,n}(t)=-{\rm{\Gamma }}(t),{A}_{n,n+1}={A}_{n+1,n}=-J/2,{B}_{n,n+1}=-{B}_{n+1,n}=-J/2$. To enforce periodic boundary conditions, the following matrix elements are nonzero as well: ${A}_{1,L}={A}_{L,1}={(-1)}^{{N}_{{\rm{F}}}}J/2$ and ${B}_{L,1}=-{B}_{1,L}={(-1)}^{{N}_{{\rm{F}}}}J/2$, where ${(-1)}^{{N}_{{\rm{F}}}}$ denotes the parity of the number of fermions ${N}_{{\rm{F}}}={\sum }_{n}{c}_{n}^{\dagger }{c}_{n}$, which commutes with H. The Hamiltonian (7) can be exactly diagonalised 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 antiperiodic 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 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 parity-dependent boundary conditions. In more 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 acts uniformly (${\kappa }_{n}=\kappa ,\ \forall n$) on each site n of the chain:

$$\begin{eqnarray}&&({\rm{i}})\ \ {L}_{n}^{(1)}={c}_{n}^{\dagger }\qquad {\rm{pumping}}\,{\rm{mechanism}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&({\rm{ii}})\ \ {L}_{n}^{(2)}={c}_{n}\qquad {\rm{decaying}}\,{\rm{mechanism}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&({\rm{iii}})\ \ {L}_{n}^{(3)}={c}_{n}^{\dagger }{c}_{n}\quad \ {\rm{dephasing}}\,{\rm{mechanism}}.\end{eqnarray} \tag{ 10 }$$

Note that while the fermionic dephasing environment (10) can be directly mapped into the dephasing environment of spins, since there is direct local mapping ${c}_{n}^{\dagger }{c}_{n}\to \tfrac{1}{2}({\sigma }_{n}^{z}+1)$, the same is not true for decay and pumping. Indeed, when mapping from a spin operator ${\sigma }_{n}^{-}$ (${\sigma }_{n}^{+}$) to a fermionic operator c_n (${c}_{n}^{\dagger }$), the Jordan–Wigner transformation (6) includes a nonlocal 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 are studying. The choice of these specific Lindblad operators is motivated by the possibility of studying—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 the fermionic model analytically (7) with the 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 only be discussed in section 4, where we will study the Ising model numerically (4) with a spin-decay mechanism provided by ${L}_{n}^{(1a)}={\sigma }_{n}^{-}$, which generalises equation (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 nonlocal terms would be required. However, they do 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 annealing. The excess energy at a given time t expresses the difference between the instantaneous energy during the annealing, $E(t)=\mathrm{Tr}[H(t)\,\rho (t)]$, where $\rho (t)$ is the solution of the master equation (2) at time t, and the ground-state energy ${E}_{0}(t)=\mathrm{Tr}[H(t)| {\psi }_{0}(t)\rangle \langle {\psi }_{0}(t)| ]$ of the instantaneous Hamiltonian system described by H(t):

$$\begin{eqnarray}&&\varepsilon (t)=\tfrac{1}{L}\{\mathrm{Tr}[H(t)\,\rho (t)]-\langle {\psi }_{0}(t)| H(t)| {\psi }_{0}(t)\rangle \}.\end{eqnarray} \tag{ 11 }$$

Using the aforementioned Bogoliubov transformation, the second term of equation (11) can be computed straightforwardly, while the first term E(t) is nontrivial (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 ${ \mathcal N }\equiv \tfrac{1}{2L}{\sum }_{n}\langle 1-{\sigma }_{n}^{x}{\sigma }_{n+1}^{x}\rangle$, which, in the case of ordered chains and at the end of the annealing, is equivalent to the excess energy $\varepsilon (0)$, apart from the trivial constants.

We first analyse a fermionic system described by the Hamiltonian in equation (7), where antiperiodic boundary conditions are imposed. A Fourier transform drastically helps in the diagonalisation of the unitary problem, since the different momentum modes decouple (see appendix 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 equation (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 appendix B. The crucial point resides in the fact that as for the Hamiltonian, the dissipative part of the Lindbladian with ${L}_{n}^{(1)}={c}_{n}^{\dagger }$ does not mix the various modes at different momenta, once a Fourier transform has been employed. As a consequence, the density matrix at time t factorises into different contributions for the various modes:

$$\begin{eqnarray}&&\rho (t)=\mathop{\displaystyle \bigotimes }\limits_{k}{\rho }_{k}(t).\end{eqnarray} \tag{ 12 }$$

The relevant Hilbert space for each positive momentum k has a dimension of 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 linearised Liouville equations for each k mode (B2). For numerical convenience, we restricted the initial point of the annealing procedure (5) to ${t}_{\mathrm{in}}=-5\tau$, and checked that the results were not appreciably affected by this choice [49]. We studied systems up to $L={10}^{3}$ sites and annealing times up to $\tau ={10}^{3},$ and a fourth-order Runge–Kutta integration procedure with a time step ${\rm{d}}{t}={10}^{-2}$ was employed.

Figure 1 shows the behaviour of the excess energy at the end of the annealing protocol, $\varepsilon (0)$, for various values of the dissipation strength κ, as a function of the annealing time τ. In the absence of dissipation ($\kappa =0$), we recover the Kibble–Zurek (KZ) scaling [46, 47]

$$\begin{eqnarray}&&\varepsilon (\tau )\sim 1/{\tau }^{\gamma }\quad {\rm{with}}\quad \gamma =1/2,\end{eqnarray} \tag{ 13 }$$

which can be obtained via knowledge of the critical Ising exponents associated with the phase transition at ${{\rm{\Gamma }}}_{{\rm{c}}}=1$ across which the system is driven. A finite dissipation $\kappa \gt 0$ induces competition between the KZ mechanism of defect generation due to the crossing of a gapless point (which is progressively reduced, with the increasing annealing time τ), and the production of defects generated by incoherent driving itself. Such competition clearly emerges in figure 1 as nonmonotonic behaviour, which generates an optimal working point for the annealing procedure in the presence of dissipation.

Zoom In Zoom Out Reset image size

Let us now have a closer look at the nonmonotonicity, and focus on the optimal (minimal) value ${\varepsilon }_{\mathrm{opt}}$ reached by the excess energy, as well as on the corresponding annealing time ${\tau }_{\mathrm{opt}}$. Figure 2 displays how such quantities depend on κ. Our numerical data nicely agrees with the power-law behaviour over more than two decades of κ values, such that ${\varepsilon }_{\mathrm{opt}}\sim {\kappa }^{1/3}$ and ${\tau }_{\mathrm{opt}}\sim {\kappa }^{-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.

Zoom In Zoom Out Reset image size

3.1. 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]:

$$\begin{eqnarray}&&| \psi ({t}_{\mathrm{fin}})\rangle =\displaystyle \prod _{k\gt 0}({\alpha }_{k}+{\beta }_{k}{\gamma }_{k}^{\dagger }{\gamma }_{-k}^{\dagger })| 0\rangle .\end{eqnarray} \tag{ 14 }$$

Here $| 0\rangle$ indicates the Bogoliubov vacuum corresponding to the final ground state of H(0), ${\alpha }_{k}$ and ${\beta }_{k}$ are complex amplitudes, while the momentum k can take $L/2$ positive values from 0 to π (see appendix A for details).

In the dissipative case, we will not only have those doubly excited states $| {1}_{k},{1}_{-k}\rangle ={\gamma }_{k}^{\dagger }{\gamma }_{-k}^{\dagger }| 0\rangle$, but also singly excited states such as $| {1}_{k}\rangle ={\gamma }_{k}^{\dagger }| 0\rangle$ and $| {1}_{-k}\rangle ={\gamma }_{-k}^{\dagger }| 0\rangle$, 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 $\langle {1}_{k}| {\rho }_{k}| {1}_{k}\rangle$. We find

$$\begin{eqnarray}&&\tfrac{{\rm{d}}}{{\rm{d}}t}\langle {1}_{k}| {\rho }_{k}| {1}_{k}\rangle =\kappa \,\langle 0| {\rho }_{k}| 0\rangle \,f({\rm{\Gamma }},k),\end{eqnarray} \tag{ 15 }$$

with

$$\begin{eqnarray}&&{[f({\rm{\Gamma }},k)]}^{-1}=1+{\left(\displaystyle \frac{{\rm{\Gamma }}-\cos k+\sqrt{1+{{\rm{\Gamma }}}^{2}-2{\rm{\Gamma }}\cos k}}{{\sin }^{2}k}\right)}^{2}\end{eqnarray} \tag{ 16 }$$

for the specific choice of ${L}_{n}={c}_{n}^{\dagger }$. 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 $\langle 0| {\rho }_{k}| 0\rangle$ can be approximated by its initial value 1. Note that since the density of defects ${ \mathcal N }$ is written as ${ \mathcal N }={\sum }_{k}{\gamma }_{k}^{\dagger }{\gamma }_{k}$ in the Bogoliubov basis, excitations of the form $| {1}_{k}\rangle$ only contribute to the positive values of k, while excitations due to coherent dynamics $| {1}_{k}{1}_{-k}\rangle$ contribute to both k and $-k$. Following this, the incoherent part of the density of defects can be estimated according to

$$\begin{eqnarray}&&{{ \mathcal N }}_{\mathrm{inc}}=\displaystyle \frac{\kappa }{L}\displaystyle \sum _{k\gt 0}{\int }_{-\infty }^{0}{\rm{d}}t\ f[{\rm{\Gamma }}(t),k]=\displaystyle \frac{1}{2}\kappa \,\tau ,\end{eqnarray} \tag{ 17 }$$

where the last equality has been obtained after a change of variables from t to ${\rm{\Gamma }}(t)=-t/\tau$, and observing that the summation over $k\gt 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:

$$\begin{eqnarray}&&{ \mathcal N }\sim {{ \mathcal N }}_{\mathrm{KZ}}+{{ \mathcal N }}_{\mathrm{inc}}=\displaystyle \frac{1}{2\pi \sqrt{2}}\,{\tau }^{-1/2}+\displaystyle \frac{1}{2}\kappa \,\tau .\end{eqnarray} \tag{ 18 }$$

From this expression for the total density of defects, the optimal annealing time minimising defect production can be thus estimated by the condition ${\partial }_{\tau }{ \mathcal N }(\tau ){| }_{{\tau }_{\mathrm{opt}}}=0$. A direct calculation gives

$$\begin{eqnarray}&&{\tau }_{\mathrm{opt}}={\left(\displaystyle \frac{1}{2\pi \sqrt{2}}\right)}^{2/3}{\kappa }^{-2/3},\end{eqnarray} \tag{ 19 }$$

with the corresponding density of defects

$$\begin{eqnarray}&&{{ \mathcal N }}_{\mathrm{opt}}={ \mathcal N }({\tau }_{\mathrm{opt}})=\displaystyle \frac{3}{2}{\left(\displaystyle \frac{1}{2\pi \sqrt{2}}\right)}^{2/3}{\kappa }^{1/3}.\end{eqnarray} \tag{ 20 }$$

The predictions given by these equations are in nice agreement with our numerical data shown in figure 2, keeping in mind that $\varepsilon =2{ \mathcal N }$.

To further highlight the role of dissipation during the annealing procedure, we also analysed the excess energy at the end of the annealing, after subtracting the corresponding excess energy in the absence of dissipation:

$$\begin{eqnarray}&&{\rm{\Delta }}(\tau )=\varepsilon (\kappa ,\tau )-\varepsilon (\kappa =0,\tau ).\end{eqnarray} \tag{ 21 }$$

Note that in order to properly define the quantity Δ, we have manifested the κ-dependence of ε in equation (21). After rescaling such a quantity as ${\rm{\Delta }}(\tau )\to {\rm{\Delta }}(\tau )/\kappa$, we observe a fairly good data collapse with τ, as plotted in figure 3. In addition, our data obeys 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 figure 1). We checked that the behaviour of $\varepsilon (\tau )-\varepsilon (+\infty )$ towards saturation decays with a power law as $\sim {\tau }^{-1}$, which is in accordance with [53].

Zoom In Zoom Out Reset image size

The observations made above point toward the substantial independence of the role played by 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.

3.2. Interplay between pumping and decay

Here, we study the interplay between the pumping and decaying mechanism, and the question of whether the steady 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}_{n}^{(2)}={c}_{n}$. The behaviour of the final excess energy $\varepsilon (\tau )$ as a function of the annealing time is shown in figure 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 incoherent pumping (see figure 1). The nonmonotonic behaviour of $\varepsilon (\tau )$ reveals the presence of an optimal working point, where the number of defects is minimal. However, for larger times τ we also recognise the appearance of an overshooting point, where the energy defects become larger than those reached for infinitely slow annealing. Here as well, we checked that the behaviour of $\varepsilon (\tau )-\varepsilon (+\infty )$, after such an overshooting point, decays with a power law as $\sim {\tau }^{-1}$, and again a linear scaling with κ [53].

Zoom In Zoom Out Reset image size

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 ${\rm{\Gamma }}({t}_{\mathrm{in}})=+\infty$, since it would tend to empty the system and thereby increase the energy there. Consequently, in the limit $\tau \to \infty$, where we can assume that the system will always be in the instantaneous steady state, its energy will be $E(t)=2{\rm{\Gamma }}(t)\gt 0$, so it will approach its final value $E({t}_{\mathrm{fin}})=0$ from above. Since for $1\ll \tau \lt \infty$ 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 a finite τ will result in the observed overshooting.

In figure 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 section 3.1. Note that the argument leading to the scaling predictions for a pumping environment holds as well for a decaying environment, only the function in equation (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\tau$ used to replace the initial value of the field by $-\infty$, which is reasonable since we have seen that this environment creates defects long before the quantum critical region is reached. As a consequence, the scaling behaviour of ${\varepsilon }_{\mathrm{opt}}$ and the corresponding τ behaves in accordance with equations (19)–(20).

Zoom In Zoom Out Reset image size

In the lower panel of figure 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 the annealing time ${\tau }_{\max }$ scales linearly with κ, while the change of the maximum excess energy ${\varepsilon }_{\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 figure 6 we show the instantaneous excess energies for different annealing times during the protocol. For very small annealing times we see that the instantaneous steady-state energy is a long way from the actual dynamics and no overshooting takes place. For long annealing times, the excess energy increases hugely at the beginning and then follows (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 overshooting can take place, and the final excess energy $\varepsilon (\tau )$ will be bigger than in the infinite-time limit $\varepsilon (\infty )$.

Zoom In Zoom Out Reset image size

To underline the difference between the two kinds of dissipation (pumping/decay), in figure 7 we plotted the instantaneous excess energy for the same parameters ($\tau ={10}^{3}$, $\kappa ={10}^{-1}$), but with a 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 completely different behaviour, rather following adiabatically the instantaneous steady state of the system.

Zoom In Zoom Out Reset image size

Now, we turn our attention to the interplay between the 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, $\eta ={\kappa }_{\mathrm{pump}}/{\kappa }_{\mathrm{decay}}$. In figure 8 we show the results for values of η ranging from 0 (no pumping), where we observe the biggest overshooting, to 1, where the overshooting completely disappears. Note that the optimal working point does not change much with a varying η since, for the given parameters, the contribution by the pumping mechanism in this regime is far smaller than the one by the decay. The scaling of the maximum value of the overshooting as a function of η is shown in figure 9. An explanation for the diminishing overshooting can be given when looking at the dependence of the instantaneous steady-state energy during the protocol: if η is smaller than 1, it decreases from an initially positive value to 0 linearly, such that an overshooting is possible. For $\eta =1$, the instantaneous steady-state energy is constant equal to 0 such that an overshooting due to adiabatic dynamics is prevented. For $\eta \gt 1$, the instantaneous steady-state energy approaches 0 linearly from below, again preventing an overshooting.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we comment on the issue of thermalisation: a single qubit subjected both to incoherent pumping (${L}^{(1)}={c}^{\dagger }$) 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) factorises into many Hamiltonians (each one describing a mode of pseudo-momentum k) whose Hilbert spaces are, each of them, essentially two-dimensional, this might raise the question of whether 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 ($\gt 98 \%$) for the complete range of physical relevant coupling strengths. However, the corresponding inverse temperature β (${k}_{{\rm{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 β.

3.3. Dephasing

Up to now, all the discussion has been based on a system-bath 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}_{n}^{(3)}={c}_{n}^{\dagger }{c}_{n}$ (which are proportional to the onsite fermionic number operator), as in equation (10). As detailed in appendix C, despite the translational invariance, in such cases the solution to the master equation (2) cannot be trivially written in a tensor structure as that in equation (12). As a matter of fact, the Lindbladian ${\mathbb{D}}[\rho ]$ now transforms into a nonlocal 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 equation (C15). By employing a fourth-order Runge–Kutta integration procedure for these equations, with a time step ${\rm{d}}{t}={10}^{-2}$, we were able to reach annealing times up to $\tau ={10}^{3}$.

The main results of our analysis are summarised in figure 10, where we plot (upper panel) the excess energy $\varepsilon (\tau )$ at the end of the annealing, as a function of the annealing time τ. Comparing this data with that of figure 1, we immediately recognise a qualitatively analogous trend as for the pumping mechanism. In particular, the nonmonotonic 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 ${\rm{\Delta }}(\tau )$ rescaled by κ (bottom panel). Its scaling with τ is completely analogous to that in figure 3, with the data growing linearly with the annealing time, and eventually deviating for sufficiently large values of κ and τ.

Zoom In Zoom Out Reset image size

Finally, we recall that the argument of section 3.1 for determining the scaling of the optimal working point for the annealing protocol as a function of κ also holds in this case. Indeed the corresponding data (with the same power-laws), shown in figure 11, is closely similar to those of figure 2.

Zoom In Zoom Out Reset image size

Summarising 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 three incoherent mechanisms we observe a competition which leads to the onset of an optimal working point for the annealing procedure at a given ${\tau }_{\mathrm{opt}}$ rate. On the other side, for larger values of τ an overshooting point only appears in the presence of a decay mechanism, due to the fact that the instantaneous energy approaches the steady-state value $\varepsilon (\tau )=1$ from below (while the opposite happens for the pumping and for the dephasing). Finally, we analysed how the final excess energy approaches the $\tau \to \infty$ limit, while for pumping and decay we observed the behaviour $| \varepsilon (\tau )-\varepsilon (\infty )| \sim {\tau }^{-1}$; for dephasing we found $| \varepsilon (\tau )-\varepsilon (\infty )| \sim \exp (-\tau )$.

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, equation (4). We first notice that dephasing can be induced by a Lindblad term ${L}_{n}^{(3a)}={\sigma }_{n}^{z}$, which is readily mapped into the local fermionic operator $(2{c}_{n}^{\dagger }{c}_{n}-1)$, through the JWT of equation (6). In such a case, one would thus recover the dephasing mechanism for free fermions (we refer to section 3.3 for details). On the other hand, incoherent pumping/decay would be induced by ${L}_{n}^{(1a)}={\sigma }_{n}^{+}$ and ${L}_{n}^{(2a)}={\sigma }_{n}^{-}$, respectively; in this case, when mapping into fermions, the appearance of the JW string operator forbids analytic treatment, such 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 an MPO [42, 43]. We expect this to be valid whenever the amount of correlations in the system is sufficiently small to satisfy 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 on the right-hand side of equation (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\approx 250$ and adopting a typical Trotter step ${\rm{d}}{t}={10}^{-2}$. We adopted the same time dependence of the field ${\rm{\Gamma }}(t)$ as in equation (5), where for practical convenience we started from ${t}_{\mathrm{in}}=-3\tau$, 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 with that previously discussed in section 3, already for small sizes $L\gtrsim 10$.

Our numerical results for annealing in the presence of incoherent decay, showing the final excess energy as a function of τ and for various dissipation strengths κ, are summarised in figure 12. Despite not being able to see the KZ power-law scaling for limited system sizes (not even in the absence of dissipation), the nonmonotonicity of the various curves for $\kappa \ne 0$ clearly emerges as a result of the open-system dynamics. We ascribe this behaviour to the emerging picture described in section 3, where we discussed much longer systems of free fermions. Indeed, in figure 13 we repeated the same analysis for the scaling of the optimal working time ${\tau }_{\mathrm{opt}}$ and of the corresponding optimal excess energy ${\varepsilon }_{\mathrm{opt}}$ with the dissipation strength, finding similar power-law behaviour. The exponents agree within 20% of the relative difference. We point out that we were not able to fully resolve the overshooting behaviour in this case, since it would have required longer annealing times. However, this is already visible in figure 12, for the curve corresponding to $\kappa =0.1$. Moreover, we checked that the scaling with τ of the final excess-energy difference ${\rm{\Delta }}(\tau )$ was again linear for sufficiently small values of κ and τ, as for the fermionic model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have 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 the adiabatic dynamics of the ground state, showing 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 consistent behaviour for all of them, with the competition between the two processes leading to an optimal working point. This can be modelled by the ansatz of the independent processes, which comes to a scaling behaviour that predicts the observed optimal working point in a fairly accurate way. For larger annealing times, we highlighted the possibility of observing an overshooting point, where defects become larger than those reached for 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 the generic nature of the observed phenomena. Within the framework of free-fermion models, a generalisation to higher dimensions is straightforward.

We thank S Barbarino and D Venturelli for fruitful discussions. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high- performance computing resources and support. SM gratefully acknowledges the support of the DFG via a Heisenberg fellowship. This research was partly supported by the EU projects QUIC and RYSQ, by SFB/TRR21 and by CRP-QSYNC.

Here we provide the technical details concerning the dynamics of the free-fermion model in equation (7), where periodic boundary conditions are imposed, namely,

$$\begin{eqnarray}&&H(t)=-\displaystyle \sum _{n=1}^{L}\{({c}_{n}^{\dagger }{c}_{n+1}+{c}_{n}^{\dagger }{c}_{n+1}^{\dagger }+{\rm{H}}.{\rm{c}}.)+2{\rm{\Gamma }}(t){c}_{n}^{\dagger }{c}_{n}\},\end{eqnarray} \tag{ A1 }$$

with ${c}_{L+1}=-{c}_{1}$ for the positive parity sector, while ${c}_{L+1}={c}_{1}$ for the negative parity sector. Here we have implicitly set the coupling strength to one.

The annealing procedure of equation (5) in this context has already been studied in [47]. The approach consists of employing a Fourier transform of the type

$$\begin{eqnarray}&&{{c}}_{{n}}=\displaystyle \frac{{{\rm{e}}}^{-{i\pi }/4}}{\sqrt{{L}}}\displaystyle \sum _{{k}}{{c}}_{{k}}{{\rm{e}}}^{{\rm{i}}{k}{n}},\end{eqnarray} \tag{ A2 }$$

where the operators ${c}_{k}^{(\dagger )}$ satisfy the canonical anticommutation relations for fermions as well, and the index k takes values (assuming L to be even, without loss of generality) $k=\pm 1\tfrac{\pi }{L},\pm 3\tfrac{\pi }{L},\,...,\,\pm (L-1)\tfrac{\pi }{L}$. The resulting Hamiltonian in Fourier space takes the form

$$\begin{eqnarray}&&H=\displaystyle \sum \,_{k}\{2{c}_{k}^{\dagger }{c}_{k}[-{\rm{\Gamma }}(t)-\cos k]+\sin k({c}_{k}^{\dagger }{c}_{-k}^{\dagger }+{c}_{-k}{c}_{k})\}.\end{eqnarray} \tag{ A3 }$$

Since H conserves the fermionic parity, the global Hilbert space ${ \mathcal H }$ can be written as a direct sum over different k subspaces: ${ \mathcal H }={\oplus }_{k\gt 0}{{ \mathcal H }}_{k}$, and indeed $H={\sum }_{k\gt 0}{H}_{k}$ as in equation (A3). Each subspace at a fixed $k\gt 0$ is built from the two states $\{| 0\rangle ,| {1}_{k},{1}_{-k}\rangle \}$ or $\{| {1}_{k}\rangle ,| {1}_{-k}\rangle \}$, 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 diagonalisation of H_k.

The Hamiltonian H_k at a given fixed time t, as extrapolated from equation (A3), can be readily diagonalised by means of a Bogoliubov transformation [50, 51]

$$\begin{eqnarray}&&{c}_{k}={u}_{k}(t){\gamma }_{k}+{v}_{-k}^{* }(t){\gamma }_{-k}^{\dagger },\end{eqnarray} \tag{ A4 }$$

so that the ground state is annihilated by all the quasiparticle operators ${\gamma }_{k}$. Once the parameter ${\rm{\Gamma }}(t)$ is varied, the system dynamics can thus be found by employing the time-dependent Bogoliubov method [54], which makes the ansatz in which the instantaneous system wave function $| \psi (t)\rangle$ is annihilated by a set of quasiparticle operators ${\tilde{\gamma }}_{k(H)}$, in the Heisenberg representation, defined though the transformation ${c}_{k(H)}={u}_{k}(t){\tilde{\gamma }}_{k(H)}+{v}_{-k}^{* }(t){\tilde{\gamma }}_{-k(H)}^{\dagger }$. This ansatz satisfies the Heisenberg equation

$$\begin{eqnarray}&&\tfrac{{\rm{d}}}{{\rm{d}}t}{c}_{k(H)}={\rm{i}}[H(t),{c}_{k(H)}],\end{eqnarray} \tag{ A5 }$$

with the constraint ${\tilde{\gamma }}_{k(H)}| \psi (t)\rangle =0$, provided the coefficients u_k(t) and v_k(t) obey the time-dependent Bogoliubov–de Gennes equations

$$\begin{eqnarray}{\rm{i}}\tfrac{{\rm{d}}}{{\rm{d}}t}{u}_{k} & = & -2{u}_{k}[{\rm{\Gamma }}(t)+\cos k]+2{v}_{k}\sin k,\\ {\rm{i}}\tfrac{{\rm{d}}}{{\rm{d}}t}{v}_{k} & = & +2{v}_{k}[{\rm{\Gamma }}(t)+\cos k]+2{u}_{k}\sin k.\end{eqnarray} \tag{ A6 }$$

These equations can be integrated starting from the initial condition ${\rm{\Gamma }}(-\infty )=+\infty$, thus mapping them into a Landau–Zener problem [47].

Let us now describe the superimposed action of an environment that induces decay in the system, so the Lindblad operator on each site n is given by ${L}_{n}={c}_{n}$. It is important to stress that when applying the Fourier transform (A2) to the dissipative part of the master equation (3), this does not mix different modes:

$$\begin{eqnarray}&&{\mathbb{D}}[\rho ]\equiv \displaystyle \sum _{k}{{\mathbb{D}}}_{k}[\rho ]=\kappa \left(\displaystyle \sum _{k}{c}_{k}\,\rho \,{c}_{k}^{\dagger }-\tfrac{1}{2}\{\rho ,{c}_{k}^{\dagger }{c}_{k}\}\right).\end{eqnarray} \tag{ B1 }$$

The reason resides in the fact that each term contains two fermionic operators ${c}_{n}^{(\dagger )}$, such as ${\sum }_{n}{c}_{n}\,\rho \,{c}_{n}^{\dagger }\to {\sum }_{n}\tfrac{1}{L}{\sum }_{k,k^{\prime} }{c}_{k}\,\rho \,{c}_{k^{\prime} }^{\dagger }{{\rm{e}}}^{-{\rm{i}}(k-k^{\prime} )n}$ for example, and thus the exponential factor, once summed over n, gives a Kronecker delta ${\delta }_{k,k^{\prime} }$.

As a consequence, the density matrix factorises into $\rho (t)={\otimes }_{k\gt 0}{\rho }_{k}(t)$, and we can decouple the problem into the same k modes as for the nondissipative case. Notice, however, that the dissipation part violates the parity conservation of fermions; therefore here the different subspaces for a given $k\gt 0$ are built up from the four states $\{| 0\rangle ,| {1}_{k}\rangle ,| {1}_{-k}\rangle ,| {1}_{k},{1}_{-k}\rangle \}$, and not simply from two. The Hamiltonian in this basis can be explicitly written as $H={\sum }_{k\gt 0}{H}_{k}$, where

$$\begin{eqnarray*}{{H}}_{{k}}=\left[\begin{array}{cccc}0 & 0 & 0 & 2\sin {k}\\ 0 & -2({\rm{\Gamma }}+\cos {k}) & 0 & 0\\ 0 & 0 & -2({\rm{\Gamma }}+\cos {k}) & 0\\ 2\sin {k} & 0 & 0 & -4({\rm{\Gamma }}+\cos {k})\end{array}\right],\end{eqnarray*}$$

with

$$\begin{eqnarray*}{c}_{k}=\left[\begin{array}{cccc}0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1\\ 0 & 0 & 0 & 0\end{array}\right]\quad {\rm{and}}\quad {c}_{-k}=\left[\begin{array}{cccc}0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0\end{array}\right].\end{eqnarray*}$$

As a matter of fact, solving the full quantum dynamics of $\rho (t)$ translates into solving the $L/2$ Lindblad equations of dimension 4 for ${\rho }_{k}(t)$. In the vectorised form, they can be written as the following linear differential equations with dimension 16 ($k\gt 0$):

$$\begin{eqnarray}\tfrac{{\rm{d}}}{{\rm{d}}t}| {\rho }_{k}\rangle \rangle & = & \{{\rm{i}}({\mathbb{1}}\otimes {H}_{k}-{H}_{k}\otimes {\mathbb{1}})\\ & & -\tfrac{\kappa }{2}[{\mathbb{1}}\otimes ({c}_{k}^{\dagger }{c}_{k}+{c}_{-k}^{\dagger }{c}_{-k})+({c}_{k}^{\dagger }{c}_{k}+{c}_{-k}^{\dagger }{c}_{-k})\otimes {\mathbb{1}}]\\ & & +\kappa ({c}_{k}\otimes {c}_{k}+{c}_{-k}\otimes {c}_{-k})\}| {\rho }_{k}\rangle \rangle .\end{eqnarray} \tag{ B2 }$$

In practice, for every linear operator $W={\sum }_{m,n}{W}_{{mn}}| m\rangle \langle n|$ acting on the four-dimensional Hilbert space ${{ \mathcal H }}_{k}$ which is spanned by the basis $\{| m\rangle \}{}_{m=1,\ldots ,4}$, we associate a vector in the 16-dimensional superoperator space ${{ \mathcal H }}_{k}\otimes {{ \mathcal H }}_{k}$, which is spanned by the basis $\{| m\rangle \otimes | n\rangle \}{}_{m,n=1,\ldots ,4}$, using the convention

$$\begin{eqnarray}&&{W}_{m,n}\to | W\rangle {\rangle }_{\phi },\quad \phi =m+(n-1)L,\end{eqnarray} \tag{ B3 }$$

with $| W\rangle \rangle \equiv {\sum }_{m,n}{W}_{{mn}}| m\rangle | n\rangle$. In this way, we have that $| {W}_{1}{W}_{2}{W}_{3} \rangle\rangle =({W}_{1}\otimes {W}_{3}^{{\rm{T}}})| {W}_{2}\rangle \rangle$, where ${}^{{\rm{T}}}$ denotes the transpose operation. With the vectorisation of the master equation for ${\rho }_{k}(t)$ using this rule, the fact that ${H}_{k}={H}_{k}^{{\rm{T}}}={H}_{k}^{\dagger }$, and that ${c}_{k}^{\dagger }={c}_{k}^{{\rm{T}}}$, we finally arrive at equation (B2). The excess energy (11) is then readily obtained, since $\rho (t)={\otimes }_{k\gt 0}{\rho }_{k}(t)$.

In the case of the dephasing Lindblad operators ${L}_{n}={c}_{n}^{\dagger }{c}_{n}$ on each site, the Fourier transform applied to ${\mathbb{D}}[\rho ]$ turns out to yield a nonlocal object, since each term now contains four fermionic operators, and thus it is not possible to decouple the different k modes. Therefore, this kind of dissipation scheme cannot be directly embedded into the Dziarmaga formalism [47] described above.

In the following, it is more convenient to reduce our study to two-point correlators, since all the relevant quantities for our purposes (such as the excess energy (11), can be expressed in terms of those correlators. This drastically simplifies the analysis into a closed set of differential equations which scale linearly (or at most quadratically, for the nonhomogeneous case) with L [40]. We define

$$\begin{eqnarray}{F}_{m,n} & \equiv & \langle {c}_{m}^{\dagger }{c}_{n}\rangle ,\qquad {G}_{m,n}\equiv \langle {c}_{m}{c}_{n}^{\dagger }\rangle ,\\ {I}_{m,n} & \equiv & \langle {c}_{m}^{\dagger }{c}_{n}^{\dagger }\rangle ,\qquad {K}_{m,n}\equiv \langle {c}_{m}{c}_{n}\rangle .\end{eqnarray} \tag{ C1 }$$

Using anticommutation relations for fermions and the fact that ${({c}_{m}{c}_{n})}^{\dagger }={c}_{n}^{\dagger }{c}_{m}^{\dagger }$, we have ${G}_{m,n}={\delta }_{m,n}-{F}_{n,m}$ and also ${K}_{m,n}^{* }={I}_{n,m}$.

Here we adopt the Heisenberg representation, where the dynamics is described by means of an adjoint Lindblad master equation for a given observable O:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}O={\rm{i}}[H,O]+\tilde{{\mathbb{D}}}[O],\quad {\rm{where}}\end{eqnarray} \tag{ C2 }$$

$$\begin{eqnarray}&&\tilde{{\mathbb{D}}}[O]=\displaystyle \frac{1}{2}\displaystyle \sum _{n=1}^{L}{\kappa }_{n}({L}_{n}^{\dagger }\,[O,{L}_{n}]-[O,{L}_{n}^{\dagger }]\,{L}_{n}).\end{eqnarray} \tag{ C3 }$$

Since in this appendix we deal with the homogeneous cases, we set ${\kappa }_{n}=\kappa$.

We first note that

$$\begin{eqnarray}&&[{c}_{m},H]=2\displaystyle \sum _{n=1}^{L}{A}_{m,n}{c}_{n}+{B}_{m,n}{c}_{n}^{\dagger },\end{eqnarray} \tag{ C4 }$$

where the matrices A and B have been defined in equation (7), with ${J}_{n}={h}_{n}=1$. Moreover, specialising to the dephasing bath ${L}_{n}={c}_{n}^{\dagger }{c}_{n}$ with uniform couplings ${\kappa }_{n}=\kappa$, we have

$$\begin{eqnarray}\tilde{{\mathbb{D}}}[{c}_{m}^{\dagger }{c}_{n}] & = & -\kappa \,{c}_{m}^{\dagger }{c}_{n}(1-{\delta }_{m,n}),\\ \tilde{{\mathbb{D}}}[{c}_{m}{c}_{n}^{\dagger }] & = & -\kappa \,{c}_{m}{c}_{n}^{\dagger }(1-{\delta }_{m,n}),\\ \tilde{{\mathbb{D}}}[{c}_{m}^{\dagger }{c}_{n}^{\dagger }] & = & -\kappa \,{c}_{m}^{\dagger }{c}_{n}^{\dagger }(1+{\delta }_{m,n}),\\ \tilde{{\mathbb{D}}}[{c}_{m}{c}_{n}] & = & -\kappa \,{c}_{m}{c}_{n}(1+{\delta }_{m,n}).\end{eqnarray} \tag{ C5 }$$

The adjoint Lindblad master equation (C2) for the Hamiltonian (A1) and the dephasing bath, referring to the operator ${c}_{m}^{\dagger }{c}_{n}$, reads:

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{c}_{m}^{\dagger }{c}_{n} & = & {\rm{i}}[H,{c}_{m}^{\dagger }]{c}_{n}+{{\rm{i}}{c}}_{m}^{\dagger }[H,{c}_{n}]+\tilde{{\mathbb{D}}}[{c}_{m}^{\dagger }{c}_{n}]\\ & = & 2{\rm{i}}\displaystyle \sum _{j}({A}_{m,j}{c}_{j}^{\dagger }{c}_{n}+{B}_{m,j}{c}_{j}{c}_{n}\\ & & -\,{A}_{n,j}{c}_{m}^{\dagger }{c}_{j}-{B}_{n,j}{c}_{m}^{\dagger }{c}_{j}^{\dagger })-\kappa {c}_{m}^{\dagger }{c}_{n}(1-{\delta }_{m,n})\end{eqnarray} \tag{ C6 }$$

and correspondingly, for the other two-point operators,

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{c}_{m}{c}_{n}^{\dagger } & = & 2{\rm{i}}\displaystyle \sum _{j}(-{A}_{m,j}{c}_{j}{c}_{n}^{\dagger }-{B}_{m,j}{c}_{j}^{\dagger }{c}_{n}^{\dagger }\\ & & +\,{A}_{n,j}{c}_{m}{c}_{j}^{\dagger }+{B}_{n,j}{c}_{m}{c}_{j})-\kappa {c}_{m}{c}_{n}^{\dagger }(1-{\delta }_{m,n}),\end{eqnarray} \tag{ C7 }$$

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{c}_{m}^{\dagger }{c}_{n}^{\dagger } & = & 2{\rm{i}}\displaystyle \sum _{j}({A}_{m,j}{c}_{j}^{\dagger }{c}_{n}^{\dagger }+{B}_{m,j}{c}_{j}{c}_{n}^{\dagger }\\ & & +\,{A}_{n,j}{c}_{m}^{\dagger }{c}_{j}^{\dagger }+{B}_{n,j}{c}_{m}^{\dagger }{c}_{j})-\kappa {c}_{m}^{\dagger }{c}_{n}^{\dagger }(1+{\delta }_{m,n}),\end{eqnarray} \tag{ C8 }$$

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{c}_{m}{c}_{n} & = & 2{\rm{i}}\displaystyle \sum _{j}(-{A}_{m,j}{c}_{j}{c}_{n}-{B}_{m,j}{c}_{j}^{\dagger }{c}_{n}\\ & & -\,{A}_{n,j}{c}_{m}{c}_{j}-{B}_{n,j}{c}_{m}{c}_{j}^{\dagger })-\kappa {c}_{m}{c}_{n}(1+{\delta }_{m,n}).\end{eqnarray} \tag{ C9 }$$

If we now set $l=n-m,l\in [0,L-1]$, for a translational invariant system we can define $\vec{{ \mathcal F }}$ such that ${{ \mathcal F }}_{l=n-m}\equiv {F}_{m,n}$ (and analogously for $\vec{{ \mathcal G }}$, $\vec{{ \mathcal I }}$, $\vec{{ \mathcal K }}$), in such a way that equation (C6) can be rewritten as

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{{ \mathcal F }}_{l} & = & 2{\rm{i}}\displaystyle \sum _{j}({A}_{m,j}\,{{ \mathcal F }}_{l+m-j}+{B}_{m,j}\,{{ \mathcal K }}_{l+m-j}\\ & & -{A}_{m+l,j}\,{{ \mathcal F }}_{j-m}-{B}_{l+m,j}\,{{ \mathcal I }}_{j-m})-\kappa {P}_{l}{{ \mathcal F }}_{l},\end{eqnarray} \tag{ C10 }$$

where ${P}_{l}=1-{\delta }_{l,1}$ ($l=1,\,\ldots L$) are the L components of the vector $\vec{P}$. The first term on the right-hand side can be manipulated as

$$\begin{eqnarray}\displaystyle \sum _{j}{A}_{m,j}{{ \mathcal F }}_{l+m-j} & = & {A}_{m,m}{{ \mathcal F }}_{l}+{A}_{m,m+1}{{ \mathcal F }}_{l-1}+{A}_{m,m-1}{{ \mathcal F }}_{l+1}\\ & = & {A}_{l,l}{{ \mathcal F }}_{l}+{A}_{l,l-1}{{ \mathcal F }}_{l-1}+{A}_{l,l+1}{{ \mathcal F }}_{l+1}\\ & = & \displaystyle \sum _{j}{A}_{l,j}{{ \mathcal F }}_{j}={(A\cdot \vec{{ \mathcal F }})}_{l},\end{eqnarray} \tag{ C11 }$$

where in the second line we used the fact that ${A}^{{\rm{T}}}=A;$ moreover, due to the 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

$$\begin{eqnarray}&&\displaystyle \sum _{j}{B}_{m,j}{{ \mathcal K }}_{l+m-j}=-{(B\cdot \vec{{ \mathcal K }})}_{l},\end{eqnarray} \tag{ C12 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{j}{A}_{l+m,j}{{ \mathcal F }}_{j-m}={(A\cdot \vec{{ \mathcal F }})}_{l},\end{eqnarray} \tag{ C13 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{j}{B}_{l+m,j}{{ \mathcal I }}_{j-m}={(B\cdot \vec{{ \mathcal I }})}_{l}.\end{eqnarray} \tag{ C14 }$$

It is possible to follow the same type of calculations for the other two-point correlators, equations (C7)–(C9), such that the dynamics of all the two-point correlators defined above can be written in a compact way as a set of time-dependent linear equations

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\left[\begin{array}{c}\vec{{ \mathcal F }}\\ \vec{{ \mathcal G }}\\ \vec{{ \mathcal I }}\\ \vec{{ \mathcal K }}\\ \end{array}\right]=[M(t)-\kappa \ \mathrm{diag}(\vec{P},\vec{P},\vec{Q},\vec{Q})]\left[\begin{array}{c}\vec{{ \mathcal F }}\\ \vec{{ \mathcal G }}\\ \vec{{ \mathcal I }}\\ \vec{{ \mathcal K }}\\ \end{array}\right],\end{eqnarray} \tag{ C15 }$$

where the $4L\times 4L$ matrix M is given by

$$\begin{eqnarray}M(t)=2{\rm{i}}\left[\begin{array}{cccc}0 & 0 & -B & -B\\ 0 & 0 & B & B\\ B & -B & 2A(t) & 0\\ B & -B & 0 & -2A(t)\end{array}\right]\end{eqnarray} \tag{ C16 }$$

and the vector $\vec{P}$ has been defined above, while the L components of the vector $\vec{Q}$ are given by ${Q}_{l}=1+{\delta }_{l,1}$ (with $l=1,\,\ldots L$).

Eventually, the instantaneous energy $E(t)\equiv \langle H(t)\rangle$ can be calculated from the quadratic Hamiltonian (A1):

$$\begin{eqnarray}E & = & -\displaystyle \sum _{n=1}^{L}({F}_{n,n+1}-{G}_{n,n+1}-{K}_{n,n+1}+{I}_{n,n+1})-2{\rm{\Gamma }}\,{F}_{n,n}\\ & = & -L({{ \mathcal F }}_{1}-{{ \mathcal G }}_{1}-{{ \mathcal K }}_{1}+{{ \mathcal I }}_{1})-2L\,{\rm{\Gamma }}\,{{ \mathcal F }}_{0},\end{eqnarray} \tag{ C17 }$$

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 generalises equation (A4) to nonhomogeneous quadratic systems [48]:

$$\begin{eqnarray}&&{c}_{i}=\displaystyle \sum _{\mu =1}^{L}({U}_{i,\mu }{\gamma }_{\mu }+{V}_{i,\mu }^{* }{\gamma }_{\mu }^{\dagger }).\end{eqnarray} \tag{ C18 }$$

The transformation satisfies the properties $\langle {\gamma }_{\mu }{\gamma }_{\nu }^{\dagger }{\rangle }_{0}={\delta }_{\mu ,\nu }$ and $\langle {\gamma }_{\mu }{\gamma }_{\nu }{\rangle }_{0}=\langle {\gamma }_{\mu }^{\dagger }{\gamma }_{\nu }^{\dagger }{\rangle }_{0}=\langle {\gamma }_{\mu }^{\dagger }{\gamma }_{\nu }{\rangle }_{0}=0$, where $\langle ...{\rangle }_{0}$ indicates the expectation value over the ground state of $H({t}_{\mathrm{in}})$, that is with ${\rm{\Gamma }}=+\infty$. This yields

$$\begin{eqnarray}{F}_{m,n}({t}_{-\infty }) & = & \langle {c}_{m}^{\dagger }{c}_{n}{\rangle }_{0}={({{VV}}^{\dagger })}_{m,n},\\ {G}_{m,n}({t}_{-\infty }) & = & \langle {c}_{m}{c}_{n}^{\dagger }{\rangle }_{0}={({{UU}}^{\dagger })}_{m,n},\\ {I}_{m,n}({t}_{-\infty }) & = & \langle {c}_{m}^{\dagger }{c}_{n}^{\dagger }{\rangle }_{0}={({{VU}}^{\dagger })}_{m,n},\\ {K}_{m,n}({t}_{-\infty }) & = & \langle {c}_{m}{c}_{n}{\rangle }_{0}={({{UV}}^{\dagger })}_{m,n}.\end{eqnarray} \tag{ C19 }$$

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 these four matrices.