Analytical model of non-Markovian decoherence in donor-based charge quantum bits

# Analytical model of non-Markovian decoherence in donor-based charge quantum bits

F Lastra, S Reyes and S Wallentowitz Facultad de Física, Pontificia Universidad Católica de Chile, Casilla 306, Santiago 22, Chile
###### Abstract

We develop an analytical model for describing the dynamics of a donor-based charge quantum bit (qubit). As a result, the quantum decoherence of the qubit is analytically obtained and shown to reveal non-Markovian features: The decoherence rate varies with time and even attains negative values, generating a non-exponential decay of the electronic coherence and a later recoherence. The resulting coherence time is inversely proportional to the temperature, thus leading to low decoherence below a material dependent characteristic temperature.

###### pacs:
03.67.Lx, 73.63.Kv, 03.65.Yz

## 1 Introduction

Implementations of quantum bits (qubits) in solid-state systems represent some of the most promising candidates for quantum-computation devices. The constituent two-level systems have been experimentally realized in a variety of systems, including quantum dots [1], impurities in diamond [2], and superconductors [3] among others. The implementation of solid-state quantum-information devices has been object of increasing interest in the last decade [4, 5, 6].

Like other types of implementations, solid-state systems also suffer from the inevitable interaction with their environment, which leads to a loss of coherence of individual qubits and of entanglement between them. Whereas the decoherence of many physical systems can be well approximated by a Markovian process, a few prominent examples exist where the coupling of the system with its environment reveals a memory effect, leading to non-Markovian behavior and the possibility of recoherence at later times. Among these are spin echo [7] — which has been recently observed with phosphorus donors in silicon [8], photon echo [9, 10], and collapse [11] and revival [12, 13, 14] in the spin-boson Jaynes–Cummings interaction [15]. Collapse and revival may occur in various physical systems, such as with atoms interacting with cavity electromagnetic radiation [16], in the vibronic dynamics of electromagnetically trapped ions [17], or in molecular vibration [18, 19]. However, also other systems may show a non-Markovian dynamics in the decay of their coherence. In particular, numerical calculations have shown that donor-based charge qubits interacting with phonons reveal a non-Markovian decay of electronic coherence [20].

Donor-based charge qubits consist of a single electron shared by two impurities in the semiconductor host. Besides being potentially scalable, their manufacturing is largely facilitated by the technological base provided by the existing semiconductor industry. One of the main sources of decoherence in this system is the coupling of the electron to acoustic phonons. This decoherence mechanism has been studied numerically for temperatures [20] and analytically for [21]. In the former reference, a non-exponential decay of coherence has been numerically obtained, which indicates a non-Markovian behavior of the electron-phonon coupling. On the other hand, in Ref. [21] decoherence has been derived even for zero temperature, where the bosonic phonon modes are in the vacuum state.

Non-Markovian coupling of a system to a reservoir, has received increasing theoretical interest in recent years [22, 23, 25, 24, 29, 30, 31, 27, 26, 28]. An important conceptual issue is to obtain a solution when decay rates attain negative values. In the literature, there exist numerical methods to take into account negative decay rates. These comprise quantum trajectory simulations in extended Hilbert spaces [32, 33, 34, 35, 36] and alternative quantum-trajectory methods [37, 38]. However, for presently studied physical systems, negative decay rates have been predicted only for extremely short time scales.

We will show in this work that donor-based charge qubits, subject to the interaction with acoustic phonons, do naturally show this feature on time scales, that may be observable in experiment. Furthermore, the model developed here contradicts the result of Ref. [21], where decoherence is obtained even for zero temperature. From our work, the coherence time is obtained as inversely proportional to temperature, leading to vanishing decoherence at zero temperature.

## 2 Model of Phonon Scattering with a Localized Electron

The Hamiltonian of the system under study is composed of three parts,

 ^H=^He+^Hph+^He−ph, (1)

where the electronic part is in principle composed by the Bloch electrons plus the electrons bounded to the impurities. The quantized Fermionic field for electrons with spin projection can be expanded as

 ^Ψς(r)=∑n,q^cn,q,ςψn,q(r)+∑a,n^da,n,ςϕa,n(r), (2)

where and are the Bloch functions of itinerant electrons and the wavefunction of the electron localized at the donor site , respectively. The operators and annihilate an electron with spin projection in the Bloch state and in the localized state , respectively. In our case the impurities consist in pairs of donors (e.g. P) embedded in a semiconductor substrate (e.g. Si). The full electronic Hamiltonian reads

 ^He = ∑n,q,ςEn,q^c†n,q,ς^cn,q,ς+∑a,n,ςℏεa,n^d†a,n,ς^da,n,ς (3) +12∑n,n′ℏΔn,n′(^d†+,n^d−,n′+^d†−,n′^d+,n),

where and are the energies of Bloch and localized electron, respectively. Here we added the third term which is responsible for tunneling between localized states at different donor sites with rates .

The phonon Hamiltonian corresponds to the energy of vibrational excitations of the crystal lattice of the host semiconductor (Si). These excitations are bosonic in nature and we can write,

 ^Hph=∑k,σℏωk,σ^b†k,σ^bk,σ, (4)

where the bosonic operator annihilates a phonon with wavevector , polarization () and energy . We will consider only acoustic phonons, as the larger energies of optical phonons make them practically unpopulated at room temperature and below.

### 2.1 Electron-phonon interaction

The electrons couple to the lattice vibrations through the third part of the Hamiltonian,

 ^He−ph=e∫d3r^ρe(r){∑n^un\textperiodcentered∇V(r−Rn)}, (5)

where the sum is performed over the lattice sites. Here is the elementary electric charge and is the electron density. Furthermore, is the Coulomb potential of a single ion of the lattice with and being the equilibrium position and quantized displacement of the ions, respectively.

In terms of the bosonic phonon operators we can write Eq. (5) as

 (6)

where is the volume under consideration and the and sums are over the first Brillouin zone and the reciprocal lattice, respectively. The coupling constant is derived as

 κk,K,σ=ie√Nℏ2m0ωk,σ(k+K)% \textperiodcenteredϵk,σVk+K. (7)

with and being the ionic mass in the unit cell and the number of unit cells in the semiconductor crystal. Here the Fourier coefficients of the ion potential are defined by

 V(r)=1V∑kVkeik\textperiodcenteredr.

Using Eq. (2), the electron density can be written as

 ^ρe(r)=^ρ(B)e(r)+^ρ(BL)e(r)+^ρ(L)e(r), (8)

where the different terms correspond to pure Bloch electrons (B), cross terms between Bloch and localized electrons (BL), and pure localized electrons (L). Inserting Eq. (8) into (6), the second term (BL) will produce phonon-induced transitions between Bloch and localized electrons. However, for the thermal energy of phonons is not sufficient to drive such transitions. Therefore, we can neglect the BL term together with the B term, since under these circumstances the latter does not affect the dynamics of the localized electron under concern. Thus, we are left with the electron density

 ^ρ(L)e(r)=∑a,b∑n,n′^d†a,n,ς^db,n′,ςϕ∗a,n(r)ϕb,n′(r), (9)

For sufficient high temperatures, where the average thermal phonon energy, , is of the order of the electronic excitation energies of individual donor sites, phonon-induced electronic transitions () within the individual donor sites () or between the donor sites () may occur. In this case, resonant phonon scattering occurs that leads to a decoherence mechanism with a strong dependence on the inter-donor distance [39, 40].

In this work we consider a different temperature range, i.e. is much smaller than the electronic excitation energies of individual donor sites. In this case only the electronic ground states of the donor sites can be populated and the density (9) is further simplified to

 ^ρ(L)e(r)=∑a,b^d†a^dbϕ∗a(r)ϕb(r), (10)

where with being the invariant spin projection of the single electron and . Inserting Eq. (10) into the interaction Hamiltonian (6) we obtain

 ^He−ph=1V∑a,b=±∑k,σ∑Kκk,K,σFa,b(k+K)^d†a^db(^bk,σ+^b†−k,σ), (11)

where the form factor is defined as

 Fa,b(k)=∫d3rϕ∗a(r)ϕb(r)eik\textperiodcenteredr. (12)

Given that the rms spread of the localized electron wavefunction is much larger than the periodicity of the semiconductor lattice, i.e. it spreads over several lattice sites, the form factor (12) is non-vanishing only within the Brillouin zone, if . Furthermore, given that the wavefunctions of the two donor sites do not overlap, for . Considering the medium as being approximately isotropic the interaction Hamiltonian (11) can then be simplified to

 ^He−ph=1V∑a∑kκkFa(k)^d†a^da(^bk+^b†−k), (13)

where and where are annihilation operators of acoustic phonons with longitudinal polarization. The simplified coupling constant is

 κk=ie√Nℏ2m0ωkkVk. (14)

with linear dispersion , sound velocity , and with being the Fourier components of the ion potential.

### 2.2 Electronic pseudo-spin

In the basis of the localized-electron states () with and for only a single localized electron being present, i.e. , one may define the electronic pseudo spin- operator,

 ^Sx = 12∑a≠b|a⟩⟨b|, (15) ^Sy = 12i∑a≠ba|a⟩⟨b|, (16) ^Sz = 12∑aa|a⟩⟨a|, (17)

and the electronic identity operator

 ^Ie=∑a|a⟩⟨a|.

With these definitions the sum over donor sites in Eq. (13) is rewritten as

 ∑aFa(k)^d†a^da=^Szfk+^IeFk, (18)

with

 Fk = 12∑aFa(k), (19) fk = F+(k)−F−(k). (20)

The term proportional to the identity operator can be shown to lead to a phonon-induced shift of the electronic ground-state energies of the two donor sites (), see A. It can be accounted for by defining the displaced phonon operator

 ^ak=^bk+βk, (21)

with the displacement being given in Eq. (84) and by redefinition of the electronic energies, see A. In this way the electronic Hamiltonian (3) can be simplified to

 ^He=ℏω0^Sz+ℏΔ^Sx, (22)

where is the transition frequency between the donor-site electronic ground states, modified by the phonon-induced energy shift, cf. Eq. (87) in A, and is the tunneling rate between the ground states of the donor sites. The interaction Hamiltonian (13), on the other hand, becomes

 ^He−ph=∑kℏ^Sz(gk^a†k+g∗k^ak), (23)

where we defined the interaction rate

 gk=κkf−kℏV. (24)

### 2.3 Interaction rate

Assuming a Yukawa potential for the screened ions of charge , Eq. (14) becomes

 κk=iZe2ϵ0√Nℏ2m0ωkkk2+k2s, (25)

with being the inverse screening length of the ionic potential. For semiconductors (Si) and the wavelengths of populated phonon modes are much larger than the screening length, i.e. , so that we may approximate

 κk≈iZe2sϵ0k2s√Nℏωk2m0. (26)

The ground-state wavefunctions at the two donor sites are displaced by the inter-donor distance vector , i.e.

 ϕ±(r)=1√πR3±exp[−|r∓d/2|R±], (27)

where () are the Bohr radii of the -wave orbitals of the two donor sites. Thus the form factors become

 F±(k)=e±ik\textperiodcenteredd/2[1+(kR±/2)2]2, (28)

so that together with Eq. (26) the interaction rate (24) becomes

 (29)

where and are the mass density and unit-cell volume, respectively, and the deformation constant results as .

## 3 Reduced Electronic Dynamics

We assume here that during the free evolution of the system, tunneling is inhibited by either an applied potential barrier between the donor sites, or by the different electronic ground-state energies of the donor sites, , which can be provided for by the application of a DC electric field. Under these circumstances the Hamiltonian is given by

 ^H=ℏω0^Sz+∑kℏωk^a†k^ak+ℏ^Sz∑k(gk^a†k+g∗k^ak), (30)

### 3.1 Electron-phonon eigenstates

We look for the eigenstates of the Hamiltonian (30),

 ^H|E⟩=E|E⟩, (31)

where from the Hamiltonian it is clear that they have the form

 |E⟩→|Em⟩=|m⟩⊗|Φm⟩, (32)

with and with being a solution of the phononic eigenvalue problem

 ^Hm|Φm⟩=Em|Φm⟩, (33)

with the spin-projected phononic Hamiltonian being

 (34)

This Hamiltonian can be diagonalized as

 ^Hm=∑kℏωk^a†k,m^ak,m+ℏmω0, (35)

where we used and discarded a constant energy term. Here the operators

 ^ak,m=^ak+mαk, (36)

are phonon operators with displacements conditioned on the electronic state . Their displacement amplitudes are given by

 αk=gkωk. (37)

The displacement of the phonon operators can be realized by the single-mode displacement operators

 ^Dk(α)=exp(α^a†k−α∗^ak), (38)

as

 ^ak,m=^ak+mαk=^D†k(mαk)^ak^Dk(mαk). (39)

Thus, the Hamiltonian can be written as

 ^Hm=^D†({mαk})^~Hm^D({mαk}), (40)

where

 ^~Hm=∑kℏωk^a†k^ak+ℏmω0, (41)

and where the multi-mode displacement operator is defined as

 ^D({mαk})=∏k^Dk(mαk)=exp[∑km(αk^a†k−α∗k^ak)]. (42)

As a consequence, the phononic eigenvalue problem can be rewritten as

 ^~Hm|~Φm⟩=Em|~Φm⟩, (43)

where the transformed eigenstates are

 |~Φm⟩=^D({mαk})|Φm⟩. (44)

Given that is the Hamiltonian of an infinite set of harmonic oscillators, these eigenstates are identified as those of the harmonic oscillators, i.e.

 |~Φm⟩→|{Nk}⟩, (45)

with being the number of phonons in modes . Consequently the energy spectrum is discrete and reads

 Em→Em,{Nk}=∑kℏωkNk+ℏmω0. (46)

Transforming back to the original frame, the complete eigenstates are therefore given by displaced number states:

 (47)

The general solution of the reduced density operator of the spin system, i.e. traced over the phonons, reads therefore

 ^ϱS(t) = (48) ×⟨{N′k}|^D({m′αk})^D†({mαk})|{Nk}⟩ ×exp[−itℏ(E{Nk},m−E{N′k},m′)].

where we used for the phononic trace and where is the initial density matrix of the complete spin-phonon system. It can be easily seen, that the diagonal matrix elements of the solution (48) are constant,

We note that the initial reduced electronic density matrix elements result from Eq. (48) as

 ⟨m|^ϱS(0)|m⟩ = ∑{Nk}ϱm,{Nk};m,{Nk}, (49) ⟨g|^ϱS(0)|e⟩ = ∑{Nk}∑{N′k}ϱg,{Nk};e,{N′k}f{N′k},{Nk} (50)

and . The non-diagonal reduced matrix element (50) contains the Franck–Condon transition amplitude [41, 42, 43], well known from molecular physics,

 f{N′k},{Nk}=⟨{N′k}|^D({αk})|{Nk}⟩. (51)

The presence of this factor is due to the fact, that the initial density matrix of the complete spin-phonon system, , is defined in the basis , that differs from the basis by the displacements of the phonon state.

### 3.2 Preparation of the initial electronic state

The preparation of a coherent superposition of electronic states and must start from a well defined state, which we assume to be the ground state . This can be experimentally realized by applying a voltage between the donor sites. When thermal equilibrium is reached after the switching on of the voltage, the occupation probability of state is very close to one. The coherent electronic superposition can then be achieved in two ways: Either by switching off the voltage to allow for coherent tunneling between the donor sites, or by applying electromagnetic radiation to perform a coherent Raman transition between states and [44, 46, 45]. To observe the decoherence of the prepared superposition state without the detrimental effects of coherent transitions between the donor sites, the coherent tunneling must be suppressed. This may be achieved either by applying an electric field creating a potential barrier between the donor sites or by switching on a sufficiently large voltage between the donor sites.

As outlined above, the initial state of the complete system is the result of a state-preparation process in which the spin originally starts in its ground state and is coherently transferred into a superposition of ground and excited states. In the preparation process, a fraction of the population is coherently transferred from the ground state to the excited state . However, as in the excited electronic state each phonon mode receives a coherent displacement , cf. the eigenstates in Eq. (47), Franck–Condon type transitions [41, 42, 43] occur for the phonons.

The probability amplitude to end in state is given by

 ψe,{Mk} = ξ∑{Nk}⟨Ee,{Mk}|^S+|Eg,{Nk}⟩ϕg,{Nk} (52) = ξ∑{Nk}f{Mk},{Nk}ϕg,{Nk},

where are the probability amplitudes for starting originally from states . On the other hand the probability amplitude to stay in state is

 ψg,{Nk}=√1−|ξ|2ϕg,{Nk}. (53)

Thus the prepared complete initial state can be written as

 |ψ(0)⟩ = ∑{Nk}(ψg,{Nk}|Eg,{Nk}⟩+ψe,{Nk}|Ee,{Nk}⟩), (54)

which by insertion of Eqs (52) and (53) becomes

 |ψ(0)⟩ = ∑{Nk}ϕg,{Nk} (55) × [√1−|ξ|2|Eg,{Nk}⟩+ξ^D({αk})|Ee,{Nk}⟩]

The corresponding initial density operator is and reads

 ^ϱ(0) = ∑{Nk}∑{N′k}ϕg,{Nk}ϕ∗g,{N′k} × {(1−|ξ|2)|Eg,{Nk}⟩⟨Eg,{N′k}| +|ξ|2^D({αk})|Ee,{Nk}⟩⟨Ee,{N′k}|^D†({αk}) +ξ∗√1−|ξ|2|Eg,{Nk}⟩⟨Ee,{N′k}|^D†({αk}) +ξ√1−|ξ|2^D({αk})|Ee,{Nk}⟩⟨Eg,{N′k}|}

We assume that originally the spin-phonon system started from the spin ground state and the phonons being in a thermal state, i.e.

 ^ϱ=∑{Nk}P{Nk}|Eg,{Nk}⟩⟨Eg,{Nk}|, (57)

with

 P{Nk}=Z−1exp(−∑kβkNk),Z=∑{Nk}exp(−∑kβkNk), (58)

where . Thus we have to substitute

 ϕg,{Nk}ϕ∗g,{N′k}→P{Nk}∏kδNk,N′k, (59)

by which the initial prepared spin-phonon density operator (3.2) becomes

 ^ϱ(0) = ∑{Nk}P{Nk} × {(1−|ξ|2)|Eg,{Nk}⟩⟨Eg,{Nk}| +ξ∗√1−|ξ|2|Eg,{Nk}⟩⟨Ee,{Nk}|^D†({αk}) +ξ√1−|ξ|2^D({αk})|Ee,{Nk}⟩⟨Eg,{Nk}|}.

In the basis of the states (47) the corresponding density matrix reads then

 ϱg,{Nk};g,{N′k} = (1−|ξ|2)P{Nk}∏kδNk,N′k, (61) ϱe,{Nk};e,{N′k} = |ξ|2∑{Mk}P{Mk}f{Nk},{Mk}f∗{N′k},{Mk}, (62) ϱg,{Nk};e,{N′k} = (63)

Furthermore, the reduced initial prepared spin density operator becomes

 ^ϱS(0) = (1−|ξ|2)|g⟩⟨g|+|ξ|2|e⟩⟨e| (64) +ξ∗√1−|ξ|2|g⟩⟨e|+ξ√1−|ξ|2|e⟩⟨g|,

which corresponds to the perfectly pure superposition state

 |ψS(0)⟩=√1−|ξ|2|g⟩+ξ|e⟩. (65)

Thus the presence of the phonons does not prevent a coherent preparation of the initial spin state.

### 3.3 Spin Decoherence

From the general solution (48), the time-dependent off-diagonal matrix elements read

 ⟨g|^ϱS(t)|e⟩ = eiω0t∑{Nk}∑{N′k}ϱg,{Nk};e,{N′k} ×

Inserting the initial state (63) this matrix element becomes

 ⟨g|^ϱS(t)|e⟩ = eiω0tξ∗√1−|ξ|2Z−1 × Tr[^D†({αk})^D({αk(t)})exp(−∑kβk^Nk)],

where

 ^D({αk(t)})=exp(i∑kωk^Nkt)^D({αk})exp(−i∑kωk^Nkt) (68)

with

 αk(t)=αke−iωkt. (69)

The displacements can be joined to obtain

 ^D†({αk})^D({αk(t)}) = ^D({αk(t)−αk}) (70) × exp[−i∑k|αk|sin(ωkt)],

so that Eq. (3.3) becomes

 ⟨g|^ϱS(t)|e⟩ = eiω0tξ∗√1−|ξ|2Z−1 (71) × Tr[^D({αk(t)−αk})exp(−∑kβk^Nk)] × exp[−i∑k|αk|sin(ωkt)].

The thermal averages can be calculated in phase space as

 Tr[^D(α)e−β^N]=∫d2ξeαξ∗−α∗ξe−β|ξ|2∫d2ξe−β|ξ|2=e−|α|2/β. (72)

 ⟨g|^ϱS(t)|e⟩=⟨g|^ϱS(0)|e⟩exp[iω0t−iϕ(t)−∫t0dτγ(τ)] (73)

 ϕ(t)=∑k|αk|sin(ωkt),

and the decoherence rate is given by

 γ(t) = 2∑k|αk|2ωkβksin(ωkt) (74) = 2kBTℏ∑k|αk|2sin(ωkt).

### 3.4 Decoherence rate

The mode sum is evaluated and after integration over the spherical angles of the mode’s wavevector, the rate becomes

 γ(t)=kBTD2(2π)2s3ℏ2ρm∫∞−∞dk{∑aksin(skt)[1+(kRa/2)2]4 +cos(|d+st|k)−cos(|d−st|k)d[1+(kR−/2)2]2[1+(kR+/2)2]2}. (75)

It can be written as

 γ(t)=∑a≠bγa,b(t), (76)

where we defined

 γa,b(t) = ΓTη2a{[13(tτdηa)3+12(tτdηa)2+14(tτdηa)]e−2t/(τdηa) (77) +(η2aη2a−η2b)2[∣∣∣1+tτd∣∣∣+ηa2η2a−5η2bη2a−η2b]e−2|1+t/τd|/ηa

where we defined the relative Bohr radii

 ηa=Ra/d,(a=±),

and the time needed for a phonon to travel between the donor sites, . The temperature dependent rate is defined by

 ΓT=(TT0)ωd,

where is the angular frequency of a phonon with wavelength . The material dependent temperature scale is defined for convenience as

 kBT0=Ndm0s2(ℏωdD)2, (78)

where is the number of unit cells within the volume .

## 4 Discussion

### 4.1 Time-dependent decoherence rate

The initial value of the decoherence rate can be shown to be . Moreover, it is easily shown that it asymptotically decays to zero for large times, . The dimensionless and temperature-independent decoherence rate is plotted in Fig. 1 as a function of the dimensionless time . Here the mean relative Bohr radius has been chosen as , whereas the relative difference of (relative) Bohr radii, with , is varied for the three different curves. Curves for relative differences of Bohr radii smaller than can be shown to converge to the blue curve with . The positive and negative peaks become more pronounced for decreasing ratio of Bohr radius over inter-site distance, see Fig. 2. It is observed that whereas the widths of the peaks decrease, their heights increase.

It is readily recognized that the maximum decoherence rate always occurs at the time , which is the time required for a phonon to cross the smallest donor site. A negative peak of the decoherence rate, i.e. a maximal recoherence rate, is observed always around the same time . This corresponds to the time needed for a phonon to travel between donor sites. Thus, the timescales can be interpreted as the time needed for the phonon to either travel within an impurity site or between the sites. In the first case, when the phonon (as quasiparticle) is traveling within one of the sites it destroys the inter-site coherence and does this at maximum rate at the moment it has propagated a distance equal to the smallest Bohr radius of the impurity states. In the latter case the negative minimum, i.e. the recoherence peak, appears when the phonon has traveled the distance between the impurities, thereby correlating their states.

The first time scale of the positive peak is produced by the first term in Eq. (77), whereas the second timescale of the negative peak comes from the third term. The second term, on the other hand, compensates the third one for times and provides for . It can be easily seen that when expanding the decoherence rate (76) with respect to the relative difference of (relative) Bohr radii we obtain

 γ(t)=γ0(t)+O(σ2).

Here is the decoherence rate for the limiting case of equal Bohr radii, , which is obtained as

 γ0(t) = ΓTη2{⎡⎣23(t/τdη)3+(t/τdη)2+12(t/τdη)⎤⎦e−2(t/τd)/η +η⎡⎣16(|1+t/τd|η)3+12(|1+t/τd|η)2+58(|1+t/τd|η)+516