Electron Spin Relaxation in Graphene Nanoribbon Quantum Dots
EElectron Spin Relaxation in Graphene Nanoribbon Quantum Dots
Matthias Droth and Guido Burkard
Department of Physics, University of Konstanz, 78457 Konstanz, Germany
Graphene is promising as a host material for electron spin qubits because of its predicted potentialfor long coherence times. In armchair graphene nanoribbons (aGNRs) a small band gap is opened,allowing for electrically gated quantum dots, and furthermore the valley degeneracy is lifted. Thespin lifetime T is limited by spin relaxation, where the Zeeman energy is absorbed by latticevibrations, mediated by spin-orbit and electron-phonon coupling. We have calculated T by treatingall couplings analytically and find that T can be in the range of seconds for several reasons: (i) lowphonon density of states away from Van Hove singularities; (ii) destructive interference between tworelaxation mechanisms; (iii) Van Vleck cancellation at low magnetic fields; (iv) vanishing couplingto out-of-plane modes in lowest order due to the electronic structure of aGNRs. Owing to thevanishing nuclear spin of C, T may be a good measure for overall coherence. These results andrecent advances in the controlled production of graphene nanoribbons make this system interestingfor spintronics applications. PACS numbers: 87.75.-d, 76.60.Es, 63.22.Rc
I. INTRODUCTION
Graphene has attracted intense scientific interest for itsmechanical, electronic, and other properties.
Withinthe plane of its two-dimensional lattice it is extremelyrigid while out-of-plane deformations are relatively softdue to the lack of a linear restoring force.
The ab-sence of a band gap leads to a quasi relativistic behaviorof the electrons that can be described by a Dirac-likeHamiltonian.
However, for typical semiconductor ap-plications such as transistors or spintronics devices, it isfavorable to work with a band gap.
Due to Klein’sparadox, a band gap is necessary to confine charge car-riers electrostatically in graphene.
There are differentsituations that lead to a band gap in graphene and someof them have already been studied in view of spintronicsapplications.
Armchair graphene nanoribbons (aGNRs) can exhibita band gap and in addition allow for coupling of qubitsin non-adjacent quantum dots (QDs).
Such a non-local coupling of qubits is ideal for fault-tolerant quan-tum computing and thus for scalability.
Over thepast years, there has been substantial progress towardsthe goal of controlling the GNR edge termination withinthe production process and the controlled production ofaGNRs might become feasible in the near future.
Spintronics applications like the Loss-DiVincenzoquantum computer require spin coherence times muchlonger than typical operation times.
When the qubitis represented by the real electron spin, carbon materialsare considered promising due to the small atomic spin-orbit coupling and weak interaction with nuclear spins incarbon.
While the curvature significantly enhancesintrinsic spin-orbit coupling and hence spin relaxationin carbon nanotubes, this effect should not occur in flatgraphene.
The natural abundance of C, the onlystable carbon isotope with a finite nuclear spin I = 1 / µ B (cid:29) µ nuc . We ex-pect that T is dominated by T and that the spin re-laxation time is a good measure for overall coherence, T ≈ T / T for electrons that are confined in an aGNR QD. The finitewidth of the quasi one-dimensional aGNR leads to con-finement in the transverse ( x -) direction. As we will dis-cuss in Sec. III, aGNRs of appropriate width have a bandgap. This allows us to avoid Klein’s paradox and confineelectrons in the longitudinal ( y -) direction by means of anelectrostatic potential V ( y ). In a perpendicular magneticfield B e z , the two possible spin states of an electron in-side the QD are split by the Zeeman energy gµ B B = ¯ hω ,where g = 2 is the electron g factor. Figure 1 shows a B1 D B2e-
FIG. 1: (Color online) Sketch of the system and definitionof the coordinate frame. The GNR has armchair termina-tions in the x direction. The width of the sketched aGNRis characterized by m = 3 and µ = −
1, which leads to aband gap that allows for electrostatic confinement in the y direction. The potential V ( y ) defines the two barrier regionsB1, B2 (shaded), and the dot region D, that lies symmetri-cally between the barrier regions. The interatomic distancein graphene is a = 1 .
42 ˚A. a r X i v : . [ c ond - m a t . m e s - h a ll ] J un sketch of the system.Due to energy conservation, the Zeeman energy mustbe transferred to the lattice upon spin relaxation. Fortypical laboratory magnetic fields B < ∼
20 T, the Zeemanenergy corresponds to low-energy acoustic phonons at thecenter of the Brillouin zone. We consider two cases sep-arately: (i) free and (ii) fixed boundaries. The electron-phonon coupling H EPC comprises the deformation po-tential as well as the bond length change and couplesin-plane vibrational modes to the electronic state. Byincluding the spin-orbit interaction H SOI , the spin thusbecomes connected to the vibrational state of the system.The coupling to the out-of-plane modes is considered, aswell. Yet such a coupling either vanishes identically dueto the electronic structure in aGNRs or appears only inhigher order.This paper is organized as follows. In Sec. II, wepresent our model and in Sec. III, we recapitulate thebound states of aGNR QDs and explain the extended,quasi continuous states. Acoustic GNR phonons areshortly reviewed in Sec. IV. The effective spin-phononcoupling mechanisms that lead to T − via Fermi’s goldenrule are clarified in Sec. V. In Sec. VI, we comment onthe actual evaluation of T . The results are presented inSec. VII and discussed in Sec. VIII. II. MODEL
The Hamiltonian of the system is H = H elec + H phon + H SOI + H EPC , (1)where H elec and H phon describe the unperturbed elec-tronic system and the unperturbed vibrational system,respectively. The spin-orbit interaction H SOI leads to anadmixture of opposite spin states such that the electronphonon coupling H EPC can induce a spin flip. Denot-ing the Fermi velocity by v F and the pseudospin by σ ,the unperturbed electronic part of the system obeys theHamiltonian H elec = − i ¯ hv F (cid:18) σ x ∂ x + σ y ∂ y − σ x ∂ x + σ y ∂ y (cid:19) + V ( y ) (2)with eigenstates | k (cid:105) . The pure vibrational modes aredescribed by H phon = (cid:88) α,q ¯ hω α,q (cid:18) n α,q + 12 (cid:19) , (3)where the summation runs over all phonon branches α and wave numbers q . The angular frequency ω α,q of avibrational mode is implicitly determined by α and q and n α,q is the occupation number operator. The eigenstatesare the occupation number states | n α,q (cid:105) .Since H EPC does not couple to the spin, the spin-orbitinteraction H SOI needs to be included in order to obtaina spin relaxing mechanism via admixture of electronic states. For this admixture, we consider both boundstates confined inside the dot and extended, quasi con-tinuous states energetically above the confinement poten-tial.As will be discussed in more detail, H SOI perturbsthe electron-spin product states | k (cid:105)| s (cid:105) = | k s (cid:105) (0) , where s = ↑ , ↓ . We denote the first order perturbed states by | k s (cid:105) . Finally, the electron-phonon coupling leads to fi-nite matrix elements (cid:104) k ↓|H EPC | k ↑(cid:105) . This allows us touse Fermi’s golden rule in order to calculate the spin re-laxation rate T − = 2 π ¯ h (cid:88) α,q |(cid:104) k ↓ , n α,q + 1 |H EPC | k ↑ , n α,q (cid:105)| × ρ states (¯ hω α,q ) , (4)where ρ states (¯ hω α,q ) is the phonon density of states atthe respective energy. The result is a function of threeparameters: (i) length-to-width ratio (aspect ratio) L/W of the QD, (ii) potential depth ∆ V of the QD, and (iii)perpendicular magnetic field B. We find that T can beas large as several seconds if ρ states is small and the twomechanisms in H EPC interfere destructively.
III. ELECTRONIC STATES
Due to the aGNR edges where the wave function van-ishes on both sublattices, electronic states in an aGNRhave transverse wave numbers q n = π ( n − µ/ /W , (5)where n = 0 , ± , ± , . . . and W = (3 m + µ ) √ a isthe ribbon width. The width depends on m ∈ N and µ ∈ {− , , +1 } . The interatomic distance is a = 1 .
42 ˚A.Due to Eq. (5) and E = ± ¯ hv F (cid:112) q n + k , where v F is theFermi velocity and k is the longitudinal electronic wavenumber, there is a band gap E gap = 2¯ hv F | q | . Since E gap = 0 for µ = 0, we assume µ = ± µ is determined by the number of atoms acrossthe GNR, Fig. 1. Spinors with different transverse quan-tum number n are orthogonal such that we shall focuson the lowest transverse wave number with | q | = π/ W .The resulting gap E gap = 2¯ hv F π/ W allows us to avoidKlein’s paradox and confine charge carriers electrostati-cally in a finite square potential V ( y ) = (cid:26) y ∈ D (dot region),∆ V : y ∈ B1 ∪ B2 (barrier regions). (6)The barrier region B1 extends from the left end of theaGNR to y = 0 and the barrier region B2 extends from y = L to the right end. The dot region D lies symmet-rically between the barrier regions. The resulting poten-tial landscape is shown in Fig. 2(a) together with a boundstate, which will be discussed in the next subsection. Thelength of the QD is denoted by L and assumed to bemuch smaller than the overall ribbon length, L (cid:28) L GNR . (a) DB1 B2 (a)(b)
FIG. 2: (Color online) Electron states in an aGNR QD. (a)Sketch of a bound state and (b) QD bound-state energy spec-trum given by roots of Eq. (12). (a) Due to armchair bound-aries, the minimum transverse wave number is q = ± π/ W for µ = ∓
1. As a consequence, the conduction band is sep-arated from the valence band by a gap of E gap = 2¯ hv F | q | .All energies shall be measured with respect to the middle ofthis band gap inside the QD region. In the barrier regions,both bands are shifted by the barrier height ∆ V . The re-sulting QD hosts at least one bound state. All bound stateshave the form given by Eq. (11) and decay exponentially for y → ±∞ . The arrows underneath the greek letters indicatethe directed character of the according part of the wave func-tion. The plotted probability density | ψ ( y ) | belongs to thelowest bound state for L/W = 5 and ∆ V = 1 . hv F q . (b)Bound states exist for roots of Eq. (12) and can be plotted ina ∆ V - E -plot. There is at least one bound state for all valuesof ∆ V . As ∆ V is increased, more bound states fit into theenergy gap until the lowest state can leave the QD via valencestates in the barrier regions. Notably, Eq. (12) has exactlyone root for every value of E ≥ ¯ hv F q . We enumerate thebound states by j = 0 , , , . . . . The circled position on the j = 0 line marks the state plotted in (a). For the shown plot,the aspect ratio is L/W = 5.
For concreteness, we assume an overall GNR length of L GNR = 50 W . The finite square potential needs tobe considered in the electronic dispersion relation, whichbecomes E = V ( y ) ± ¯ hv F (cid:113) q + k . (7) Provided that the barrier height ∆ V does not exceed acritical value 2¯ hv F | q | + ∆ V , we can easily order boundstates and extended states by their energies. The criticalvalue and ∆ V will be explained in the next subsection- for now, we only assume that ∆ V does not exceed it.Then, a state with energy E ∈ [¯ hv F | q | , ¯ hv F | q | + ∆ V ]is bound since its longitudinal wave number k is real inthe dot region and complex in the barrier regions, thusleading to an evanescent behavior. For E > ¯ hv F | q | +∆ V ,the longitudinal wave number is real in all regions. Thisleads to extended waves. Both bound and extended statescontribute to the admixture mechanism and thus shall bediscussed in more detail. A. Bound states
To describe bound states in aGNRs, one can assumean infinite ribbon. On one hand, L GNR will always befinite in reality. On the other hand, bound states aremainly localized in the dot region 0 ≤ y ≤ L and decayexponentially in the barrier regions, as shown in Fig. 2(a).As mentioned above, we assume L (cid:28) L GNR , such thatthe overall ribbon still appears approximately infinite forbound states. This allows us to follow the descriptionwith L GNR → ∞ for bound states. Accordingly, we denote the four-component envelopewave function by ψ = ( ψ ( K ) A , ψ ( K ) B , − ψ ( K (cid:48) ) A , − ψ ( K (cid:48) ) B ) (8)and assume plane waves along the ribbon, ψ ( ± ) n,k ( x, y ) = χ ( ± ) n,k ( x ) e ± iky , where χ (+) n,k = a (+) n (1 , z n,k , , e iq n x + b (+) n ( − z n,k , , , e − iq n x + c (+) n (0 , , − z n,k , e iq n x + d (+) n (0 , , , z n,k ) e − iq n x (9)and χ ( − ) n,k = a ( − ) n ( z n,k , , , e iq n x + b ( − ) n (1 , − z n,k , , e − iq n x + c ( − ) n (0 , , , − z n,k ) e iq n x + d ( − ) n (0 , , z n,k , e − iq n x . (10)With z n,k = ± ( q n + ik ) / (cid:112) q n + k , and longitudinalwave numbers k D = (cid:112) ( E/ ¯ hv F ) − q n (dot region), κ B = k B /i = (cid:112) q n − (( E − e ∆ V ) / ¯ hv F ) (barrier re-gions), bound states have the form ψ = α n χ ( − ) n,κ B e κ B y : y ∈ B1, β n χ (+) n,k D e ik D y + γ n χ ( − ) n,k D e − ik D y : y ∈ D, δ n χ (+) n,κ B e − κ B ( y − L ) : y ∈ B2. (11)The matching conditions at the interfaces B1/D andD/B2 (that is, at y = 0 , L ) are discussed in Ref. [12]and can be met for roots of the transcendental equationtan( k D L ) = k D κ B ± (cid:112) q n − κ (cid:112) q n + k − q n . (12) B1B1 D B2
FIG. 3: (Color online) Sketch of an extended state. Thepotential landscape, the aspect ratio
L/W , and the barrierheight ∆ V are the same as in Fig. 2(a) for bound states. Theplotted probability density belongs to an extended state thatis incident from the left as described by Eq. (14) and for which k EB = 20 π/L GNR . The arrows underneath the greek lettersindicate the direction of propagation of the according part ofthe wave function.
For | q | = π/ W and L/W = 5, Fig. 2(b) shows theseroots as a function of the barrier height ∆ V . There isa finite number of longitudinal excitations for any given∆ V . The different bound states can be enumerated by j = 0 , , , . . . and have distinct coloring in our figure.The j -th bound state has j nodes inside the dot re-gion. For a given excitation, ∆ V can be increased un-til the valence band reaches the energy of the loweststate, which can now leave the QD via valence statesin the barrier regions. Note that this occurs exactlywhen the argument on the left-hand side of Eq. (12)equals a multiple of π . For the ground state, this means k D ∈ [0 , π/L ] such that the maximum ground state en-ergy is E , max = ¯ hv F (cid:112) q + ( π/L ) . States of higher en-ergy belong to the j -th longitudinal excitation ( j > V j = ¯ hv F ( jπ/L ). For ∆ V < ∆ V ,the ground state is the only bound state. This will beimportant for the evaluation of T , see Secs. VI and VII.The critical value for ∆ V mentioned before is ∆ V =2¯ hv F | q | +∆ V . If the barrier height surpasses this value,the lowest state inside the QD can leave it via valencestates in the barrier region. That is the state becomesextended thus affecting the ordering of bound and ex-tended states. Throughout this paper we assume that∆ V does not exceed this threshold such that the groundstate belongs to j = 0. B. Extended states
We assume L GNR = 50 W for the overall length ofthe GNR such that possible wave numbers are k EB =0 , ± π/L GNR , . . . , ± π/a with lattice constant a . Since energy is conserved, the wave number becomes k ED = (cid:115)(cid:18)(cid:113) q n + k + ∆ V / ¯ hv F (cid:19) − q n (13)in the dot region. Depending on the sign of k EB , thestate is incident from y <
0, leading to ψ = (cid:15) n χ (+) n,k EB e ik EB y + α n χ ( − ) n,k EB e − ik EB y : y ∈ B1, β n χ (+) n,k ED e ik ED y + γ n χ ( − ) n,k ED e − ik ED y : y ∈ D, δ n χ (+) n,k EB e ik EB ( y − L ) : y ∈ B2,(14)see Fig. 3, or it is incident from y > L , which is describedby ψ = α n χ ( − ) n,k EB e − ik EB y : y ∈ B1, β n χ (+) n,k ED e ik ED y + γ n χ ( − ) n,k ED e − ik ED y : y ∈ D, δ n χ (+) n,k EB e ik EB ( y − L ) + (cid:15) n χ ( − ) n,k EB e − ik EB ( y − L ) : y ∈ B2.(15)The matching conditions at y = 0 , L can always be met.In contrast to bound states, extended states are propa-gating waves in the barrier regions. IV. ACOUSTIC GNR PHONONS
The phonon energies we are interested in need tomatch the Zeeman splitting, ¯ hω = gµ B B , where ω isthe phonon frequency, g the electron g factor, and µ B denotes Bohr’s magneton. For typical laboratory mag-netic fields B < ∼
20 T, this implies low-energy acousticphonons at the center of the Brillouin zone, which can bemodeled by continuum mechanics.
In this model, de-formations are described by the displacement field u ( r ).While the components u xz and u yz of the strain tensor u ik = ( ∂ i u k + ∂ k u i ) / u zz must vanish, as well. With u iz ≡
0, the elas-tic Lagrangian density of monolayer graphene is givenby L = T − V = ρ u − κ (cid:52) u z ) − B + µ u ii + µu ik , (16)where (cid:52) = ∂ x + ∂ y , the sum convention with u ii = u xx + u yy + u zz and u ik = u xx + u xy + · · · has been used, ρ isthe mass density, and κ is the bending rigidity. The bulk( B ) and shear ( µ ) moduli can be expressed by Poisson’sratio σ and Young’s modulus E . The numerical values ofelastic and other constants we use are listed in Table I.Equation (16) shows that in-plane vibrations u (cid:107) decouplefrom out-of-plane vibrations u ⊥ . By assuming u i ( x, y ) = f i ( x ) exp[ i ( qy − ωt )] for a single mode and imposing freeor fixed conditions as discussed in Ref. [30], we obtainthe according phonon dispersions (Fig. 4) and the explicit (a) (b) free boundaries n m fixed boundaries FIG. 4: (Color online) Phonon dispersion for (a) free and (b)fixed boundaries. The dimensionless frequency ¯ ω xy is con-nected to the physical frequency by ¯ ω xy = ω (cid:112) ρ/ E W , wherethe radicand contains elastic constants listed in Table I. Werestrict our interest to the frequency range ¯ ω xy ∈ [0 ,
5] sincefor W = 30 nm, the upper bound already relates to a mag-netic field of 20 T. The scale on the right-hand side showsthe magnetic field for W = 30 nm. Due to parity with re-spect to x , only the labeled branches ( α i for free and ˜ α forfixed boundaries) assist in spin relaxation. (a) The phononspectrum is gapless for free boundaries. The branch α hasa minimum and hence a diverging density of states for finite¯ q . Its constituent parts α , and α , shall be treated sepa-rately. (b) Fixed edges lead to gapped phonon spectrum. For W = 30 nm, this gap corresponds to 8 .
25 T. In our range ofinterest, the branch labeled ˜ α provides the only channel forspin relaxation. displacement fields. The latter can be quantized and thentake the form u || = (cid:88) α,q r α,q ( f xα,q e x + f yα,q e y ) e iqy , (17) u ⊥ = (cid:88) α,q r α,q f zα,q e z e iqy , (18)where q is the phonon wave number, α labels the phononbranch, and r α,q = (cid:113) ¯ h/ (2 ρLW ω α,q )( b α,q + b † α, − q ) (19)is the normal coordinate. The operator b α,q ( b † α,q ) anni-hilates (creates) a phonon on branch α with wave num-ber q . As discussed in the following section, we can ne-glect coupling to out-of-plane modes and focus on in-plane modes. The dimensionless frequency of in-planemodes ¯ ω xy is related to the physical frequency by ¯ ω xy = ω (cid:112) ρ/ E W . In the continuum model, all branches ex-tend to infinity but we are only interested in the range¯ ω xy ∈ [0 ,
5] since for a typical GNR width of W = 30 nm,¯ ω xy = 5 relates to a magnetic field of 20 T. The di-mensionless wave number ¯ q is obtained from the physicalwave number via ¯ q = qW . For symmetry reasons ex-plained in Sec. VI, not all branches contribute to spinrelaxation but only those which are explicitly labeled inFig. 4.In the case of free boundaries, the branches α , α ,and α are relevant. While α extends throughout theconsidered interval, α only exists above ¯ ω xy = 4 .
16, and α needs further discussion. Its minimum is 3 .
05 andoccurs at a finite value ¯ q where the density of stateshas a Van Hove singularity. We consider two parts: for¯ q < ¯ q , we label the branch α , and its label for ¯ q > ¯ q is α , . The range of α , is ¯ ω xy ∈ [3 . , .
18] and α , extends from its minimum to infinity.For fixed boundaries, we only need to consider thebranch ˜ α . It extends from ¯ ω xy = 2 .
06 to infinity. Weemphasize that for a typical GNR width of 30 nm, single-phonon processes do not occur up to 8 .
25 T as there areno phonons below ¯ ω xy = 2 .
06 for fixed boundaries.
V. COUPLING MECHANISMS
Phonons do not couple to the electron spin directly.The relevant mechanism usually involves the spin-orbitinteraction.
In graphene, the spin-orbit interactionis given by H SOI = H I + H R = λ I τ z σ z s z + λ R ( τ z σ x s y − σ y s x ) (20)and we will consider it in order to obtain an indirect spin-phonon coupling. The intrinsic (Dresselhaus) term H I has coupling strength λ I and the Rashba (or extrin-sic) term H R couples with strength λ R . The valley isdenoted by τ z , pseudospin by σ , and real spin by s . Inthe following, we show how the real electron spin canbe connected to the vibrational state of the system bytaking the spin-orbit interaction into account. A. Coupling to in-plane modes
In first order perturbation theory, H R corrects theelectron-spin product states | k (cid:105)| ↑(cid:105) = | k ↑(cid:105) (0) to | k ↑(cid:105) = | k ↑(cid:105) (0) + (cid:88) k (cid:48) (cid:54) = k | k (cid:48) ↓(cid:105) (0) (0) (cid:104) k (cid:48) ↓ |H R | k ↑(cid:105) (0) E k − E (cid:48) k + gµ B B , (21)and | k ↓(cid:105) accordingly. We emphasize that the summationindex k (cid:48) runs over both bound and extended states. Thepotential depth and the aspect ratio determine how manybound states exist, Fig. 2(b). For extended states, weconsider all wave numbers inside the first Brillouin zone, k EB = 0 , ± π/L GNR , . . . , ± π/a .The second term in Eq. (21) admixes states with op-posite spin such that the electron-phonon coupling H EPC can induce a spin flip (cid:104) k ↓ |H EPC | k ↑(cid:105) = (cid:88) k (cid:48) (cid:54) = k (cid:34) ( H EPC ) kk (cid:48) ( H R ) ↓↑ k (cid:48) k E k − E k (cid:48) + gµ B B + ( H EPC ) k (cid:48) k ( H R ) ↓↑ kk (cid:48) E k − E k (cid:48) − gµ B B (cid:35) , (22)where we denote the numerator in Eq. (21) as ( H R ) ↓↑ k (cid:48) k and the spin-conserving transitions of H EPC accordingly.We find that for a given k (cid:48) , the two terms in Eq. (22) ex-actly cancel each other at B = 0. This effect is known asVan Vleck cancellation and is expected for time-reversal-symmetric systems. Moreover, ( H R ) ↓↑ kk (cid:48) vanishes if both k and k (cid:48) represent bound states and the longitudinal ex-citation indices j k , j k (cid:48) (see Fig. 2) are both even or bothodd. In the electron phonon coupling Hamiltonian H EPC ,we consider the deformation potential H DP as well asbond length change H BLC : H EPC = H DP + H BLC , (23) H DP = g ∇ · u (cid:107) , H BLC = g A x − iA y A x + iA y A x + iA y A x − iA y , where g , are coupling constants, ( A x , A y ) = ( u xx − u yy , − u xy ), and the basis of Eq. (8) has been used. B. Vanishing out-of-plane deflection coupling
Low-energy acoustic phonons at the center of the Bril-louin zone have a wavelength much larger than the latticeconstant and produce a local tilt of the GNR. In the lo-cal ribbon frame Σ (cid:48) where n = e (cid:48) z is the vector normalto the ribbon plane the local spin matrix is describedby s (cid:48) z = s z − ∂ x u z s x − ∂ y u z s y . As a consequence, theintrinsic spin-orbit interaction H I = λ I τ z σ z ( s z − ∂ x u z s x − ∂ y u z s y ) (24)becomes dependent on out-of-plane phonons such thatthese could flip the spin. This is known as deflectioncoupling. However, there is a proportionality to τ z in Eq. (24). Since the electronic states we use have theproperty that the wave function has equal weight on eachvalley, the contributions from K and K (cid:48) add up to zeroand the deflection coupling between spin and out-of-planemodes vanishes.If the spin-orbit admixed states of Eq. (21) are used,there is a finite overlap only between both admixed parts such that the resulting mechanism is proportional to λ I λ R and hence negligible.Compared to in-plane phonons, both deformation po-tential and bond length change appear only in higherorder such that we neglect these mechanisms for out-of-plane phonons. VI. EVALUATION OF T Using Eq. (4), we calculate the spin relaxation rate forthe electron in the lowest bound state (ground state) ofthe QD. For concreteness, we assume µ = − x = 0, H BLC and H DP are even or odd in x , depending on what phonon branchthey belong to. Due to their similar form, the mecha-nisms are either both even or both odd for a given branch.The x dependencies of the electronic states in the ma-trix element ( H EPC ) k (cid:48) k cancel out, e i ( q − q ) x = 1, suchthat the x integral vanishes if Eq. (23) is odd in x . Thebranches α , α , and α in Fig. 4(a) and ˜ α in Fig. 4(b)have couplings H EPC that are even in x and hence canrelax the spin.For a given relaxation channel ( α, q ), both mecha-nisms H DP and H BLC are combined in Eq. (4) coherently.Moreover, the couplings via bound states and extendedstates in Eq. (21) are added up in a coherent way. Weare interested in the relaxation of the spin in the groundstate, which corresponds to j = 0 in Fig. 2(b) and hencerestrict the barrier height to ∆ V ∈ [0 , hv F q + ∆ V ].If ∆ V exceeds this upper bound, valence states becomeavailable in the barrier regions and the lowest state in-side the QD can leave the dot region. For ∆ V < ∆ V onthe other hand, the ground state is the only bound statesuch that the perturbation in Eq. (21) comes about onlydue to extended states, which fully determine the spinrelaxation in this case.For spin relaxation, Eq. (4) is proportional to n α,q + 1and we assume n α,q = 0, i. e., sufficiently low tempera-ture, k B T (cid:28) ¯ hω = gµ B B . By k B we denote Boltzmann’sconstant. Assuming a magnetic field of B = 1 T, thismeans T (cid:28) . T > ∼
15 K, spontaneous emissioncan be neglected since n α,q (cid:29) (cid:104) n α,q ( B, T ) (cid:105) = (cid:18) e gµ B Bk B T − (cid:19) − . (25)The spin relaxation time T is a good measure for over-all coherence when pure dephasing, which comes fromcoupling to nuclear spins, is negligible. Due to the lowdensity of nuclear spins in natural carbon and the verydifferent magnetic moments µ B (cid:29) µ nuc , we expect thatflip-flop processes with nuclear spins can be neglected formagnetic fields above 10 mT. For a typical GNR widthof W = 30 nm, 10 mT correspond to ¯ ω xy = 0 . ω xy ∈ [0 . , VII. RESULTS
To calculate T , we need to use the specific valuesof the elastic constants that define the phonon spec-trum. Young’s modulus for the two-dimensional latticeof graphene is obtained by multiplying the bulk valuewith the thickness associated with graphene, E = E h ,where h = 3 . σ = 0 . [39–41] g = 30 eV [25,34,35] E = 3 . [39–41] g = 1 . [25,35] B = 12 . λ R = 40 × − eV [25,28,29] µ = 9 . v F = 8 . × m/s [8,25,29] ρ = 7 . × − kg/m TABLE I: Numerical values of the parameters we use in ourcalculation.
The spin relaxation time T depends on three parame-ters: (i) the aspect ratio L/W of the QD, (ii) the poten-tial depth ∆ V of the QD, and (iii) the applied perpen-dicular magnetic field B ∝ ¯ ω xy . Moreover, the phononspectrum and hence the spin relaxation depends on themechanical boundary conditions. We discuss free bound-ary conditions separately from fixed boundaries. A. Free boundary conditions
For symmetry reasons explained above, only thephonon branches with labels α , α (consisting of parts combined bDP bBLC eDP eBLC (a)(b)(d) (c) nm FIG. 5: (Color online) Partial rates for various relaxationchannels. For
L/W = 5 and ∆ V = 1 . hv F q , all contributionsto the four relaxation channels α (a), α , (b), α , (c), and α (d) are shown. The contributions stem from H DP withadmixture of bound (labeled “bDP”) states or extended states(“eDP”) and from H BLC with the same admixtures (“bBLC”and “eBLC”, respectively). These contributions are addedup coherently to the “combined” relaxation of the respectivechannel. α , and α , ), and α in Fig. 4(a) need to be consid-ered. The respective rates of these relaxation channelsare shown in Fig. 5 for an aspect ratio of L/W = 5and a barrier height of ∆ V = 1 . hv F q . The allocationof panels to branches is as follows: Fig. 5(a) belongs tobranch α , 5(b) to α , , 5(c) to α , , and 5(d) to α . Eachpanel shows four separate contributions to T − that comeabout from the two mechanisms in Eq. (23) and the ad-mixture of bound states or of extended states in Eq. (21)for each mechanism. The coherent sum of all four con-tributions is displayed by the gray line. The deformationpotential usually dominates over the bond length changesince its coupling constant is 20 times larger, Table I.For ∆ V = 1 . hv F q , extended states are energeticallyfar away from the ground state such that the contribu-tion from the deformation potential with admixture ofbound states dominates in Fig. 5. Oscillations in indi-vidual rates may be due to the phonon phase e iqy thatis integrated with the matrix elements ( H EPC ) k (cid:48) k androtates according to the phonon dispersion when ω ischanged. Figures 5(c) and 5(d) show that the matrixelements ( H DP ) k (cid:48) k and ( H BLC ) k (cid:48) k may interfere destruc-tively, thus decreasing T − by several orders of magni-tude, yet typically not to zero.In all these plots, the bottom scale shows ¯ ω xy and thetop scale shows the magnetic field B that corresponds to¯ ω xy , assuming a width of W = 30 nm. Note, that T − does not depend on B and W separately, but only on theproduct BW ∝ ωW ∝ ¯ ω xy .Figure 6(a) shows the full spin relaxation rate for thesituation of Fig. 5, that is, the combined rates of all relax-ation channels α , α , , α , , and α (gray lines in Fig. 5)are summed up to the full relaxation rate T − [gray linein Fig. 6(a)]. The rate with the label “bound” (“ex-tended”) is obtained in a similar fashion, but only con-tributions with admixture of bound (extended) states areconsidered, here. For ∆ V = 1 . hv F q , the admixture ofbound states dominates the admixture of extended statesby several orders of magnitude. Yet by lowering ∆ V , theinfluence of extended states can be close to [Fig. 6(b)]or even surpass the influence of the bound states. Fig-ure 7 shows two cases for an aspect ratio of L/W = 2.In Fig. 7(a), the barrier height is ∆ V = 2 . hv F q andextended states are basically irrelevant compared to therelaxation via bound states. However, Fig. 7(b) showsthat for ∆ V = 0 . hv F q , the major contribution comesfrom the extended states.Figure 8(a) shows T − as a function of parameters ∆ V and ¯ ω xy ∝ B , and for a fixed aspect ratio of L/W = 5.In contrast to ¯ ω xy , the barrier height hardly changes thequalitative picture. The orange cut at ∆ V = 1 . hv F q is repeated in Fig. 8(b) in a doubly logarithmic plotthat highlights the B dependence in the range ¯ ω xy ∈ [0 . , . α is availableand has a linear dispersion B ∝ ω ∝ q . The matrix ele-ments ( H EPC ) k (cid:48) k have one power in B due to (i) the gra-dients ∝ q in Eq. (23), (ii) dipole approximation ∝ q , and(iii) Van Vleck cancellation ∝ B , each. Because of theprefactor ∝ ω − . in Eq. (19), we find ( H EPC ) k (cid:48) k ∝ B . .As α is linear and hence ρ states ∝ B for this low-fieldregime, this explains T − ∝ B . Destructive interferenceof matrix elements ( H DP ) k (cid:48) k and ( H BLC ) k (cid:48) k can lead toa very small but nonzero relaxation rate. full bound extended (a) nm (b) FIG. 6: (Color online) The relaxation rates for different dotdepths ∆ V . By summing up the combined relaxation rates(see Fig. 5) of all channels available for a certain ¯ ω xy , thefull relaxation rate (gray line) is obtained. The lines labeled“bound” and “extended” are obtained in a similar way byconsidering only bound or extended states, respectively. At¯ ω xy = 3 . T − is discontinuous due to the advent of therelaxation channel α that has a diverging density of statesat this point, Fig. 4(a). (a) accords to parameters L/W = 5and ∆ V = 1 . hv F q as in Fig. 5. Clearly, the energeticallyfar off extended states play a negligible role for such a deepdot. In (b), the barrier height is reduced to 0 . hv F q suchthat extended states are about as important as bound states. Figure 9(a) shows a plot similar to Fig. 8(a), yet for
L/W = 2. The qualitative picture is much different fromthe aspect ratio
L/W = 5. Figures 6 - 9(a) all show dis-continuities at ¯ ω xy = 3 .
05 that stem from the branch α ,for which the density of states has a Van Hove singularityat ¯ q while the coupling H EPC remains finite, Fig. 4(a).
B. Fixed boundary conditions
Most importantly, fixed boundaries result in a gappedphonon spectrum. This means that spin relaxation in-volving only one phonon cannot occur for magnetic fieldsthat correspond to ¯ ω xy < .
06. Note, that for a typicalwidth W = 30 nm, ¯ ω xy = 2 .
06 corresponds to a mag-netic field of 8 .
25 T. However, phonon scattering maystill take place below this threshold. In contrast to ourclaim in Ref. [30], only the branch ˜ α contributes to the full bound extended (a) nm (b) FIG. 7: (Color online) This plot shows the same quantitiesas Fig. 6, yet for the aspect ratio
L/W = 2. Again, theinfluence of extended states depends on the barrier height:∆ V = 2 . hv F q in (a) and ∆ V = 0 . hv F q in (b). The ex-tended states dominate in the latter case. spin relaxation rate. Gradients, dipole approximation,and Van Vleck cancellation play the same role as for freeboundaries, yet due to the gap the frequency ω is not pro-portional to some power of q such that there is no powerlaw that connects T − and B as for free boundaries.Figure 9(b) shows an analog to Fig. 9(a), yet for fixedboundaries. For aspect ratios larger than in Fig. 9(b),oscillations occur, which can again be explained with thephonon phase e iqy that rotates according to the phonondispersion when ω changes. These oscillations arise onlyif the dot length is large enough. VIII. DISCUSSION
The spin relaxation times we find in our work rangefrom 10 − seconds to beyond the second range. Forcases where T is very long, it can be expected that othermechanisms not considered here will dominate. Our re-sults depend on the aspect ratio L/W , the barrier height∆ V , and the Zeeman splitting gµ B B ∝ ¯ ω xy but also onthe mechanic boundary conditions that lead to differentphonon dispersions. By choosing/adjusting these degreesof freedom properly, T can be in the range of seconds.We attribute such long relaxation times to several effects: (a) nm nm (b) FIG. 8: (Color online) Spin relaxation rate T − for an aGNRwith aspect ratio L/W = 5 and free edges. (a) The rate isshown as a function of barrier height ∆ V and phonon fre-quency ¯ ω xy . The orange cut corresponds to the gray linein Fig. 6(a) and is repeated in (b) with a doubly logarith-mic scale that highlights the B dependence in the interval¯ ω xy ∈ [0 . , . (i) GNRs are quasi one-dimensional systems similarto carbon nanotubes. Both the phonon and the elec-tron density of states are thus limited compared to bulkgraphene. (ii) Destructive interference between the deformationpotential and the bond length change as well as oscilla-tions due to the phonon phase e iqy that rotates accordingto the phonon dispersion when ω changes both reducethe relaxation rate T − by several orders of magnitudefor specific magnetic fields.(iii) In contrast to other graphene QD systems,the electronic states in aGNRs are invariant undertime-reversal symmetry, which leads to Van Vleckcancellation. As a result, Eq. (22) vanishes for B = 0.(iv) Deflection coupling to out-of-plane modes vanishesas the evenly distributed weights on K and K (cid:48) spinorcomponents cancel out. As a result, only the very rigid0 (a)
16 20 nm (b) FIG. 9: (Color online) The relaxation rate T − for (a) freeand (b) fixed mechanic boundaries. (a) This case is similar toFig. 8(a) yet with aspect ratio L/W = 2. (b) Fixed boundaryconditions and
L/W = 2. Due to the gapped phonon spec-trum, the rate T − vanishes below ¯ ω xy = 2 .
06 for our modeland with fixed boundaries. Moreover, T − is not discontin-uous in ¯ ω xy as the branch ˜ α never becomes flat for finite ¯ q ,see Fig. 4(b). in-plane modes need to be considered. This rigidity leadsto a generally small density of phonon states ρ states . (v) Phonons do not couple to spin directly so that spin-orbit coupling needs to be included. However, spin-orbitcoupling in graphene is weak compared to other systems(e.g. carbon nanotubes). (vi) The admixture of electronic states in Eq. (21) in-cludes bound and extended states. However, only everysecond bound state contributes; for parity in y direction, j k and j k (cid:48) may not be even or odd at the same time,Fig. 2(b). States that are energetically far apart from theground state play a small role in the sum which is usu-ally the case for extended states, depending on ∆ V . Asa consequence, the admixture of these electronic states issuppressed. (vii) Due to parity in the x direction, not all phononbranches contribute to spin relaxation but only thosewith explicit labels in Fig. 4, for which Eq. (23) is evenin x . This limits the number of relaxation channels. It is an open question how strong the avoided relaxationchannels contribute to T − if this symmetry is broken.(viii) We assume phonon vacuum in Eq. (4). A fi-nite temperature can be included by multiplying the rate T − with the expectation value of the Bose distribution (cid:104) n α,q ( B, T ) (cid:105) as explained in Sec. VI.The carbon isotope C has no nuclear spin and thenatural abundance of C, which has spin 1 /
2, is only1%. Thus, pure dephasing, which comes from coupling tonuclear spins is likely to play a minor role in graphene de-vices and T ≈ T / IX. ACKNOWLEDGEMENTS
We thank the European Science Foundation and theDeutsche Forschungsgemeinschaft (DFG) for supportwithin the EuroGRAPHENE project CONGRAN andthe DFG for funding within SFB 767 and FOR 912. P. R. Wallace, Phys. Rev. , 622 (1947). K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M.I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A.Firsov, Nature (London) , 197 (2005). J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S.Novoselov, T. J. Booth, and S. Roth, Nature (London) , 60 (2007). R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov,T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim,Science , 1308 (2008). A. Fasolino, J. H. Los, and M. I. Katsnelson, Nature Mater. , 858 (2007). D. Gazit, Phys. Rev. B , 113411 (2009). M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, NaturePhys. , 620 (2006). A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Rev. Mod. Phys. , 109(2009). K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y.Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science , 666 (2004). D. B. Farmer, H.-Y. Chiu, Y.-M. Lin, K. A. Jenkins, F.Xia, and P. Avouris, Nano Lett. , 4474 (2009). R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha,and L. M. K. Vandersypen, Rev. Mod. Phys. , 1217(2007). B. Trauzettel, D. V. Bulaev, D. Loss, and G. Burkard,Nature Phys. , 192 (2007). O. Klein, Z. Physik , 157 (1929). P. Recher and B. Trauzettel, Nanotechnology , 302001(2010). L. Brey and H. A. Fertig, Phys. Rev. B , 235411 (2006). M. Y. Han, B. ¨Ozyilmaz, Y. Zhang, and P. Kim, Phys.Rev. Lett. , 206805 (2007). M. Braun, P. R. Struck, and G. Burkard, Phys. Rev. B K. M. Svore, B. M. Terhal, and D. P. DiVincenzo, Phys.Rev. A , 022317 (2005). L. Jiao, L. Zhang, X. Wang, G. Diankov, and H. Dai, Na-ture (London) , 877 (2009). X. Wang, Y. Ouyang, L. Jiao, H. Wang, L. Xie, J. Wu, J.Guo, and H. Dai, Nature Nanotechnol. , 563 (2011). X. Zhang, O. V. Yazyev, J. Feng, L. Xie, C. Tao, Y.-C.Chen, L. Jiao, Z. Pedramrazi, A. Zettl, S. G. Louie, H.Dai, and M. F. Crommie, ACS Nano , 198 (2013). D. Loss and D. P. DiVincenzo, Phys. Rev. A , 120(1998). D. P. DiVincenzo and D. Loss, Superlattice Microst. ,419 (1998). D. V. Bulaev, B. Trauzettel, and D. Loss, Phys. Rev. B , 235301 (2008). P. R. Struck and G. Burkard, Phys. Rev. B , 125401(2010). F. Kuemmeth, S. Ilani, D. C. Ralph, and P. L. McEuen,Nature (London) , 448 (2008). C. L. Kane and E. J. Mele, Phys. Rev. Lett. , 226801(2005). H. Min, J. E. Hill, N. A. Sinitsyn, B. R. Sahu, L. Kleinman,and A. H. MacDonald, Phys. Rev. B M. Gmitra, S. Konschuh, C. Ertler, C. Ambrosch-Draxl,and J. Fabian, Phys. Rev. B , 235431 (2009). M. Droth and G. Burkard, Phys. Rev. B , 155404 (2011). A. V. Khaetskii and Y. V. Nazarov, Phys. Rev. B ,125316 (2001). M. S. Rudner and E. I. Rashba, Phys. Rev. B , 125426(2010). L. D. Landau and E. M. Lifshitz,
Theory of Elasticity (Pergamon Press, New York, 1986). H. Suzuura and T. Ando, Phys. Rev. B , 235412 (2002). E. Mariani and F. von Oppen, Phys. Rev. B , 155411(2009). J. H. Van Vleck, Phys. Rev. , 426 (1940). See Sec. III in this paper and Eqs. (34), (35) in the sup-plementary information of Ref. [12]. Both mechanisms contain derivatives ∂ x,y . Due to the plateequation used in continuum mechanics (see Eq. (3) inRef. [30]), f x ( x ) is odd when f y ( x ) is even and vice versa.While the derivative ∂ x applied to an even function re-turns an odd function and vice versa, the derivative ∂ y corresponds merely to a multiplication with iq . C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science ,385 (2008). R. Faccio, P. A. Denis, H. Pardo, C. Goyenola, and `A. W.Mombr´u, J. Phys.: Condens. Matter , 285304 (2009). K. N. Kudin, G. E. Scuseria, and B. I. Yakobson, Phys.Rev. B , 235406 (2001). This value follows directly from the atomic weight ofnatural carbon, 12 .
01 u, and the interatomic distance ingraphene, 1 ..