A Supplemental Material

# Valley-based Cooper Pair Splitting via Topologically Confined Channels in Bilayer Graphene

## Abstract

Bilayer graphene hosts valley-chiral one dimensional modes at domain walls between regions of different interlayer potential or stacking order. When such a channel is brought into proximity to a superconductor, the two electrons of a Cooper pair which tunnel into it move in opposite directions because they belong to different valleys related by the time-reversal symmetry. This is a kinetic variant of Cooper pair splitting, which requires neither Coulomb repulsion nor energy filtering but is enforced by the robustness of the valley isospin in the absence of atomic-scale defects. We derive an effective model for the guided modes in proximity to an -wave superconductor, calculate the conductance carried by split and spin-entangled electron pairs, and interpret it as a result of local Andreev reflection processes, whereas crossed Andreev reflection is absent.

###### pacs:
72.80.Vp,74.45.+c,03.65.Ud

Creating mobile nonlocal spin-entangled electrons in a transport experiment with the help of superconductor–normal junctions has attracted a lot of attention in theory Recher et al. (2001); Lesovik et al. (2001); Recher and Loss (2002); Bena et al. (2002); Recher and Loss (2003); Levy Yeyati et al. (2007); Cayssol (2008); Sato et al. (2010) and experiment Hofstetter et al. (2009); Herrmann et al. (2010); Schindele et al. (2012); Das et al. (2012); Tan et al. (2015); Deacon et al. (2015) because the spin degree of freedom of the electron could serve as a solid-state qubit Loss and DiVincenzo (1998). In the existing experiments, the envisaged process where a Cooper pair is split over two normal leads is crossed Andreev reflection (CAR) Torrès and Martin (1999); Falci et al. (2001), which is enhanced by the repulsive electron-electron interaction on two quantum dots weakly coupled to the superconductor Recher et al. (2001) or by energy filtering Lesovik et al. (2001); Sadovskyy et al. (). The basic mechanism of these entanglers is not very sensitive to the specific material used, i.e., the underlying band structure. It has been shown that characteristic features of new materials exhibiting Dirac-cones like graphene or topological insulators can be useful for splitting Cooper pairs Cayssol (2008); Chen et al. (2011); Reinthaler et al. (2013); Nilsson et al. (2008); Wang et al. (2015). In these proposals, the efficiency of the splitting process, in the absence of interactions, relies on non-protected resonance conditions or the split Cooper pair is not spin-entangled due to spin-helicity or spin-polarization of the leads. Helical edge states of the quantum spin Hall regime have, however, been proposed to detect spin entanglement Sato et al. (2010); Chen et al. (2012).

Here, we propose to exploit the valley degree of freedom in bilayer graphene (BG), where valley-chiral, spin-degenerate one-dimensional (1D) channels are formed at domain walls. Such domain walls can be engineered by switching the sign of an interlayer voltage or by reversing the stacking order Martin et al. (2008); Zhang et al. (2013). If brought into proximity to a superconductor, the pairs emitted into the channel are split, i.e., two electrons propagate in different directions but remain spin-entangled since, as required by time-reversal symmetry, the two electrons forming the Cooper pair in the superconductor are from different valleys Beenakker (2008). As long as the valley degree of freedom is robust, the splitting efficiency is unity, independent of resonance conditions. The device extends the upcoming “valleytronics” in graphene Rycerz et al. (2007) to nonlocal Einstein-Podolsky-Rosen pairs. A 1D channel defined by opposite stacking order has recently been created experimentally in BG Ju et al. (2015) with mean free paths over several 100 nm, demonstrating weak intervalley scattering. In this scenario, the normal reflection of an incoming hole (or electron) and CAR are absent. In the limit of a weak proximity effect, where the normal transmission through the proximity region of the channel is almost perfect, the spin-entangled pair emission with electrons moving to opposite normal terminals (Fig. 1) is equivalent to local Andreev reflection (LAR) processes, opposite to the normal reflection-dominated case, where CAR produces entangled pairs.

We analyze the setup of Fig. 1 in two steps: first, we investigate the influence of the superconductor on the 1D channel by solving a Bogoliubov-de Gennes (BdG) equation, and derive an effective 1D model to describe the proximity effect in the channel. Second, we calculate the subgap conductance when applying a bias voltage between the superconductor and the channel using a rate equation approach. We interpret the subgap transport in a scattering matrix picture and show that the outgoing scattered state is a two-particle spin-entangled state on top of a filled normal-state Fermi sea with a chemical potential lowered by the bias voltage. To leading order its weight is given by the LAR amplitude.

Model.—We consider a BG sheet with Bernal stacking in the presence of an interlayer voltage McCann and Fal’ko (2006). We model the superconductor region as BG in which the bands are shifted by a scalar potential due to doping and which has an induced -wave pairing amplitude . We employ the low-energy approximation for BG, valid at energies and (inter)layer voltages smaller than the interlayer hopping  eV. Without the superconductor, the valley index distinguishing the two -points and the electron spin are good quantum numbers and we write the Bogoliubov-de Gennes equation as

 Hχv,sBdG= [αℏ2(2χv∂x∂yσy+(∂2x−∂2y)σx) +U(r)+V(r)σz]τz+Δ(r)τx, (1)

where . The Pauli matrices act in the pseudospin space and in electron-hole space and we set the Fermi energy . The 4-component spinor is where we have introduced the electron and hole components () on the two sublattices. Excitations with energy are then expanded as with the vector of field-operators .

In the absence of the superconductor () and assuming the modes to propagate along the -direction along a domain wall at , i.e., with , the topologically confined modes can be found analytically Martin et al. (2008). The electron and hole sectors in Eq. (1) decouple. In the electron sector, the solutions in each half space have the form where with

 u0χv,s=(−ε−Vα2(px+iχvpy)2). (2)

For any fixed energy and momentum there are four allowed values which become complex when , i.e. there are no propagating modes in the bulk at energies below . Matching the wavefunctions decaying away from the domain wall and their derivatives, one obtains the two electronic subgap solutions in each valley, , with the dispersion relation

 ε0,1/2χv,s,py=±√2V0−αp2y2−χvpy2√2√2αV0+α2p2y, (3)

with velocities opposite in the two valleys Martin et al. (2008). The solutions for the hole sector where at energy are obtained from Eqs. (2) and (3) by setting .

The relevant momenta are close to the points: taking , we obtain from Eq. (3) the momentum scale , on which the K-points are located at for . The guided modes decay into the bulk on a length scale of , which then is on the order of several  nm. This sets the scale of the separation between the guided mode and a superconductor required to obtain a proximity effect.

Perturbation theory for superconducting pairing.—Assuming a superconductor/normal interface with translational invariance along the -direction, there are three distinct areas: in the superconductor area, , the pairing amplitude is finite and is negative. The area is in the normal state as before, , but the interlayer voltage is finite, . This region is a tunnel barrier between the superconductor and the domain wall at the interface to the third region , where and . In this situation guided modes exist at because states above can propagate in the normal regions and states above can propagate in the superconductor. Because of the tunnel barrier the guided modes are only weakly affected by the superconductor and we can apply standard quasidegenerate perturbation theory Winkler (2003); sup (), for which the unperturbed Hamiltonian is obtained from by setting everywhere and so the perturbation , which adds the missing parts, is finite only at . As a result of the perturbation the electron and hole states of the channel acquire a finite overlap , where label the subgap bands. This allows for particle number non-conserving processes, i.e., Cooper pair transport. To first order only the electron and hole states belonging to the same subgap band mix, , . This agrees with the result one expects when introducing superconductivity phenomenologically by constructing the BdG equation directly from the guided modes with a uniform pairing . The second order corrections, which take into account the modification of the wavefunctions due to the superconductor, however, reveal that the situation is different in the geometry we consider. The electron hole overlap differs in both bands, , and band mixing is finite, (Fig. 2) ?. This is confirmed by the full dispersion relation of we obtain by matching the 4-component spinor and its derivatives at both interfaces numerically (Fig. 2, inset): two gaps of different size open at zero energy ( and ) and two gaps open at zero momentum where electron and hole states from different subgap bands cross (). This means that there is Cooper pair transport at zero energy as well as at the finite energies .

Cooper pair transport.—We use Fermi’s golden rule to calculate the Cooper pair current , where is the transition rate from an initial state with probability at energy to the final state with 2 more (less) electrons at energy . The tunnel Hamiltonian comprises the particle number non-conserving terms of the second-quantized perturbative model with electron operators and hole operators , where ,

 HT=∑χvnn′ks~Δnn′(k)scnχv,s(k)cn′−χv,−s(−k)+H.c. (4)

Because the superconductor interface has a finite width , we restrict the pairing amplitude in real space to . In momentum space (suppressing all indices) this amounts to with from the microscopic calculation and

 Δkk′ =~Δ(k−k′2)L2π∫dlsin[(l−k)w2](l−k)L2sin[(l+k′)w2](l+k′)L2,

where is the total length of the system, which does not enter the final results, and we have exploited that the integrand is peaked around . With this, the rates for removing (adding) a Cooper pair, , become , where at low temperatures the occupation probability with the voltage applied between the superconductor and the channel. Rewriting the sum over momenta as energy integrals, the current becomes

 I =32eℏπL2(2π)2∑nn′∫δμ−δμdε∣∣∂kK+n(ε)∂ε∣∣∣∣∂kK−n′(−ε)∂ε∣∣ ×|Δnn′(ε,−ε)|2. (5)

The combination of energy conservation and approximate momentum conservation implies that the pair tunneling probability has a single peak as a function of for each pair . Injection into the same subgap band, , happens near , and into different subbands, , near (Fig. 1b). Linearizing the dispersion (3) around these points sup (), , the tunnel amplitude becomes and we obtain the conductance

 G≈4G0∑nn′Tnn′[δw(δμ−εnn′0)+δw(δμ+εnn′0)], (6)

where is the conductance quantum, becomes the delta function for , and is the effective tunneling strength. Note that the conductance grows with the length of the interface. This is in contrast to conventional Cooper pair splitters, which suffer from an exponential suppression in the spatial size. The reason is that here Cooper pairs are split kinematically only after having tunneled locally into the channel, a process which can happen simultaneously along the whole interface. The conductance contains a central zero-bias peak and two characteristic side peaks [Fig. 3(a)], which arise because of the special subgap band structure and which correspond to the injection points marked in Fig. 1(b). The peak height is proportional to the induced superconducting pairings. A factor of 4 arises due to the spin and valley degeneracy and a factor of 2 due to pair transport.

Local Andreev reflection and Cooper pair splitting.—In Eq. (4) the singlet nature of the injected Cooper pairs is manifest. It is well established that Cooper pair splitting is closely related to CAR Prada and Sols (2004); Samuelsson et al. (2005). This applies if the dominant process is ordinary reflection. In our device, the 1D channel with a proximity-induced superconducting region is a NSN junction, in which only transmission through the S-region with amplitude and local Andreev reflection (an incoming quasiparticle in valley is reflected as an outgoing antiparticle with opposite velocity in valley ) with amplitude are possible. We consider the general scattering problem with finite CAR and normal reflection in the Supplemental Material sup (). When a voltage bias is applied between the superconductor and the channel to extract Cooper pairs, the incoming modes are filled with holes up to . Without the superconductor, all are transmitted and fill the outgoing modes up to , which is equivalent to a Fermi sea for electrons with the Fermi energy Samuelsson et al. (2003, 2005)

 \Ketδμ≡∏0<ε<δμαshα†s(ε)\Ket≡∏0<ε<δμαscαs(−ε)\Ket. (7)

Here, is the quasiparticle vacuum with respect to the Fermi level of the superconductor and creates outgoing holes with spin at energy in the left/right lead, which is the same as annihilating an electron with opposite spin at energy . We drop the valley index which is fixed by the requirement that outgoing modes move away from the superconducting region and the band index for simplicity sup (). Due to the proximity effect LAR becomes finite. The key observation is that when LAR occurs, no hole with spin at energy is transmitted to the other side. The outgoing mode is therefore occupied by a spin electron at energy [Fig. 3(b)]. To see this, we use Eq. (7) to write the outgoing state in terms of sup (),

 ∏0<ε<δμαs(stc¯α−s(−ε)+rcα†s(ε))c¯α†−s(−ε)\Ketδμ = ∏0<ε<δμα[t2+r2cα†↑(ε)c¯α†↓(−ε)cα†↓(ε)c¯α†↑(−ε) +rt(cα†↓(ε)c¯α†↑(−ε)−cα†↑(ε)c¯α†↓(−ε))]\Ketδμ. (8)

If is small, it becomes , where the desired nonlocal singlet state is explicit. This corresponds to a situation, where individual splitting events are well separated and it is meaningful to talk about pairs. In this regime of interest the perturbative result from the previous section holds. Only the emitted pairs contribute to the shot noise of the scattering state. In the opposite limit of perfect LAR with , , the outgoing state is a nonentangled product state. LAR is most pronounced at energies and [Fig. 3(a)] where the superconductor opens gaps in the spectrum for the case of an infinitely long () tunnel-junction (Fig. 2). The LAR process becomes weak for all energies, when falls below the coherence lengths .

Conclusion.—Our setup allows for highly efficient creation of nonlocal spin-entangled electrons without the need for repulsive interaction or energy filters. We note that the topological channel can be created electrically in the bulk of the BG sample, completely avoiding sharp sample edges, the main source of intervalley scattering Das Sarma et al. (2011), which could reduce the splitting efficiency. Moreover, using an electrically tunable channel geometry ballistic beamsplitters could be created to prove the spin entanglement via noise Burkard et al. (2000), so far an elusive goal. The spin relaxation and decoherence in BG are expected to be weak due to the small spin-orbit coupling van Gelderen and Smith (2010); Guinea (2010); Konschuh et al. (2012) and the sparsity of nuclear spins.

###### Acknowledgements.
We thank A. Baumgartner, P. Samuelsson, C. Schönenberger, and A. Levy Yeyati for helpful discussions and acknowledge support from the EU-FP7 Project SE2ND, No. 271554, the DFG, Grant No. RE 2978/1-1 and Research Training Group GrK1952/1 “Metrology for Complex Nanosystems”, and the Braunschweig International Graduate School of Metrology B-IGSM.

## Appendix A Supplemental Material

### a.1 Perturbation theory for the proximity effect

Here we give the details of the calculation of the proximity induced amplitudes . Because the superconductor is only weakly coupled through the tunnel barrier , we can derive an effective 1D model via low-order quasi-degenerate perturbation theory. We split into two parts, where is diagonal in the eigenbasis , and . Since is diagonal in spin and valley, we suppress the indices in the following. To first order in ,

 H(1)nn′(py)=∫dxΦ0,n†py(x)H1(x)Φ0,n′py(x) (A1)

and to second order,

 H(2)nn′(py) =12∑px∫dxΦ0,n†py(x)H1(x)~Φ0px,py(x) ×[(ε0,npy−~ε0px,py)−1+(ε0,n′py−~ε0px,py)−1] ×∫dx~Φ0†px,py(x)H1(x)Φ0,n′py(x), (A2)

where are the unperturbed free states above the gap with real and at energy . We impose the quantization condition and normalize the extended wavefunctions according to . The quantization length and the highest momentum are increased until the second order matrix elements converge. To study Cooper pair transport only the parts of , are relevant which are proportional to , i.e., they mix electron and hole states and therefore change the particle number. The relevant momenta are close to the crossing of the respective electron and hole band (see the discussion on approximate momentum conservation in the main text). This can involve one band, and , or both, , where are the Fermi points of the unperturbed dispersion, Eq. (3).

### a.2 Linearized subgap dispersion

The linearized subgap dispersion, Eq. (3), reads

 E≈−4321/4√V0α(χvpy∓2−3/4√V0/α) (B1)

around the Fermi points and

 E≈±V0√2−χv2−1/4√V0αpy (B2)

around zero momentum. The coefficients and used in the transport calculation, e.g., Eq. (6), can be read off immediately.

### a.3 Local and crossed Andreev reflection

In the most general case the incoming holes in a NSN junction can be transmitted , reflected , or undergo local ( or crossed () Andreev reflection. The outgoing state is

 ∏εsα(thhh¯α†s(ε)+tehc¯α†s(ε)+rhhhα†s(ε)+rehcα†s(ε))\Ket. (C1)

Rewriting the hole operators in terms of electron operators , and the Fermi sea in terms of the lowered Fermi sea as explained in Eq. (7) in the main text, we arrive at

 ∏εs[ (t2hh+r2hh)+(r2eh−t2eh)cR†s(ε)cR†−s(−ε)cL†s(ε)cL†−s(−ε) + s(rehrhh−thhteh)(cL†s(ε)cL†−s(−ε)+cR†s(ε)cR†−s(−ε)) + s(tehrhh−thhreh)(cR†s(ε)cL†−s(−ε)+cL†s(ε)cR†−s(−ε))]\Ketδμ. (C2)

The first line contains the product state contributions, the second line local pairs, and the third line nonlocal pairs. In the conventional reflection-dominated case, , realized in Y-junction Cooper pair splitters, the leading order contributions are

 [1+∑εs( rehs(cL†s(ε)cL†−s(−ε)+cR†s(ε)cR†−s(−ε)) + Missing or unrecognized delimiter for \Big (C3)

i.e., LAR produces local pairs and CAR produces nonlocal pairs. In the transmission-dominated situation, , the situation is reversed: the leading order is

 [1−∑εs( tehs(cL†s(ε)cL†−s(−ε)+cR†s(ε)cR†−s(−ε)) + Missing or unrecognized delimiter for \Big (C4)

so LAR produces nonlocal pairs and CAR produces local pairs. In the situation discussed in the main text, both CAR and reflection are forbidden, ruling out local pairs to all orders, as long as the valley symmetry is obeyed. Generally speaking it is undesirable to have simultaneously strong ordinary reflection and strong LAR or to have simultaneously strong transmission and strong CAR to build a Cooper pair splitter useful to create spin entanglement.

### a.4 Multiband Andreev reflection

The notation becomes more cumbersome, when both subgap bands are considered but the considerations are completely analogous. Without superconductivity the outgoing scattering state is

 ∏0<ε<δμαsnhα†ns(ε)\Ket≡∏0<ε<δμαsncαns(ε)\Ket≡\Ketμ, (D1)

where is the band index. In the presence of the superconductor, the transmitted holes can change the subgap band from to with an amplitude . Like in the one-band case, whenever the energy of an incoming electron is such that the spectrum of the S region has a gap, the transmission amplitude is exponentially suppressed with the length of the proximity region, and due to unitarity there is a finite amplitude for the spin- hole to be Andreev reflected locally as a spin- electron at energy :

 \Ketout =∏0<ε<δμαns∑m(stmnc¯αm,−s(−ε)+rmncα†ms(ε))\Ket =∏0<ε<δμαns∑m(stmnc¯αm,−s(−ε)+rmncα†ms(ε))∏m′c¯α†m′,−s(−ε)\Ketδμ =∏0<ε<δμαs[(t12t21−t11t22)+(r11r22−r12r21)cα†1s(ε)cα†2s(ε)c¯α†1,−s(−ε)c¯α†2,−s(−ε) +∑nm(−1)m(rn1t¯m2−rn2t¯m1)scα†ns(ε)c¯α†m,−s(−ε)]\Ketδμ. (D2)

Writing out the spin part of the product explicitly, the nonlocal singlet nature of the injected Cooper pairs becomes obvious:

 \Ketout =∏−δμ<ε<δμ[(t12t21−t11t22)2 +(t12t21−t11t22)∑nm(−1)m(rn1t¯m2−rn2t¯m1)(cL†n↑(ε)cR†m↓(−ε)−cL†n↓(ε)cR†m↑(−ε)) +O(r2)]\Ketδμ. (D3)

The higher order terms in contain multiple Cooper pairs and are not necessarily entangled, e.g., the contribution is a pure product state in which all states in the left/right lead at energy are occupied.

### References

1. P. Recher, E. V. Sukhorukov,  and D. Loss, Phys. Rev. B 63, 165314 (2001).
2. G. Lesovik, T. Martin,  and G. Blatter, Eur. Phys. J. B 24, 287 (2001).
3. P. Recher and D. Loss, Phys. Rev. B 65, 165327 (2002).
4. C. Bena, S. Vishveshwara, L. Balents,  and M. P. A. Fisher, Phys. Rev. Lett. 89, 037901 (2002).
5. P. Recher and D. Loss, Phys. Rev. Lett. 91, 267003 (2003).
6. A. Levy Yeyati, F. S. Bergeret, A. Martin-Rodero,  and T. M. Klapwijk, Nat. Phys. 3, 455 (2007).
7. J. Cayssol, Phys. Rev. Lett. 100, 147001 (2008).
8. K. Sato, D. Loss,  and Y. Tserkovnyak, Phys. Rev. Lett. 105, 226401 (2010).
9. L. Hofstetter, S. Csonka, J. Nygård,  and C. Schönenberger, Nature (London) 461, 960 (2009).
10. L. G. Herrmann, F. Portier, P. Roche, A. Levy Yeyati, T. Kontos,  and C. Strunk, Phys. Rev. Lett. 104, 026801 (2010).
11. J. Schindele, A. Baumgartner,  and C. Schönenberger, Phys. Rev. Lett. 109, 157002 (2012).
12. A. Das, Y. Ronen, M. Heiblum, D. Mahalu, A. V. Kretinin,  and H. Shtrikman, Nat. Commun. 3, 1165 (2012).
13. Z. B. Tan, D. Cox, T. Nieminen, P. Lähteenmäki, D. Golubev, G. B. Lesovik,  and P. J. Hakonen, Phys. Rev. Lett. 114, 096602 (2015).
14. R. S. Deacon, A. Oiwa, J. Sailer, S. Baba, Y. Kanai, K. Shibata, K. Hirakawa,  and S. Tarucha, Nat. Commun. 6, (2015).
15. D. Loss and D. P. DiVincenzo, Phys. Rev. A 57, 120 (1998).
16. J. Torrès and T. Martin, Eur. Phys. J. B 12, 319 (1999).
17. G. Falci, D. Feinberg,  and F. W. J. Hekking, Europhys. Lett. 54, 255 (2001).
18. I. A. Sadovskyy, G. B. Lesovik,  and V. M. Vinokur, ArXiv:1412.8145.
19. W. Chen, R. Shen, L. Sheng, B. G. Wang,  and D. Y. Xing, Phys. Rev. B 84, 115420 (2011).
20. R. W. Reinthaler, P. Recher,  and E. M. Hankiewicz, Phys. Rev. Lett. 110, 226802 (2013).
21. J. Nilsson, A. R. Akhmerov,  and C. W. J. Beenakker, Phys. Rev. Lett. 101, 120403 (2008).
22. J. Wang, L. Hao,  and K. S. Chan, Phys. Rev. B 91, 085415 (2015).
23. W. Chen, R. Shen, L. Sheng, B. G. Wang,  and D. Y. Xing, Phys. Rev. Lett. 109, 036802 (2012).
24. I. Martin, Y. M. Blanter,  and A. F. Morpurgo, Phys. Rev. Lett. 100, 036804 (2008).
25. F. Zhang, A. H. MacDonald,  and E. J. Mele, Proc. Nat. Acad. Sci. USA 110, 10546 (2013).
26. C. W. J. Beenakker, Rev. Mod. Phys. 80, 1337 (2008).
27. A. Rycerz, J. Tworzydlo,  and C. W. J. Beenakker, Nat. Phys. 3, 172 (2007).
28. L. Ju, Z. Shi, N. Nair, Y. Lv, C. Jin, J. Velasco Jr, C. Ojeda-Aristizabal, H. A. Bechtel, M. C. Martin, A. Zettl, J. Analytis,  and F. Wang, Nature (London) 520, 650 (2015).
29. E. McCann and V. I. Fal’ko, Phys. Rev. Lett. 96, 086805 (2006).
30. R. Winkler, Spin-Orbit Coupling Effects in Two-Dimensional Electron and Hole Systems (Springer-Verlag, Berlin, 2003).
31. See Supplemental Material at … for further details .
32. F. Parhizgar and A. M. Black-Schaffer, Phys. Rev. B 90, 184517 (2014).
33. E. Prada and F. Sols, Eur. Phys. J. B 40, 379 (2004).
34. P. Samuelsson, E. V. Sukhorukov,  and M. Büttiker, New Journal of Physics 7, 176 (2005).
35. P. Samuelsson, E. V. Sukhorukov,  and M. Büttiker, Phys. Rev. Lett. 91, 157002 (2003).
36. S. Das Sarma, S. Adam, E. H. Hwang,  and E. Rossi, Rev. Mod. Phys. 83, 407 (2011).
37. G. Burkard, D. Loss,  and E. V. Sukhorukov, Phys. Rev. B 61, R16303 (2000).
38. R. van Gelderen and C. M. Smith, Phys. Rev. B 81, 125435 (2010).
39. F. Guinea, New Journal of Physics 12, 083063 (2010).
40. S. Konschuh, M. Gmitra, D. Kochan,  and J. Fabian, Phys. Rev. B 85, 115423 (2012).
104161