Theory of open quantum dynamics with hybrid noise

We develop a theory to describe dynamics of a non-stationary open quantum system interacting with a hybrid environment, which includes high-frequency and low-frequency noise components. One part of the system–bath interaction is treated in a perturbative manner, whereas the other part is considered exactly. This approach allows us to derive a set of master equations where the relaxation rates are expressed as convolutions of the Bloch–Redfield and Marcus formulas. Our theory enables analysis of systems that have extremely small energy gaps in the presence of a realistic environment. As an illustration, we apply the theory to the 16 qubit quantum annealing problem with dangling qubits (Dickson et al 2013 Nat. Commun. 4 1903) and show qualitative agreement with experimental results.


Introduction
The theory of open quantum dynamics [1][2][3][4][5] is an important and active area of physics with applications in nanotechnology, chemical physics, quantum biology, and quantum information. In open quantum theories, the system under consideration is assumed to interact with an environment that has many degrees of freedom. Because details of the environmental Hamiltonian are usually unknown, measurable quantities such as temperature T or noise spectrum S(ω) are used to describe the average statistical behavior of the environment. An open quantum model, therefore, provides a set of differential equations that describe the statistical dynamics of the quantum system, taking the temperature and the spectrum of the bath as input parameters.
Increased research and development of technology in quantum computing and simulation is renewing interest in open quantum modeling. One such promising computation scheme is quantum annealing (QA) [6][7][8] (in particular, adiabatic quantum computation [9]). In QA, the system is evolved slowly so that it stays at or near the ground state throughout the evolution. At the end of the evolution, the system will occupy a low-energy state of the final Hamiltonian, which may represent a solution to an optimization or a sampling problem. Another important direction of research is related to the simulation of large quantum systems with various quantum hardware [10] including quantum dots [11], ion traps [12,13], and lattices of coupled superconducting qubits [14,15].
Open quantum dynamics of a QA processor have been studied theoretically [16][17][18][19]. These models assume weak coupling to an environment, which is typically taken to have ohmic spectrum with large high-frequency (HF) content. This limit is well described by the Bloch-Redfield theory [3-5, 17, 20]. Realistic qubits [21], however, suffer from strong interaction with low-frequency (LF) noise (in particular, noise with 1/f-like spectrum). Incoherent dynamics of a qubit coupled to such an environment are described by the Marcus theory [22][23][24]. A complete (hybrid) open quantum model should account for both LF and HF environments. Such a model for a single qubit has been developed and agreement with experiment has been demonstrated [25][26][27]. A generalization of this theory to multiqubit systems has also been developed and compared with experimental observation [28,29]. Several attempts to combine the Bloch-Redfield and Marcus methods have been undertaken in chemical physics and quantum biology (see, for example, [30][31][32]).
In this paper, we expand the work of [28,29]. We provide a systematic and detailed derivation of a hybrid open quantum model, which agrees with the results of [28,29] for problems with large spectral gaps. Our theory, however, can also be applied to small-gap problems with nonstationary Hamiltonians, for which the model in [28,29] is not applicable. We provide an intuitively appealing and computationally convenient form for the transition rates in terms of a convolution between Redfield and Marcus formulas. As an example, we investigate a dissipative evolution of a 16 qubit system strongly interacting with LF noise and weakly coupled to a HF environment [33]. The problem is characterized by an extremely small gap in the energy spectrum of qubits in the middle of annealing. Solving the problem requires the right combination of Bloch-Redfield and Marcus approaches as well as a proper consideration of the calculation basis, which takes into account the nonstationary effects. The results of the present paper can also be applied to any other open quantum system within the validity range of the approximations.
The paper is organized in the following way. Section 2 describes a single-qubit system to provide the necessary intuition before moving to more complicated multiqubit problems. Section 3 formulates the Hamiltonian and introduces important definitions and notations for the system and the bath. Master equations for the probability distribution of the quantum system are derived in section 4. Section 5 presents the relaxation rates as convolution integrals of the Bloch-Redfield and Marcus envelopes. In the section 6 we show that in equilibrium the master equations obey the detailed balance conditions. We also demonstrate that the convolution expression for the relaxation rates turns into the Marcus or to Bloch-Redfield formulas in the corresponding limits. Dissipative dynamics of a 16 qubit system with an extremely small energy gap is considered in section 7. A brief compilation of the commonly encountered notations is presented in appendix A. In other appendixes we provide a detailed derivation of many important formulas.

Single-qubit system
We begin with a single-qubit system interacting with a heat bath. The total Hamiltonian H of the problem is a sum of three terms: where H S is the Hamiltonian of the single-qubit system, H B is the Hamiltonian of the bath, and the Hamiltonian H int describes an interaction between the qubit and the bath. The qubit has a Hamiltonian

W = D +
We assume that the free bath (with no coupling to the system) has a Gaussian statistics [2,20] determined by a spectrum of fluctuations: is a quantum-mechanical operator of the bath. The specific form of the bath Hamiltonian H B is not important here, although frequently the Gaussian bath is represented as a collection of quantum harmonic oscillators [3] (see appendix B for more information). The system-bath interaction is commonly determined by the Hamiltonian In many realistic situations (see, for example, [27]), the noise may come from different sources, some dominating at low frequencies (such as 1/f noise) and others dominating at high frequencies. As such, we consider S(ω) to be a sum of two terms: where S L (ω) and S H (ω) are functions that are peaked at low and high frequencies, respectively. Each function may tail into the other function's region. Hereafter we refer to the noise with the spectrum (5) as the hybrid noise. The formula of HF spectrum S H is given in section 3.5. For the explicit expression for S L (ω) we refer to equation (B3) shown in the supplementary information section of [8], although there is no need for these formulas here. Notice also that in the present paper we operate with the experimentally-measured parameters of the LF bath, such as the noise intensity W 2 and the reorganization energy ε L , which are defined below. The relaxation dynamics of the qubit become simple in two situations. First, when the qubit is weakly coupled to only a HF bath and the energy splitting of the qubit is larger than the Bloch-Redfield relaxation rate [3][4][5]: The Bloch-Redfield approach is valid when Γ = Ω 0 . Herein, we make the following assumptions for the Boltzmann and Planck constants: The second case is when the qubit is coupled only to a LF bath, and its tunneling amplitude Δ is much smaller than the energy broadening W caused by LF noise: W D  . The linewidth W and the reorganization energy ε L (the shift of the bath energy due to the change of the qubit state) can be found from the formulas The qubit dynamics therefore becomes incoherent and the resulting macroscopic resonant tunneling (MRT) rate is given by [25,26] We note that frequency integrals in equations (7) are taken between -¥ and +¥ unless otherwise specified. The reorganization energy ε L is defined as a principal value integral. In the case of pure 1/f noise, where S L 1 w~w ( ) , the integrals in (7) diverge. To avoid this problem we prefer to use parameters W and ε L , which are directly extracted from experiments such as performed in [27,34]. We also note that in the real systems (see, e.g. [35]) the LF spectrum is described as S L w w a -( ) | | where 1 a~. Even a small deviation of α from one removes divergencies in integrals (7).
The fluctuation-dissipation theorem leads to W 2 =2ε L T, where T is the equilibrium temperature of the bath [25]. Equation (8), commonly known as the Marcus formula [22,23,30,31], is valid when the tunneling amplitude Δ and the rate Γ (8) are much smaller than the MRT linewidth W: Δ = W, Γ = W. This equation has been successful in explaining experimental data from flux qubits [27,34].
In practice, low and HF noises coexist and both have to be considered in the dynamics of the qubit. In [25,27] the formula (8) has been generalized to include effects of HF noise on the MRT rate. For small tunneling amplitudes Δ the modified rate Γ is described by the following integral: We notice that the integrand of equation (9) The HF factor is characterized by the more complicated integral: We notice that both functions, G L (ω) and G H (ω), satisfy the normalization condition, where μ=L, H. The rate (9) can be represented as a convolution of the Gaussian envelope G L (ω) and the function G H (ω), Here we need not to assume that the qubit-bath coupling is small. Equation (14) is valid at frequencies 0 1 , t is a correlation time of the HF fluctuations described by the function S H (ω). Later we introduce a spectral density S H of the ohmic noise characterized by the correlation time τ H ∼1/T.
The Bloch-Redfield limit is described by equation (11) with the frequency ω, which is much larger than the coupling to the environment given by the spectrum S H w ( ): 1.
With this small parameter, we can expand the dissipative factor in equation (11). Now the function G H (ω) turns into the form Equations (14) and (15) can be approximately combined into one Lorentzian formula that has a frequencydependent numerator, The Markovian and Bloch-Redfield expressions follow from this formula in the corresponding limits. Notice also that the function (16) is normalized according to equation (12). Thus, the single-qubit relaxation rate (9) can be conveniently represented as a convolution of the Gaussian and Lorentzian line shapes, The convolution integral in (17) has a simple interpretation. One can think of the LF noise as a random shift in energy bias: h→h+h noise , where h noise has Gaussian distribution with variance of W 2 . The HF relaxation rate, given by the Lorentzian lineshape, will therefore be shifted by h noise . Ensemble averaging over LF fluctuations will lead to a convolution integral similar to (13) and (17). The reorganization energy ε L is a result of the action of the qubit on the environment. Equation (17) is valid in the limit of small tunneling amplitude Δ = W. In this case the rate Γ is also much less than the MRT line width W: It should be mentioned that the intuitive description outlined above holds beyond the validity of equation (17). In the next sections, we will generalize this approach to multiqubit systems without resorting to the small tunneling amplitude approximation.

The Hamiltonian
We are interested in dissipative evolution of a quantum annealer [8,14,15,33,37]   The energy functions s ( ) and s ( ) determine the annealing schedule with s=t/t f being the dimensionless annealing parameter ( s 0 1   ), and t is and t f being the running time and the total annealing time, respectively. Details of the annealing schedule are unimportant for the current discussion as long as the timedependent Hamiltonian changes slowly, which is exactly the case for QA algorithms.
We assume an interaction with a bath of the form: with operators Q α characterized by Gaussian statistics with zero average values, Q 0. á ñ = a We also suppose that different qubits, labeled as α and β, are coupled to statistically independent environments, such that ¹ This has been experimentally confirmed for flux qubits [38].

Schrödinger picture
The total system-bath Hamiltonian H written in the Schrödinger representation has the form given by equation (1). Here, the time evolution of the system-bath can be described by the density matrix The evolution of the density matrix can now be defined in terms of the unitary operator Hereafter, for simplicity of notation, we remove time dependences from unitary matrices. A consecutive application of the operators U B and U produces the system-bath density matrix ρ SB (t) at time t, This time-dependent matrix presents the solution of the von Neumann equation (21). The average value of an arbitrary Schrödinger operator , which describes a physical variable of the qubits or of the bath, is determined by the density matrix ρ SB (t) (27) taken at time t, Here the total trace Tr includes the trace Tr B over free-bath variables and also the trace Tr S , over a full set of qubit states k . ñ {| }

Heisenberg picture
In the theory of open quantum systems interacting with a heat bath [39] the state of the system can be described by the reduced density matrix, which is obtained by averaging the system-bath matrix ρ SB over the bath fluctuations. Some information about quantum fluctuations is lost after the averaging. This limits the method to calculations of only the averages and same-time correlation functions. Other properties such as different-time correlations remain beyond the reach of this approach. In the Heisenberg picture, the equations are written in terms of the operators without taking averages. This allows calculations of correlation functions to any order as long as the equations can be solved.
In the Heisenberg representation, the average value of an arbitrary operator  in (28) can be written as

The bath
We assume that the bath coupled to α-qubit is described by Gaussian statistics [2,20]. These statistics are characterized by a correlation function K t t , ¢ a ( ) Here Q α (t) is a free-evolving bath operator (25). The brackets á¼ñ denote the average of  over the free-bath fluctuations, unless otherwise specified. For stationary processes, K t t , ¢ a ( )depends on the time difference, hence allowing the spectral density to be defined as In addition to the correlator (34), we introduce dissipative functions f α (t) and g α (t) defined as These definitions can be compared to similar functions introduced by equation (3.36) in [3]. Notice that The total reorganization energy of the bath is defined as The response of the bath to an external field is described by the retarded Green function The causality is provided by the Heaviside step function t t . q -¢ ( ) The response function j α is related to the susceptibility of the bath defined through According to the fluctuation-dissipation theorem, in equilibrium S α (ω) is proportional to the imaginary part c w  a ( ) of the bath susceptibility, where T is the temperature of the equilibrium bath.

Hybrid noise
Hereafter we assume that the dissipative environments coupled to different qubits, although uncorrelated, have the same spectral density of bath fluctuations: S α (ω)=S(ω). The same is true of the functions f f g , , and χ α (ω)=χ(ω). In the case of hybrid noise, S w ( ) is given by equation (5). The dissipative functions f and g can be split into low and HF components, . 4 2 For the LF part of the function f, one can expand e iwt in equation (37), assuming ωτ=1. Keeping up the second order in ωτ, we obtain with ε L and W defined in equation (7).
To treat the HF parts, we assume Ohmic noise Here η is a small dimensionless coupling constant and ω c is a large cutting frequency of the HF noise. This assumption is justified experimentally [27] and also theoretically [3]. The dissipative functions f H and g H are calculated in appendix C. The total reorganization energy, , e e º a is defined by (38), so that .
L H e e e = + Here ε L is defined in (7), and H e is the HF component of the reorganization energy (38). For the Ohmic spectrum (44) of the bath, we have .
We also introduce the following notations, which we will use later:

System evolution in the interaction representation
Properties of the system of qubits are determined by the reduced density matrix where the system-bath density matrix ρ SB is given by equation (27). A system-bath average of an arbitrary operator S  of the system is written as In the basis introduced in section 3.6, the matrix ρ S has the form n m , with the matrix elements defined as Our goal is to derive a set of master equations for the probability distribution of the qubits, P n , over the states nñ The time evolution of the matrix (50) is determined by the unitary operator U defined by (26) where the Hamiltonian H¢ is given by equation (45). The objective is to go beyond the perturbation theory in the systembath coupling. This can be done by treating Q n exactly, but Q mn perturbatively. The interaction representation is best suited for this goal. A transition to the interaction picture, although straightforward for time-independent bases, becomes more involved if the basis changes in time. Let us introduce a unitary operator where nñ | is a time-dependent basis of the system, and is written in terms of average energies E t n t H t n t n S = á ñ ( ) ( )| ( )| ( ) . We also introduce the S-matrix: with  being the time-ordering operator for t 1 . Notice that the time-dependent matrix element t n s a ( ) is taken out of the integral over t 1 . This becomes necessary when we want to express correlation functions in terms of dissipative functions.
The interaction Hamiltonian is given by the expression

This Hamiltonian defines the unitary evolution operator
U t e . 5 9 The Hamiltonian H I does not contain the nonperturbative diagonal terms Q n . To calculate the time-derivative In the interaction picture, the system-bath Hamiltonian H I (58) takes the form The modified bath operator Q mñ has diagonal terms Here we assume that the spectrum E n is nondegenerate and that H S (18) is characterized by real parameters. In this case we have n n 0.
We show that the only terms that survive during the annealing run are characterized by the relation , , In addition, we demonstrate that, during annealing, correlations between diagonal operators Q kk (63) and off-diagonal bath variables Q mñ (64) rapidly disappear in time, such that Q t Q t , 0 .
The same is true for the average values of the operators (64): Q t 0. mn á ñ( )

Time evolution
The evolution of the matrix ρ S is determined by the unitary matrix (26), which can be written as: U=U 0 U I , where U 0 and U I are given by equations (55) Here, we have used (22), (50), (53) and have introduced the interaction picture operator U m n U , 7 3 which will play an important role in our theory. The bath average ... á ñ is defined in (35). Equation (72) takes a simpler form for the diagonal elements: In the following, we consider time evolution of the operators Λ nn instead of working with the elements (72) and (74) of the system density matrix. Working with operators instead of averages allows derivation of more accurate master equations.
Taking the derivative of (73) and using (62), we obtain Here we use the fact that mn L and Q kl I , taken at the same moment of time t, commute: ] for any set of indexes m, n, k, l. The evolution of the diagonal elements nn L is of prime interest since these elements determine the probabilities (74): This equation is exact and difficult to solve without approximations.
To simplify equation (78), in appendix E we use perturbation expansion assuming that the combined operator Q mñ is small. It does not mean that all system-bath interactions are assumed to be small. A part of these couplings described by the operator Q n (t) is treated precisely (see equations (47), (63), (64) for definitions). In appendix E we take into account the above-mentioned perturbative approximation and express the right-hand side of equation (78) in terms of products of bath and system correlation functions, such as K t t t t , In appendix D we find that the bath correlator K t t , mn ¢ ( )(71) has a short correlation time τ mn given by equation (D36). With the assumption that the time evolution of the operator Λ mn (t) is slow on the time scale τ mn , in appendix E we derive the following master equation for the probability distribution P n of the system (74): mn mn mn are defined by equation (49) and where ε is the total reorganization energy described in section 3.5, T T .
mn nm * =¯All matrix elements of the system operators in equation (80) are taken at the running moment of time t. The master equation (79) with the rates (80) can be considered as a generalization of the NIBA (non-interacting blip approximation) approach [3,29,40] to the case of the multi-level quantum system and also to the case of not-so-small tunneling amplitudes and energy gaps. The rate Γ nm can be written in a form similar to the single-qubit expression (9) and also to the multiqubit rate 1 0 G  given by equation (5) from [28] and by equation (68) from [29], Here we use equations (43) (5) and (6) in [28] and equations (43), (52), (54), and (68) in [29] with two important differences. First, T mn in (81) includes the term m n á ñ |˙which can take care of non-adiabatic transitions such as those in the Landau-Zener effect. Second, the rate (82) does not contain any polaron shifts in the definition of ω mn . While polaron shift has little effect on the results obtained in [28,29] due to the relatively large energy gap and small number of qubits, it can have significant effect on small gap multiqubit problems, such as the one studied in section 7, rendering a large disagreement with the experiment.

Applicability conditions
The master equation (79) have been derived in appendix E with the proviso that 1, 84 where Γ nm is the relaxation rate (80 It should be emphasized that we work in an arbitrary basis formed by the time-dependent states n t ñ | ( ) . The specific basis is chosen in view of the requirement that the perturbation parameter Γ nm τ mn is much less than one. In the systems with ultra-small energy gaps the characteristic time mn t (85) calculated in the instantaneous energy basis goes to infinity since both energy parameters, the gap E m −E n and the linewidth W mn , tend to zero. In this basis the parameter Γ nm τ mn is extremely large, and our perturbation treatment is not valid. The situation can be improved by going to another, so-called pointer basis [41], as we do in section 7.

Relaxation rate as a convolution of Bloch-Redfield and Marcus envelopes
In this section we show that, in addition to the expression (82) of the rate nm G as the integral over time, the same rate can be conveniently represented as a convolution integral over frequencies of the Gaussian envelope multiplied by the Lorentzian function. As in the case of a single qubit described in section 2, the Gaussian curve is produced by LF bath noise. The Lorentzian factor is due to effects of the HF environment. Here, our aim is a generalization of the MRT single-qubit formula (17) to the multiqubit case where both, LF noise and the singlequbit tunneling, can be large.

Convolution form of the rate nm G
Using integration by parts and f a f a f mn mn

Let us introduce Fourier transformations
Substituting in (86) and taking the integral over τ, we obtain:

Special cases
In this section we verify detailed balance conditions for the equilibrium distribution of qubits and also consider the Bloch-Redfield and Marcus limits of the rates (88). In addition, we apply the results of the previous section to a single qubit interacting with a hybrid environment.

Equilibrium condition
We conclude from equations (88) where ω mn is defined by equations (46) and (66). It follows from equation (79) that the equilibrium probabilities P n eq and P m eq to observe the qubits in the states nñ | and mñ | , respectively, obey the equation: The solution of this equation follows the detailed balance condition: with the local energy levels E m and E n (46) and the bath temperature T. Although there seems to be a freedom of choice in selecting the basis nñ | , in reality equilibrium is established in the pointer basis [41], the basis which is most stable under environmental perturbations. We get back to this point in section 7.

Bloch-Redfield and Marcus limits
For the qubits weakly interacting with the HF noise in the absence of LF noise, the parameters of the LF bath go to zero: W=0, ε L =0. The Gaussian envelope G mn L w ( )(91) is approximated by the function 2πδ(ω), and the rate Γ mn (88) as it follows from equations (85) and (101). The condition (103), which is written in the energy basis, fails near the degeneracy point, h=0, at a small tunneling amplitude Δ. In order to keep our theory alive we have to move to the computation basis formed by eigenstates of σ z -matrix.

Dissipative evolution of a 16 qubit system
In this section we analyze dynamics of the 16 qubit system with the set of couplings and biases determined by the Dickson instance (see figure 1). This instance was proposed in [42] and investigated in details in [33]. The energy spectrum of the problem features an extremely small gap between the ground and first excited states. The existence of such a gap presents a computational bottleneck for QA. An experimental technique to overcome this difficulty by individual tuning qubit's transverse fields has been demonstrated in [43]. Nevertheless, a theoretical analysis of dissipative dynamics in this system presents a real challenge. The probability distribution of the qubits is governed by the master equation (79) with the relaxation matrix given by equation (88). The qubits are described the Hamiltonian H S (18). In the problem Hamiltonian H P (19) we have ferromagnetic couplings between qubits, J ij =−1, for every pair of coupled qubits. Two internal qubits have zero biases, h h 0   figure 2. This annealing schedule, which is taken from [33], has minor variations between the qubits. These variations, although not so important, allow us to reproduce the same extremely small energy gap, E 2 −E 1 =0.011mK, between the ground and the first excited states as in experiments performed in [33]. The gap is located at s 0.6396. * = In figure 3 we show the four lowest energy levels of the system near the anticrossing. The most interesting annealing dynamics happen in the interval s s s ,  More details can be found in [33] and in the supplementary information for that paper. It follows from figure 2(c) of [33] that, before the anticrossing at s s* < , the instantaneous eigenstates of the 16 qubit system coincide with the diabatic states: 1 After the anticrossing point at s s* > , we have the reverse situation, with 1 GM ñ = ñ | | and 2 . ñ = Sñ | | Although the experimental results provided in [33] were in accordance with the physical intuition given in the paper, no theoretical analysis was provided. This was due to the lack of an open quantum theory that takes into account both LF and HF noises. Here, we apply our approach to provide a theoretical explanation of the experimental results of [33].
The presence of a very small gap and the time-dependence of the system Hamiltonian, which becomes nonadiabatic near the minimum gap, make the problem instance in figure 1 difficult to analyze within one theoretical framework in all regions during the annealing. As such, some tricks are necessary to choose the proper basis as we discuss next.

Rotation of the basis
The dissipative dynamics of the qubits coupled to a heat bath is described by the master equation (79). These equations are derived with the proviso that the rate Γ nm of the relaxation (88) between the states nñ | and mñ | is much less than the inverse time scale mn   Both of these break the applicability condition (84) in the instantaneous energy basis. Moreover, in this basis the time dependence of the Hamiltonian can create nonzero off-diagonal elements of the density matrix near the minimum gap due to nonadiabatic transitions. These terms do not decay quickly as required by our theory. As we shall see, all these issues can be resolved by rotating the basis. This is equivalent to the introduction of the pointer basis as described in [28,29,41]. We rotate the two anticrossing states as:

| | ( )
Notice that (106) assures minimum quantum transitions between the two states 1¢ñ | and 2¢ñ | near the anticrossing and therefore minimum generation of off-diagonal elements of the density matrix. As we show in section 7.2, it also resolves all issues with the applicability condition (84) discussed above. It should be mentioned that here the rotation angle Θ is different from the rotation angle * q of the pointer basis in [28,29], where * q maximizes the parameter W 2 1 ¢ ¢ calculated for the states under consideration (see, for example, equations (30), (32) in [29]). In our case the angle Θ minimizes the Landau-Zener transitions between two anticrossing states. Recall that in [28,29] the Landau-Zener transitions are just neglected.
A solution of equation (107) is shown in figure 4(c). In the beginning of annealing Θ ; 0, so that the rotated basis coincides with the instantaneous basis. Near the anticrossing point, at s s* = , the angle Θ rapidly switches to π/2. We notice that, with the condition (106), the renormalized matrix element T 2   it small, within the applicability range of our model. For this example we rotate only the two lowest-energy states that have an anticrossing. However, the rotation can be applied to anticrossing excited states as well, if necessary.

Rates and probabilities for the 16 qubit system
In this section we calculate the rate 1 2 G ¢ ¢ (88) of the system relaxation between the states 1¢ñ | and 2¢ñ | . These states are defined as superpositions     Figure 7 shows that, at the smaller coupling to the HF bath, where 0.1 h = , and at the quite strong interaction of qubits with the LF noise, with W=10mK, the hybrid rate Γ 12 almost coincides with the Marcus rate .

M 12
G The validity of these results is verified by the small parameter 0.008 12 21 t G < shown in figure 7(b). In figure 8, in parallel with the energy spectrum, we plot a time dependence of the probabilities P 1 and P 2 to find the 16 qubit system in the instantaneous eigenstates 1ñ | and 2ñ | of the Hamiltonian H S (18). To calculate these probabilities, we obtain the numerical solution of the master equation (79) for the probabilities P 1¢ and P 2¢ to observe the system in the states 1¢ñ | and 2¢ñ | (105). After that, we rotate the basis back, to the instantaneous energy eigenstates 1ñ | and 2 . ñ | The contribution of the off-diagonal elements of the density matrix to the probabilities P 1 and P 2 rapidly disappears as it follows from equations (96) and (D19). We also show the The evolution of the probabilities P 1 , P 2 to find the system in the instantaneous energy eigenstates during the annealing process at the total anneal time t f =2 ms and at W 0.1, 10 h = = mK, T=10 mK. We also show the time evolution of the probability P GM for the system to be in the GMñ | .
evolution of the probability P GM (blue curve) to find the system in the state GMñ | defined in (104). Here, we have a qualitative agreement with the experimental results shown in figure 2(d) of [33]).

Thermal enhancement of the success probability
The goal of annealing is to reach the ground state at the end of the evolution. It follows from figure 2(c) of [33] that at the end of annealing, at t=t f , the ground state of the 16 qubit system coincides with the state GMñ | shown in equation (104). We therefore define the success probability as the probability P GM to observe the system in GMñ | at t=t f . Figure 3 of [33] demonstrates the temperature dependence of P GM . It is clear from this figure that, at sufficiently fast annealing (t f 100 ms), P GM grows with increasing temperature from 20 to 40mK and decreases after. Our goal is to reproduce this non-monotonic behavior with our open quantum model. It should be noted that we do not aim to precisely fit the experimental data to the results of our model.
We solve numerically the master equation (79) with the relaxation rates given by the convolution formula (88) written in the rotated basis (105). We perform simulation for the total anneal time t f between 0.04 and 4 ms. An example of P GM as a function of time is shown in figure 8(b) in the section 7.2. Figure 9 plots the success probability, P GM (t f ), as a function of temperature T for different speeds of annealing characterized by the anneal time t f . In the theoretical calculations we assume that 0.1 h = and W=20mK. Figure 9 reproduces the results shown in figure 3 of [33], including the enhancement P GM at low temperatures and its reduction at T 40  mK.
As mentioned in [33], this decrease may be related to the excitement of the high-energy levels separated from the two lowest states by a gap of order 40mK (see the spectrum in figure 3).

Conclusions
In this paper, we have derived a set of master equations describing a dissipative evolution of an open quantum system interacting with a complex environment. The environment has low-frequency and high-frequency components, as in the case of realistic qubits affected by the hybrid bath, which includes 1/f and Ohmic noise. A part of the system-bath interaction is treated in a nonperturbative way. This mathematical treatment allows us to combine the Bloch-Redfield and Marcus approaches to the theory of open quantum systems and obtain the relaxation rates, which are well-suited for the description of dissipative dynamics of many-qubit quantum objects, such as quantum annealers, in the whole annealing range. The relaxation rates are expressed in the convenient convolution form clearly showing the interplay between the low-and high-frequency noise. The main results of the paper are given by the master equation (79) with the relaxation rates (88). As an illustration, we apply the theory to the 16 qubit quantum annealer investigated in [33]. The instance studied there features an extremely small gap between the ground and first excited states. With the proper rotation of the basis, we have solved the master equations and theoretically confirmed the main experimental findings of [33]. The results of the paper may be useful for understanding a dissipative evolution of various systems, from chromophores in quantum biology [30][31][32] to qubits in real-world quantum processors [14,15,44].

Appendix A. Notations
In this appendix we assemble notations used throughout the paper, so that the main part of the paper becomes easier to follow. In a chosen basis nñ {| }the matrix elements of the system Hamiltonian H S (18) and those of the Pauli matrix z s a of the α-qubit are denoted as is proportional to the imaginary part of the bath susceptibility, c w ( ), which defines dissipative properties of the bath,

Appendix C. HF dissipative functions
For the ohmic HF bath characterized by the spectrum (44) provided that the cutting frequency ω c is much higher than the temperature, k T.
, all the three functions vanish. In equation (44) we introduce η as a small dimensionless coupling constant [3], and ω c as a large cutting frequency of the HF noise. For the Ohmic bath, the spectral density H c w  ( ) is defined as We assume independent noise sources coupled to every qubit, each described by the above-mentioned formulas with identical parameters.
This correlator can be represented as a sum of four components: .
The bath variable Q mn is defined in (47), and n  is the S-matrix of the bath given by equation (57).
where the modified S-matrix of the bath, t, , ) is defined by equation (61). We notice that The functional t t , ; , mn m n  t t ¢ ¢ ¢ ¢ ( )obeys two differential equations: The solution of these two equations is given by the expression Here the matrix elements m s a and n s a are taken at time t, whereas the elements m s a ¢ and n s a ¢ depend on the time t¢. The total reorganization energy ε α is defined by equation (38).
With equation ( We expect now that the time interval t t -¢ is short enough, so that t t t t . ¢ -¢  | | However, we have t∼t f and t t , f ¢~where t f is the total annealing time. Notice also that the functions f t a ( ) and f t¢ a ( ) are growing with time. It follows from the formulas in appendix C that the HF parts of the functions f α (t) and f t¢ a ( ) are linearly increasing with time: f H (t)∼ηTt at t .
The function t t t t , , Taking into account the integrals, such as .
where g α is shown in (37). We also notice that Taking into account that the dissipative function g α (t) goes to zero at t  ¥, we obtain the steady-state expression for the function t t , , D We presume that all matrix elements of qubit operators are taken at time t.
For the case of qubits coupled to environments described in section 3.5, we obtain the following expression for the bath correlator (D27): We notice that the polaron shift (D35) contains contributions of both, LF, ε L , and HF, ε H , parts of the bath reorganization energy since ε α =ε L +ε H . The polaron shift introduced in equations (5) and (6) which rapidly decays during the annealing process according to equation (D17). Therefore, the cross-correlators (D37) between diagonal and off-diagonal operators of the bath give no contribution to the evolution equation (78). áL L ¢ ñ ( ) ( ) As the next step in the derivation of master equations, we calculate the correlator of system operators in the right-hand side of equation (E6). Of interest is when the moments of time t and t¢ are separated by the short time interval τ mn such that For equation (E6) the parameter mn t corresponds to the correlation time given by (D36). It follows from (75) that the evolution of the operator Λ mn is quite slow. The rate of this evolution is determined by the bath operators, such as Q I km and Q nk I , which are proportional to the off-diagonal elements of qubit Pauli matrices, km s a and nk s a , and also to the off-diagonal terms such as T km and T ; nk see equations (46) , (48), (64), (68) for definitions. At first glance, this fact allows us to ignore the variation of t nm L ¢ ( ) in time during the interval τ mn . In this case, the correlator of system operators can be easily calculated: which follows from (75). We show that Markovian fluctuations of the bath characterized by zero-frequency susceptibility χ α (0) contribute to the right-hand side of equation (E9). This contribution is significant because it is proportional to the sums k n nk 2 s å a ¹ | | given by equation (E8). We choose the system basis where the offdiagonal elements, such as nk s a , are small, whereas the diagonal elements n s a can take any values from the interval [−1, 1]. The components of the Pauli matrix z s a are defined by equation (48). Notice that the Markovian contribution to the right-hand side of equation (E9) can be traced without resorting to the perturbation theory in the system-bath interaction. In the process, we drop perturbative terms, which are proportional to individual matrix elements of the Pauli matrices and also to the off-diagonal elements T mñ of the system Hamiltonian H S . For this reason, in the bath operators Q I mn involved in (E9) and defined by equation (76), we assume that From here on, matrix elements, such as km s a , and system operators, such as k n ñá | |, are taken at time t. According to the Wick theorem [45], in equation (E10) we pair the Gaussian bath operator Q α (t) with the other operators, which contain bath variables. Notice that, as follows from the definition in (57), pairings of Q α (t) in (E10) with the matrices t k  ( ) † and t m  ( ) give rise to diagonal matrix elements, such as k s a and m s a , and, therefore, to products such as km k s s a a and . Taking pairings of Q α and the operators U I † and U I , we obtain The bath correlator K t t K t t , 1 1 = a a ( ) ( ) is defined by equation (34). The Markovian part of this correlator has a sharp peak at t=t 1 . Its contribution to the function t t , mnk X ¢ a ( )is described by equation (E18) where in all functions of t 1 we have to put t 1 =t and, after that, remove these functions from the integrals over time. We notice that operators t kn L ( ) (73) and t k l  ¢ ¢ ( ) (E16) taken at the same moment of time commute at any sets of indexes. This fact allows us to calculate the same-time products of the system operators with relations outlined in (E7). In particular, we have E.3. Equations for system operators nn áL ñ As the next step, we substitute the correlation functions (E26) to equation (E6) taking into account that t t 1 ¢ = .
In the process, the system operator t mm 1 áL ñ ( ) is replaced by t mm áL ñ ( ) since this variable is practically unchanged