Diabatic errors in Majorana braiding with bosonic bath
Abstract
Majorana mode based topological qubits are potentially subject to diabatic errors that in principle can limit the utility of topological quantum computation. Using a combination of analytical and numerical methods we study the diabatic errors in Majoranabased topological Yjunction that are coupled to a Bosonic bath in the Markovian approximation. From the study we find analytically that in the absence of a bath, the error rate can be made exponentially small in the braiding time only for completely smooth pulse shapes. Thus, pristine topological systems can reach exponentially small errors even for finite braid times. The presence of a dominantly dissipative Markovian bath is found to eliminate this exponential scaling of error to a powerlaw scaling as with being the braiding time. However, the inclusion of relaxation imroves this scaling slightly to go as . Thus, coupling of topological systems to Bosonic baths can lead to powerlaw in braiding time diabatic errors that might limit the speed of topologically protected operations.
I Introduction
NonAbelian states are the building blocks of topological quantumcomputation as these states carry nonlocally encoded quantum information immune to decoherence Chetan Nayak, Steven H. Simon, Ady Stern, Michael
Freedman, and S. Das Sarma. (2008). One class of such nonAbelian states are Majorana zero modes (MZMs) that emerge naturally as the groundstate of (effective) pwave superconductors Alexei Kitaev (2000); C. W. J. Beenakker (2013).
MZMs (and nonAbelian states in general) appear as degenerate ground states of a stationary Hamiltonian protected by a gap wherein its simplest form involves appearance of two localized MZMs macroscopically separated in one dimensional geometry (a nanowire in an experimental setup) forming a twofold topologically protected degenerate ground state subspace Alexei Kitaev (2000); C. W. J. Beenakker (2013). In such systems quantum information can be encoded in the fermion parity of the degenerate ground state which depend on the occupation/nonoccupation of nonlocal Dirac fermion mode formed out of the pair of localized MZMs. While the gap stabilizes (to exponential accuracy) the ground state subspace at finite temperature, the nonlocally encoded quantum information is protected against environmental decoherence Chetan Nayak, Steven H. Simon, Ady Stern, Michael
Freedman, and S. Das Sarma. (2008); Cheng et al. (2012).
Concrete theoretical proposals to realize MZMs in solidstate systems Roman M. Lutchyn, Jay D. Sau, and S. Das
Sarma. (2010); Yuval Oreg, Gil Refael and Felix von
Oppen. (2010) have led to substantial experimental efforts to realize MZMs in a laboratory. Encouraging experimental progress in this direction in the recent years V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P.
A. M. Bakkers and L. P. Kouwenhoven. (2012); Leonid P Rokhinson, Xinyu Liu, and Jacek K
Furdyna. (2012); A. D. K. Finck, D. J. Van Harlingen, P. K. Mohseni, K.
Jung, and X. Li (2013); Anindya Das, Yuval Ronen, Yonatan Most, Yuval Oreg,
Moty Heiblum and Hadas Shtrikman. (2012); H. O. H. Churchill, V. Fatemi, K. GroveRasmussen, M.
T. Deng, P. Caroff, H. Q. Xu, and C. M. Marcus. (2013); Albrecht et al. (2016); Deng et al. (2016); Zhang et al. (2018) has motivated new set of theoretical proposals that outline basic architecture of quantum gates to be realized following specific protocols to move and manipulate MZMs in order to braid and exchange a set of MZMs with each other Alicea et al. (2011); Sau et al. (2011); Van Heck et al. (2012); Hyart et al. (2013).
Despite the excitement regarding the possibility of exponentially small error in topological qubits realized from MZMs, most of the studies of the protection of MZMs has been in the equilibrium phase. A less studied question is the protection of quantum information stored in MZMs to diabatic errors resulting from the finite speed of operations Cheng et al. (2011); Karzig et al. (2013); Scheurer and Shnirman (2013); Knapp et al. (2016). Specifically, one might worry that since the topological phase of matter is a ground state property, finite gate speed might have an effect of taking the system out of the ground state in a way that introduces errors. This issue has been raised in some studies of the dynamics of braiding Knapp et al. (2016) that suggest the use of a measurement based protocol as a possible way to avoid such diabatic errors. In contrast keyboard based braiding protocols Alicea et al. (2011) seem to reduce some of the diabatic errors and find error rates that scale exponentially in the rate of the process Bauer, Bela and Karzig, Torsten and Mishmash, Ryan V. and Antipov, Andrey E. and Alicea, Jason (2018); Pedrocchi and DiVincenzo (2015). Another potentially critical ingredient in MZM braiding is the interaction of the MZMs with a Bosonic bath. While the effects of Bosonic bath on stationary MZMs have been studied Cheng et al. (2012); Goldstein and Chamon (2011); Rainis and Loss (2012); Konschelle and Hassler (2013), its effect of diabatic braiding have been not been investigated until recently Pedrocchi and DiVincenzo (2015) (also see Ref. Pedrocchi et al. (2015)) where the combined effect of (Bosonic) environmental noise and diabatic braiding is found to result in powerlaw scaling errors even for the keyboardlike braiding schemes Pedrocchi and DiVincenzo (2015).
Our work focuses on the question of diabatic errors in nominally the simplest braiding protocol in a Yjunction type Majorana architecture Karzig et al. (2016) that is coupled to a Bosonic bath (see Fig. 1). The braiding protocol is based on the tunneling induced transport of MZMs Sau et al. (2011) where the splitting between pairs of MZMs is used to exchange the decoupled MZMs (see Fig. 1) that are used to store quantum information. Fig. 1 describes a particular example of a braiding protocol that leads to exchange of and . Since the system is isolated (apart from the Bosonic bath) the total Fermion parity of the system is conserved throughout the braiding protocol. Using the fact that at least one of the MZMs are isolated from the rest of the MZMs at any time in the protocol, it can be shown Karzig et al. (2016) that the two Fermion parity states remain topologically degenerate. Both the initial and final state leave decoupled from everything else. Therefore, the conservation of Fermion parity equates the conservation of the MZM parity , which is used to store quantum information, to the conservation of the of the Fermion parity of the coupled pair . The latter is associated with excitations of the system, so that the bitflip error is directly related to the rate of exciting the system out of the ground state into the excited state. The ground and excited states are the only states in a fixed Fermion parity sector, so that the problem of determining the bitflip error maps to the excitation probability of spin in a timedependent magnetic field. Later we will show that the system in Fig. 1 is topologically protected against dephasing errors in the Fermion parity basis. Therefore the problem of bitflip errors in the system in Fig. 1 can be mapped entirely to the relatively wellstudied problem of diabatic errors for a spin in a timedependent magnetic field (see for e.g. Childs et al. (2001); Yao et al. (2007)).
Despite the mapping of the system in Fig. 1 to a spin in a magnetic field, the topological nature of the setup leads to certain unique features when considering interactions of the system with a Bosonic bath. Microscopically, we assume the Bosonic bath couples to the spin as a magnetic field noise similar to the classic spinBoson model Leggett et al. (1987); Weiss (2012). To simplify our treatment, we assume that the coupling to the Bosonic bath is weak compared to temperature (much smaller than the gap) so that the bath can be modeled within the Markovian approximation using the Davies prescription Brian Davies (1974); Breuer and Petruccione (2002). However, unlike a conventional spin, the vanishing of a component of the magnetic field also implies a vanishing of the noise. This is because such a vanishing of a component of the magnetic field is assumed to occur because of isolation of one of the MZMs from the rest of the system. This leads to conservation of the associated MZM operator, which in turn leads to conservation of the associated excitation. Specifically, this means that for the setup in Fig. 1, the Bosonic bath is forced to decouple from the topologically protected quantum information at the end of the process. However, this also means that any excitation generated during the dynamics of the effective magnetic field cannot relax away at the end of the process. This is in contrast to the dynamics in the spinBoson model, where the system in contact with a zerotemperature Bosonic bath would always relax back to the ground state once the magnetic field becomes static at the end of the process. The absence of such relaxation leads to the a finite excitation probability in the braiding setup in Fig. 1, which leads to the possibility of the bitflip error studied in this work.
Motivated by the mapping of the setup of Fig. 1 to a spin in a magnetic field, in this work we study the probability of excitation of a spin in a timedependent magnetic field that is coupled to a Bosonic bath. The coupling to the Bosonic bath is assumed to be small enough so that it can be studied within the Markovian approximation using the Davies prescription Brian Davies (1974); Breuer and Petruccione (2002). This leads to a timedependent master equation which is further reduced to a Bloch equation describing spin1/2 particles with relaxation and dephasing Breuer and Petruccione (2002); Preskill (1999). To simplify the parameter space of possibilities, we assume that the temperature is low enough so that thermal excitation can be ignored. We show analytically that the excitation probability in such as spinsystem (corresponding to the error rate in Fig. 1) vanishes only polynomially as the time, , within which the braid in completed.
The paper is arranged as follows. In Sec.II, we map the time dependent braiding Hamiltonian describing Fig. 1 to an effective twolevel spin1/2 system coupled to a bath described by a Bloch equation. In Sec.III, following the analytical method developed in Ref. Hagedorn and Joye (2002), we introduce the framework to calculate diabatic corrections to timeevolved groundstate and study diabatic corrections in the absence of bath coupling and present an analytical formula to calculate diabatic correction for a purely dephasing bath (i.e. in presence of bathinduced dephasing but complete absence of absence of relaxation). In Sec.IV we study our system coupled to a general bath where both dephasing and relaxation mechanisms are present. We summarize our results in Sec.V and provide details of our calculations in the appendix.
Ii Braiding Hamiltonian and bosonic bath
In this work, we describe the system coupled to the bath by a Master equation that describes the effect of the bath on the timeevolution of the density matrix through a sequence of socalled jump operators Preskill (1999). Below we first show how the Majorana braiding system shown in Fig. 1 can be described by an effective spin1/2 Hamiltonian in a fixed Fermion parity sector. We then write down the jump operators in the Master equation that describe the coupling to the thermal bath.
ii.1 Spin representation of braiding Hamiltonian
The ideal system (i.e. apart from the bath) comprise four Majoranas labeled shown in Fig. 1. The unitary time evolution generated by the timedependent braiding Hamiltonian over a time cycle results in exchange of two specific Majoranas. Such a braiding Hamiltonian is given by
(1) 
with . The Hamiltonian of the system coupled to the bath conserves Fermion parity . We use the lowest energy state of each parity to define the qubit. Choosing describes coupling among Majoranas such that there is atleast one uncoupled Majorana for at any time ensures a topological degeneracy Karzig et al. (2016). To understand this, we label the Majorana as the isolated Majorana at a particular time and check that
(2) 
If we define as the instantaneous ground state with even parity , it follows from Eq. 2 that the is an eigenstate of , which is degenerate with and has odd Fermion parity (i.e. ).
The time evolution of the Majoranas is governed by the Hamiltonian where the timedependent magnetic field shows a three step time dependence shown in Fig. 2. These steps are defined by the vanishing of one of the components and , respectively, which implies the isolation of the corresponding Majorana that is required to guarantee the topological degeneracy. Using the relations Eq. 2, we can use the isolation of (i.e. ) at time , to constrain the timeevolution of a pair of states related by:
(3) 
at time . Eq. 2 lets us extend this relation from to the entire first segment .
To determine the timeevolution for the rest of the braid we must assume that the dynamics is slow enough so that the Majorana part of the system remains in the ground state. Most of the paper that follows is dedicated to determining the conditions for validity of this assumption. So we use this assumption for the rest of this subsection only. Using this assumption at time , where are the only coupled Majoranas, the parity of the gapped wires maybe constrained as
(4) 
Since , this also implies
(5) 
Repeating the arguments from Eq. 3,4,5 two more times leads to the relation
(6) 
at . Comparing this equation to Eq. 3, we realize that within the lowenergy subspace spanned by , the unitary timeevolution (ignoring an overall phase) can be written as , which is the standard representation Chetan Nayak, Steven H. Simon, Ady Stern, Michael Freedman, and S. Das Sarma. (2008) for exchanging the Majoranas . Note that the only assumption that was crucial for this argument was that the Majorana part of the system remains in the ground state. In principle, this allows us to consider to be a wavefunction of a Majorana system interacting with a Bosonic bath that we will consider later. Therefore, we can conclude that the qubit error (both bitflip and dephasing) probability for the setup in Fig. 1 is given by the probability of excitation of the system out of the ground state in a fixed parity sector.
The flexibility of being able to focus on the excitation probability in a fixed parity sector allows us to map the Majorana dynamics problem to that of a spin system. To see this we note that the Majorana operators may be expressed in terms of Pauli matrices as:
(7) 
(where and are two sets of Pauli matrices) so that Eq. 1 can be rewritten as,
(8) 
The Hamiltonian commutes (i.e. ) with the fermionic parity operator . Within the parity sector, the dynamics of the system may be described by the reduced 2level Hamiltonian,
(9) 
Derivative discontinuities in the timedependence of are known to introduce diabatic errors in Majorana qubits that scale as a powerlaw in Knapp et al. (2016). To avoid such diabatic errors the profile for in Fig. 2 is chosen to be a piecewise continuous functions where all derivatives are continuous (i.e. ). The specific form in Fig. 2 that accomplishes this level of smoothness is given by
(10) 
where
(11) 
and
(12) 
is the dimensionless time parameter.
ii.2 Master Equation
Here we describe the Master equation that we use to model the interaction of the system with a parityconserving thermal bath within the BornMarkov approximation. As discussed such a bath can only generate errors by flipping the system to the excited state of the spin Hamiltonian in Eq. 9. Therefore, similar to the previous subsection, it suffices us to derive a Master equation for the spin1/2 system.
The total Hamiltonian describing the system interacting with a bosonic bath is,
(13) 
where, are system Hamiltonian, bosonicbath Hamiltonian and Hamiltonian describing systembath coupling, respectively. In general, the systembath Hamiltonian have the following form,
(14) 
where and are system and bath operators, respectively, satisfying and .
Following the prescription introduced by Davies Brian Davies (1974), under the assumptions of weak systembath coupling (Born approximation) and memoryless bath (Markov approximation: bath correlation time vanishes on system time scales), the master equation describing the time evolution of density matrix of the system can be expressed in the Lindblad form Brian Davies (1974); Breuer and Petruccione (2002) as
(15) 
where the bath term is written as:
(16) 
The above bathinduced dissipation term is written in terms of projections of the system operators , referred to as ”jump operators” that are written as
(17) 
with being the projection operator to an eigenspace spanned by eigenstates of with eigenvalue . The coefficients in Eq. 16 are correlators of the bath operators that are written as
(18) 
where, denotes thermal average, . Depending on the value of , the jump operators are classified as excitations (), relaxation () or dephasing ().
The bath operators couple to each of the components of (where ) in a way where the bath coupling to a particular wire in Fig. 1 vanishes when the corresponding vanishes. This is essential to represent the fact that the vanishing results from decoupling of from . To satisfy this condition, we choose the system operators for the bath coupling as
(19) 
Using Eq. 17, the jump operators corresponding to the above choice of operators are expressed as,
(20) 
where,
(21) 
being the timeindependent systembath coupling strength with and and denoting the instantaneous ground and the excited states of , respectively.
We wish to focus on only those batheffects that arise due to timedependence of the Hamiltonian. In other words we want to explicitly avoid any bathinduced effect that do not vanish in the adiabatic limit provided the system is initialized in the ground state. The only such environmental effect is thermal excitation associated with (energy gap above ground state is 2 for ) operator. We explicitly set temperature to zero to completely suppress these thermal excitations or equivalently,
Since the braiding process involves sequence of three identical clockwise rotations of along and axes respectively, as shown in Fig. 1, we focus on the first sequence where is restricted to XY plane without loss of generality, setting . Influence of the bath is captured by the strength of the systembath couplings (captured by and ), the relaxation strength governed by (since the gap in the system is 2) and dephasing strength governed by . Since the timedependence of the Hamiltonian does not affect the gap in the system and are fixed parameters determined by entirely by the microscopic properties of the bath Hamiltonian. From here onwards, for the sake of brevity, we relabel the relaxation strength, and the dephasing strength, with new symbols, and .
ii.3 Bloch equation
The density matrix appearing in the master equation (Eq. 15 can be parametrized in terms of a Bloch vector defined by
(23) 
Rewriting the Master equation (Eq. 15) in terms of leads to the socalled Bloch equation
(24) 
where,
(25) 
with,
(26) 
and . Note that the time derivatives in the above equation refer to the rescaled time .
The Bloch equation can be more compactly expressed in a matrix form as,
(27) 
with, being a null vector of , and where,
(28) 
Iii Adiabatic expansion for Bloch Equation
As discussed in Sec. I, the diabatic error for the Majorana qubit in Fig. 1 is related to the probability of transition out of the ground state. Such a transition can occur in any one of the three steps of the braiding protocol, so that the order of magnitude of the error can be estimated from the transition probability in any one step. We consider the first step of the braiding process where and calculate the diabatic error incurred in the process which involves rotating along axis by in the clockwise direction from the initial orientation along to final orientation along . The system is initialized in the ground state of to which the corresponding Bloch vector is . Given this initial condition, the solution of the Bloch equation (i.e. Eq. 27) in the extreme adiabatic limit (i.e. ), is written as:
(29) 
At nonzero , the timedependent solution of the Bloch equation 27, referred to as timeevolved Bloch zerovector (TBZV), differs from the adiabatic limit . The magnitude of the difference at
(30) 
quantifies the diabatic error for the braiding.
Let us now consider the solution of the Bloch equation Eq. 27 for infinitesimal . We first focus on the case of vanishing relaxation (i.e. ) where the Bloch equation reduces to
(31) 
The finite relaxation situation will be dealt with in Sec. IV using a slightly different technique. The solution of the Bloch equation in the above form admits an adiabatic series expansion written as:
(32) 
Substituting the above ansatz in Eq. 31 and equating both sides of the equation at each order in we formally solve for order correction to for (see App. A for details),
(33) 
with given by,
(34) 
The matrix is singular so that it is crucial to note that Eq. 31 implies that . Additionally this implies that the generalized inverse can be defined so that is orthogonal to . Motivated by this, it is convenient to split the expansion (Eq. 32) into two parts
(35) 
where, and and are expanded as:
(36) 
Taking repeated derivatives of Eq. 33(see App. A for details) we can show that vanishes when all derivatives of vanish, at . This implies that vanishes as well only when the magnetic field goes to 0 smoothly at , which is a condition that is ensured by the specific timedependence of the magnetic field chosen (i.e. Eq. 10). More precisely, considering the expansion for (Eq. 36) to a finite order , one can show that
(37) 
for every integer . This shows that vanishes faster than any powerlaw in the adiabatic limit. A numerical analysis shows that
(38) 
which is similar to (but slower than) exponential scaling with .
In the absence of any bath (i.e. ) the dynamics of is unitary and . Thus the bound on (Eq. 38) also implies a bound on the diabatic error for the isolated qubit given by
(39) 
This is completely consistent with the numerical results in Fig. 3 and also with those derived in Ref. Hagedorn and Joye (2002), which suggests an exponentiallike dependence of diabatic error on .
In the presence of dephasing , the unitarity constraint does not apply, but one can use the essential vanishing of (Eq. 37) in combination with Eq. 35 to write
(40) 
Evaluating the recursive relation Eq. 34, one finds (see App. A for details) a nonvanishing contribution to that is written as
(41) 
This leads to a powerlaw diabatic error estimate,
(42) 
Note that such error term (as well as the other powers of ) dominate over estimate for the component that was ignored but would vanish in the absence of a bath () restoring the exponential estimate for the isolated system situation.
In Fig. 4, we plot as a function of for different systembath coupling strengths and compare it against numerical curve obtained by solving Eq. 31. We find that numerically exact calculation of validates the asymptotic (large ) limit of the diabatic error (Eq. 42) obtained for truncation at firstorder in epsilon. We find that the scaling region that follows Eq. 42 begins rather abruptly above a value of that depends on the dissipation strength. The behavior below this threshold coincides with the behavior of the dissipationless system shown in Fig. 3. Thus, the strength of the dissipation determines the critical braiding rate above which the exponential behavior of diabatic error (Eq. 39) crosses over into powerlaw error described by Eq. 42.
Iv General systembath coupling
We now consider the error in the presence of finite relaxation. In order to solve for in Eq. 24 for the case of finite we switch to a new vectorvariable ,
(43) 
where, is the rotation matrix such that ( being the initial condition on ) given by
(44) 
From Eq. 43 it follows, the initial condition on is . Since is the solution to the Bloch equation in limit, must have perturbatively small norm in .
The equation of motion for is found to be,
(45) 
On the right hand side of the equation, the operator acting on can be diagonalized along the basis unit vectors with being the corresponding eigenvalues where we have defined,
(46)  
(47) 
Representing in this basis set,
(48) 
leads to the following coupleddifferentialequation,
(49a)  
(49b) 
Finite relaxation ensures that at sufficiently long time (or sufficiently small ), the memory of the initial conditions at are exponentially suppressed. In this case, Eq. 49a suggests is so that, to lowest order in ,
(50) 
This implies as (see Eq. 11). Therefore, we conclude . Substituting expressions in Eq. 49a, it can be solved for resulting in
(51) 
where we have defined,
(52) 
The behavior of in the largeT limit is analyzed using saddlepoint method detailed in the App. B to obtain the following asymptotic form for the diabatic error,
(53) 
Note that Eq. 51 reduces to Eq. 42 in absence of relaxation, i.e. , showing that the approach in this section is not inconsistent with that in the previous section for the case without relaxation. An examination of the exponential factor in Eq. 51 tells us that for the case where relaxation is weak compared to dissipation (i.e. ), we should expect an intermediate regime of where might still be ignored so that the error would scale as as concluded in the previous section.
Therefore as is increased we expect the error rate to crossover from the isolated system limit (Eq. 39) to the dominantly dissipative system(Eq. 42) to the asymptotics with relaxation (Eq. 53). The numerical plot of the error rate in the general case, Fig. 5 indeed shows this expected pattern of crossovers.
These expectations turn out to be quite accurate as seen from the numerical simulations in Fig. 5, which shows the numerical estimate of vs total time in loglog scale for dephasing strength , and relaxation strength , both set equal to 1. It is interesting that the crossover behavior expected for the limit seems to apply even to chosen for Fig. 5. Finally, we see in Fig. 5 that as we decrease the bath coupling strength , the crossover scale modes to higher as expected. While Eq. 52 predicts the correct crossover behavior, it turns out to be off by a independent scale factor in the preasymptotic dissipation dominated regime. We have remedied this by adapting the analysis in Sec. III to finite (but small ) as described in App. B. This analysis leads to a somewhat different approximation for which is written as
(54) 
This more exact theoretical form, which is the analytic approximation used in the dashed lines in Fig. 5, can be seen from the figure to fit the numerical results very well in both the dissipative and relaxation dominated regimes.
V Summary
In this work we have studied the diabatic error rate as a function of braiding time (in units of inverse gap) of the Yjunction braiding protocol (shown in Fig. 1) by mapping it to the problem of excitation probability of a spin in a timedependent magnetic field. Consistent with previous work Hagedorn and Joye (2002) we find (see Eq. 39) that one can reduce diabatic errors Knapp et al. (2016) in isolated Majorana systems to be exponentially (scaling as ) small in the braiding time (in units of inverse gap) by choosing the timedependence of the Hamiltonian to be completely smooth (including the beginning and end of the protocol). This analytic result explains the recent obeservation of exponentially suppressed diabatic errors apparent from numerical simulation of some braiding protocols Bauer, Bela and Karzig, Torsten and Mishmash, Ryan V. and Antipov, Andrey E. and Alicea, Jason (2018) In fact, while the requirement of a completely smooth timedependent Hamiltonian of the form of Eq. 8 (also see Eq. 10) seems finetuned, it is quite natural in a topological system where MZM splitting is tuned either by chemical potential Alicea et al. (2011); Sau et al. (2010); Bauer, Bela and Karzig, Torsten and Mishmash, Ryan V. and Antipov, Andrey E. and Alicea, Jason (2018) or by screening charging energy through tunable Josephson junctions Hyart et al. (2013); Clarke et al. (2016). In both these cases, the tunneling is exponentially suppressed as one tunes the Majorana wire deep into the topological phase or introduces a strong Josephson coupling between the Majorana wire and a bulk superconductor. The main focus of our work is to consider the effect of interaction of the topological system with a Fermion parity conserving Bosonic bath, such as phonons or plasmons on the diabatic error rate for braiding. Similar to previous work Pedrocchi and DiVincenzo (2015), we assume the bath to be weakly coupled so that we can treat the bath within the Markovian approximation. We describe the dynamics of this system within the Markovian approximation by a Bloch equation. We find that the coupling to such a Bosonic bath generically changes the asymptotic of the diabatic error from exponential in the braiding time to power law (i.e. error scaling as (see Eq. 53). as expected, the general result including relaxation leads to much lower excited state probability than in the absence of relaxation. We find from a controlled analytic solution that the error rates in the presence of dephasing from the bath but no relaxation decreases the slowest scaling as .
In addition to the asymptotic forms we study the dependence of the diabatic error rate on the systembath coupling strength through direct numerical simulations. For purely dephasing bath (i.e. no relaxation), we find (see Fig. 4) that the dephasing strength determines the crossover from the scaling error rate that is expected in the absence of a bath to the scaling of the error rate expected for a purely dephasing bath. Increasing the dephasing strength leads to the crossover occuring at shorter braiding time , in turn leading to error rates that increase with increasing dephasing. As seen from our results in Fig. 5, adding relaxation in addition to dephasing leads to an additional relaxation strength dependent crossover timescale beyond which the error rate changes its scaling from to .
In our work, we have ignored finite temperature effects despite the fact that a finite temperature is needed to justify the Markovian approximation. Such finite temperature effects can be expected to be negligible for temperatures substantially below the gap. Moreover, such finite temperature excitations are expected to only increase the excitation probability and thus the error rate. In summary, we find that while a Bosonic bath does not directly interfere with topological protection of quantum information, finite (but low) temperature Bosonic baths can lead to powerlaw in time diabatic errors that might limit the speed of topologically protected operations.
This work was supported by the Microsoft Station Q, NSFDMR1555135 (CAREER), JQINSFPFC (PHY1430094) and the Sloan research fellowship.
Appendix A Diabatic expansion of Bloch vector
Expanding the Bloch vector in powers of ,
(55) 
we seek solution to
(56) 
for initial condiction where, and ( being the matrix representation of operation) with and being the timedependent functions defined by Eq. 25. The presentation here follows the work of Hagedorn et. al. in Ref. Hagedorn and Joye (2002).
, being the antisymmetric matrix representation of operator has eigenvalues . Consequently, the instantaneous eigenvalues of are given by with,
(57) 
Denote and thereby, is the zero eigenvector with eigenvalue. can be inverted in the eigenvector subspace with nonzero eigenvalues,
(58) 
Expanding and substituting in , we get,
(59) 
with evaluated using the condition which follows from ,
(60) 
Now we show that the series expansion is welldefined in small epsilon limit. Consider the partial sum of the series expansion,
(61) 
If the series expansion is welldefined, the partial sum (as defined above) must converge to the actual solution . Let the actual solution . Consider,
(62) 
where we have used which follows from . Now,
(63) 
where we simply used the definition given in Eq. 61 and the relation given in Eq. 59 to arrive at the second line above. Using the above relation we get,
(64) 
Therefore it follows from Eq. 62,
(65) 
We conclude that converges to the actual solution provided and are bounded. We refer the reader to Ref. Hagedorn and Joye (2002) for the proof of boundedness of in absence of bath. It seems likely that a similar proof holds for boundedness of in presence of bath. In the limiting , is unitary, so clearly when for , is bounded by 1. Without speculating about case, we restrict our following discussion to that corresponds to provided .
Our goal now is to arrive at a bound for the diabatic error defined by,
(66) 
where is the Bloch vector that corresponds to instantaneous zero eigenvector of at . Note as consistency check that when follows from the series expansion of , Eq. 34. Since we have shown that the series expansion of is well defined, we can replace in by its corresponding series expansion. However to make progress towards arriving at a bound for we need a useful result stated and proved below.
We show that if all derivative of ( derivative denoted by , for some then where is any vector such that . The proof follows in three steps:

We first show that .
Let be a positive integer. . Therefore, it must be . 
Next we show that .
Again, let be a positive integer. . Therefore, it must be . 
Finally we prove our original assertion, .
We will prove this assertion by induction. Assume . Now,Using the the induction hypothesis it follows,