Large twisting angles in Bilayer graphene Moire quantum dot structures
LLarge twisting angles in Bilayer graphene Moire quantum dot structures
Jozef Bucko
1, 2 and Frantiˇsek Herman
1, 3 Institute for Theoretical Physics, ETH Zurich, CH-8093, Switzerland Institute for Computational Science, University of Zurich,Winterthurerstrasse 190, 8057 Zurich, Switzerland Department of Experimental Physics, Comenius University,Mlynsk´a Dolina F2, 842 48 Bratislava, Slovakia (Dated: January 29, 2021)Recent exploration of the commensurate structure in the turbostratic double layer graphene showsthat the large angle twisting can be treated by the decrease of the effective velocity within the energyspectra of the single layer graphene. Within our work, we use this result as a starting point, aimingtowards understanding the physics of by a large angle twisted double layer graphene (i.e. Moire)quantum dot systems. We show that within this simple approach using the language of the firstquantization, yet another so far unrealized (not up to our knowledge), illustrative property of thecommutation relation appears in the graphene physics. Intriguingly, large twisting angles show to bea suitable tunning knob of the position symmetry in the graphene systems. Complete overview ofthe large angle twisting on the considered dot systems is provided.
I. INTRODUCTION
Graphene, elegant honeycomb structured atomicmonolayer material, has already shown to have inter-esting mechanical as well as electronic propertiesof the ideal semi-metal. Due to its linear dispersion,related mathematical description as well as resultingphysical properties, understanding and describinggraphene unquestionably belongs to the showcase of thetwenty-first century physics.Adding another graphene layer shows to be animportant step towards the simple way of creating thegap in the energy spectrum . This natural progress inthe development of graphene systems also opens the fieldof twistronics by introducing another degree of freedomthrough the possibility of twisting the layers towardseach other . New macroscopically tunable length scaleemerges, resulting from the appearing structure of theMoire pattern. Its presence leads to lowering of the slopeof the energy band close to the valley explained bythe reduction of the Fermi velocity . Also, additionalDirac points as well as flattening band at the magicangle θ ∗ ≈ . ◦ occur. Described tuning mechanismshows to lead towards interesting and novel propertiesof the underlying density of states in the normal state ,as well as to influencing e.g. the value of T c in thesuperconductive state .Area of twistronics becomes also attractive from thepoint of view of dot and flake (mesoscopic) systemsbounded by the additional confining potential, whichkeeps the electron confined on the suitable length scaleunder ∼ nm . Dot twistronics might be even moreattractive from the application point of view, since we donot run into the dimension vs. quality of the preparedsystems kind of problems . Simple exploration of theenergy scale of the created bound states (theoreticallyand experimentally), their behavior with changing char- acteristic length, twisting angle, as well as application ofthe external magnetic field create standard tools in thedeveloping area of the dot twistronics .In what follows, we introduce reduced velocity in theconnection with the single layer graphene quantum dotmodel, in order to simulate the effect of twisting layers ofthe dot about large angles. In the next sections, we willhave a look at the effect of twisting angles on the scalingof the states, bending of the energy levels and effect of thelarge twisting angles on the wave-function of the twistedbilayer graphene Moire quantum dot (BGM-QD) system. II. (LARGE ANGLE) TWISTED BILAYERGRAPHENE QUANTUM DOTS
FIG. 1. Moire pattern situation close to the boundary ofthe cartoon 7 nm wide BGM-QD system considering differentvalues of the large twisting angle θ . In the following sections, we focus on the semianalytic a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n approach towards the explanation of the generic behaviorof the energy bands of the twisted bilayer graphenequantum dot, considering large twisting angles θ . Recentnumerical studies published in Ref. 10 and 11 arebased on tight-binding approximation and suited as anenlightening motivation of our work. Mentioned purelynumerical approach has its limitations in computationalpower and therefore usually allows one to focus ondot-systems considering radiuses up to ∼ nm .Our semianalytical approach, which allows us toexamine general physical phenomena considering largersystems, is motivated by the results of the secondmentioned study . Its authors claim that consideringthe interval of large twisting angles 10 ◦ (cid:46) θ (cid:46) ◦ ,interaction between the two single sheets of graphenelayers is small and the band structure basically copiesthe one of the two single layer sheets.In the Fig. 1, we plot the Moire pattern structurefor few chosen values of θ , together with two differentimportant length scales considered in our situation. R marks the radius of our dot, while D corresponds to theMoire length scale. From now on, we properly define theterm BGM-QD (bilayer graphene Moire quantum dot)as regime when R (cid:38) D (i.e. when the Moire pattern isclearly recognizable in the region of the dot). A. Effective velocity
In fact, our simple approach combines the single layerquantum dot model solved in the Ref. 12 together withthe effective Fermi velocity close to the minimum of eachvalley (already measured in the Ref. 6) within the secondorder of the perturbation theory suited to the crystal andelectronic structure of the turbostratic graphene : v F = v SLF (cid:16) − α ∆ K (cid:17) , (1)where α represents the binding constant element betweenstates in individual graphene layers. Single layer Fermivelocity is represented by v SLF . The distance of valleysfrom the individual graphene layers ∆ K , displayed in theFig. 2, is in the Moire structures tunable variable anddepends on the twisting angle θ by :∆ K = 43 a sin θ , (2)where a = 0 . nm is the shortest distance between theCarbon atoms in the graphene.Naturally, since the Moire period D can be easily con-nected together with the twisting angle θ according to : D = √ a csc θ , (3) FIG. 2. Brillouin zone of the unrotated (red), rotated (blue)single layer graphene, as well as twisted bilayer graphene(black). ∆ K relates with the Moire period by simple equation∆ K.D = 2 / √
3. For illustration, we plot the standardpicture of two mutually rotated Brillouin zones consider-ing two graphene sheets in the Fig. 2. One (light red)corresponds to the unrotated BZ of the graphene layer,while the other (light blue) corresponds to the rotatedone. The picture describes the situation corresponding tothe large twisting angles θ ≈ ◦ . In the Fig. 2, we alsoplot the individual valleys of the unrotated ( K U , K ∗ U )as well as the rotated layer ( K R , K ∗ R ) in the Brillouinzone corresponding to the twisted bilayer graphene (andtherefore also BGM-QD) according to the mappingdescribed in the Ref. 7.Within our treatment, we consider interval of twistingangles 5 ◦ (cid:46) θ (cid:46) ◦ , where the approximation for theeffective reduced velocity in the form of the Eq. (1) holds.This assumption is also reasonable, once we considerour condition R (cid:38) D , together with the Eq. (3) and R min = 10 nm . After the simple calculation we findthe minimal value of the twisting angle θ min = 1 . ◦ ,which is way bellow our analytical approximation ofthe reduced velocity. Larger values of the radius R lead to even smaller values of the θ min . Notice, thatconsidering cases of the large twisting angles θ ≈ ◦ ,∆ K is close to its maximum and the correction in thereduced Fermi velocity corresponding to the secondorder of the perturbation theory is therefore fully justified.On the contrary, considering small twisting angles, thesystem starts to be dominated by the flat dispersion be-havior coming from the merging of the two Van Hovesingularities , which makes it unsuitable for our de-scription by the model with the linear dispersion. Also,it is good to emphasize that the idea of using effectivevelocity in the form of the Eq. (9) allows us (in the ap-proximation of large θ ) to solve two problems (two layers)at once. Therefore, our solution for the energy bands willautomatically have 2 × × spin. B. Hamiltonian
After the proper introduction of the efective velocity v F ( θ ) from the Eq. (1), now it is time to introduce itto our effective Hamiltonian valid close to the Diracpoints. Since the energy scale that we are focusing onis much smaller than the Fermi energy for graphene, wecan allow ourselves to work within the assumption of thelinearized tight binding approximations close to the Diracpoints. The resulting differential equation of motion isthen coming from the usage of the first quantized language,where the momentum coordinates are exchanged by thedifferential operators of derivatives in position. Next, wewill also consider rotational symmetry of the confiningpotential U ( r ) = 0, assuming r ≤ R resp. U ( r ) = U ,assuming r > R : H τ ( θ ) = v F ( θ ) ( p + e A ) . σ + τ ∆ σ z + U ( r ) , (4)where p is the momentum operator, σ = ( σ x , σ y ) isthe vector of the Pauli matrices and A is the vectorpotential related to the magnetic field in the z -directionby A = B/ − y, x, , where τ differentiates thevalleys.To avoid any confusion, let us discuss monolayerstructure of the Hamiltonian in the Eq. (4). As wehave already mentioned, in the regime of large twistingangles, we are in the situation, where the limit of theBGM-QD being represented as two single layer sheetsof graphene as well as the introduced treatment ofthe reduced velocity are valid. Therefore, the Eq. (4)models twice the same circular monolayer graphene flakestructure (corresponding to the twisted bilayer) andthus, effectively, we solve the same problem twice (twolayers, two Hamiltonians, two sets of Dirac equations).In principle, one can also solve the problem as set of twoDirac equations (each corresponding to separate layer),while considering two different gap terms. First suited toone layer and the second to another layer. The differencewould be at the end in the splitting of the degeneracy ofthe states coming from two layers in our case. Otherwise,the qualitative features of the provided solution would bethe same as the ones described in the next chapters.Therefore, our model corresponds to the (theoretically)simplest realization where we assume using the substratefor one graphene flake (applied from one side of thegraphene layer) and then the same substrate from theother side of the second graphene layer. Our two twistedgraphene flakes would be therefore sandwiched by thesame substrate, which would create the same gap in bothof them. It is fair to say, that even though this modelrepresents the simplest theoretical realization, it may be experimentally challenging.Important thing to mention and also make clear is thedifference between the geometry of the gap and boundingpotential of the dot in our model compared to already dis-cussed numerical tight-binding studies. In our approach(which is close to the one used in the Ref. 12) the gap isrealized in the whole region of the dot and the confiningpotential U ( r ) is nonzero only for r > R . However in thenumerical study in Ref. 11, the dot is being created bythe gating of the region outside of the dot. C. Uncertainty principle
In the following, let us consider following commutatorof the Hamiltonian (4) (considering value of the vectorpotential A = ) together with the position operator r :[ H τ ( θ ) , r ] = v F ( θ ) [ p , r ] · σ , = − i (cid:126) v F ( θ ) σ , (5)where the well known canonical commutation relation[ p a , r b ] = − i (cid:126) δ ab was used. In general, results of thecommutators are interesting in two cases. First, oncethey are zero and we start to think about commoneigenstates of the considered operators (together withtheir corresponding eigenvalues), symmetries andconservation laws (if one of the operators is Hamiltonianof the system). Or second, once the commutator isconstant and therefore we have uncertainty principleavailable.Notice, that the Eq. (5) offers a compromise of thesetwo properties. Notice also, that the twisting angle θ entering the reduced effective velocity works as a tuningknob of the deviance from the exact symmetry of theposition operator of the Hamiltonian. In such a way, dueto the unique properties of the graphene Hamiltonian,we have non-trivial commutation relation where themacroscopically accessible parameter tunes the quantumproperties of the system. If (assuming theoretical limit) v F ( θ ) = 0, then the mentioned operators have commoneigenstates with their corresponding eigenvalues. Onewould expect that the rising constant on the right sideof the Eq. (5) will measure deviance from this case.Since the Hamiltonian will still have an exact energyeigenvalue due to Dirac equation, the state will be nolonger precisely set in position. In simpler words, Eq. (5)means that the slower Dirac electron is the one which ismore localized.If we think about the overall effect of the twisting angle θ in the bilayer graphene, we can easily qualitatively sumit up by saying, that the large angle (macroscopicallytunable) twist corresponds to the creating additionalpotential rising from the new periodic length scale D .Treating this additional potential on the level of thesecond order perturbation theory assuming two sheetsof graphene shows, that the overall result leads towardsslowing down (reduction of the velocity) of the Diracelectron once we assume decreasing twisting angles fromthe value of θ = 30 ◦ .On the top of this interesting property, we can stillexplore underlying uncertainty principle related to theEq. (5). Let us remind the Schr¨odinger uncertaintyrelation suitable for our situation:∆ E ∆ r ≥ (cid:115)(cid:18) (cid:104){ H, r }(cid:105) − (cid:104) H (cid:105)(cid:104) r (cid:105) (cid:19) + (cid:18) i (cid:104) [ H, r ] (cid:105) (cid:19) , (6)where ∆ E = (cid:104) ψ ( H − E ) | ( H − E ) ψ (cid:105) and ∆ r = (cid:104) ψ ( r − (cid:104) r (cid:105) ) | ( r − (cid:104) r (cid:105) ) ψ (cid:105) . Now, thanks to the Hamil-tonian being one of the operators, the equation of motion(Dirac equation) H | ψ (cid:105) = E | ψ (cid:105) and also the fact that inour considered system we assume radial symmetry (cid:104) r (cid:105) = 0,we get: ∆ E ∆ r ≥ . (7)Notice, that this uncertainty principle allows us inprinciple to have well localized states together with theexact values of energy. All that due to the symmetry ofour system and the equation of motion.To be more general and also a bit pedagogical, Eq. (7)is exactly the reason why all of the particles in the boxes(where (cid:104) r (cid:105) = 0) etc. are so grateful objects of our focusin quantum mechanics. All of these systems allow us tohave discrete energy levels, together with the localizedstates on the reasonable length scales. D. Reduced energy scale
However enlightening we can find the Eq. (7) to be inthe strictly mathematical sense, in the real life experimen-tal conditions we would expect ∆ r to be on the same scaleas the radius of the dot R . Together with the Eq. (5), wecan get an estimation for the energy level resolution:∆ ER ≈ (cid:126) v F ( θ ) . (8)This equation makes complete sense together withthe Eq. (5), under which we would expect that oncethe velocity reduces, the system will go closer to theclassical one, with more dense energy levels (closer tothe continuum since the position operator is alreadycontinuous).So, the other idea of this part can be summarized bythe statement that the reduced Fermi velocity, whichleads to the reduction of the energy scale occurring in thegraphene physics: (cid:126) v SLF R → (cid:126) v F ( θ ) R , (9) affects the behavior of the allowed energies as well asstates in the dot problem with the twisting angle θ . III. FORMULATION OF THE EIGENPROBLEM
After the short discussion focused on the physics of theconsidered system, based on the analysis of the uncer-tainty principle, let us now introduce equations of motiontogether with the corresponding boundary condition.Mesoscopic system of the BGM-QD considering largetwisting angles is on the scale of ∼ nm in space and ∼ meV in energies.Based on our previous discussion related to the Hamil-tonian defined by the Eq. (4), Dirac equation of motiondescribing our BGM-QD system reads: H τ ( θ )Ψ τ ( θ ) = E τ ( θ )Ψ τ ( θ ) , (10)where the Brillouin zone of the system has two sites andwe order the components of the two-site envelope wave-function in the following way :Ψ τ ( θ ) = (cid:18) Ψ τA ( θ )Ψ τB ( θ ) (cid:19) . (11)Let us remind, that we explicitly wrote down the depen-dence on the new external parameter of the twisting angle θ . The two-site envelope wave function (due to cylindricalsymmetry) can be factorized as:Ψ τ ( r, ϕ ) = e imϕ √ r (cid:18) e − iϕ (cid:19) Ψ τ ( r ) , (12)where r and ϕ are the polar spatial coordinates andΨ τ ( r ) is the remaining part of the envelope wave-function.In the Appendix A, we provide further explanationand simplification of the introduced eigenvalue problemdescribed by the Eq. (10) up to the complete recipe ofthe half analytical (angular part) and half numerical(radial part) solution.After the introduction of the considered eigenvalueproblem, let us emphasize, that the final eigenenergies aswell as final form of the allowed eigenvectors can be foundfrom the solution of the boundary condition (at r = R )using states inside ( < ) and outside ( > ) of the dot:Ψ τ< ( R, ϕ ) = Ψ τ> ( R, ϕ ) . (13)The difference between solutions inside and outside of thedot is of course coming from the bounding potential U ( r )defined by the Eq. (4), which is zero (nonzero) inside(ouside of the dot). IV. RESULTS AND DISCUSSION
After introducing realized model as well as describingits solution, we can finally address the effect of the reducedeffective velocity on the considered energy levels, theirbehavior in the magnetic field and also their angulardependence. At the end, we look at the behavior ofthe considered wave function evolving together with thetwisting angle θ . For clarity, we will focus only on thecase considering m = 0. A. State scaling
In the Fig. 3, we plot the R -scaling of the states consid-ering two different values of the twisting angle θ . Notice,that the twisting angle θ does not change the R -evolutionof the states, but works as a tuning knob of the distancebetween the energy levels. Provided numerical solutiontherefore agrees nicely with our understanding gainedfrom the uncertainty principle formulated in the Sec. II C. FIG. 3. Energy levels of the BGM-QD considering m = 0,potentials U = ∆ = 10 δ , where δ = (cid:126) v F ( θ ) /R , two differentvalleys τ = ± − θ = 5 ◦ (red and blue) and θ = 30 ◦ (green and orange). In the very rough approach the energy levels can beapproximated by the formula: E = (cid:126) v F ( θ ) g ( m, τ ) /R (wewill see that the situation is slightly more complicated).Function g ( m, τ ) is discrete function of m as well as τ obeying all symmetries already discussed in the Ref. 12and energy states have clear 1 /R scaling. Notice also,that in the limit of the small radiuses, our rather generalsemi-analytical analysis qualitatively agrees and also gen-eralizes the one realized and numerically calculated in theRef. 11. B. Effect of the magnetic field on the states andenergy levels
Let us start by the Fig. 4, where we compare evolutionof the energy levels in the magnetic field consideringtwo different values of the twisting angles θ = 5 ◦ and θ = 30 ◦ . In the overall paragraph, we use dimensionlesssuitable units of the fraction R/l B , where l B = (cid:112) (cid:126) /eB .For clarity, we focus just on the states with m = 0 and τ = ±
1. As we can clearly see, the main difference resultsfrom the already discussed scaling of the states with thetwisting angle θ . This effect then represents itself in thesmaller distances between corresponding energy levelsassuming lower values of the twisting angle θ . a)b) FIG. 4. Energy levels in the magnetic field considering twistingangles a) θ = 5 ◦ and b) θ = 30 ◦ .
1. Wave-function
Let us now focus on the effect of the magnetic field onthe wave-function of the underlying state. In the Fig. 5,we plot the BGM-QD wave-function at different valuesof ρ = R/ √ l B . These units are motivated by the ideaof counting the allowed states N inside the dot, sinceaccording to Ref. 18 N = R / l B . We can immediatelynotice, that considering rising ρ the state becomes morelocalized in the centre. This property is in completeagreement with the logic that the increasing magneticfield allows us to have a better localized state inside thedot with the radius R . We know, that this is true dueto the uncertainty principle based on the commutationrelation between coordinates of the orbit in the magneticfield with the center at coordinates ( X, Y ) :[ X, Y ] = 2 il B , ∆ X ∆ Y = l B . (14)I.e. the weight of the wave-function lies therefore moreand more inside the dot. Considering increasing valueof ρ going towards the limit ρ (cid:29)
1, the creation of theLandau levels (plotted already in the Fig. 4) allows us tohave another area of higher probability density localizedaround the boundary of the dot.Fig. 5 also opens discussion in the direction of theapplicability of the graphene quantum dots as promisingcandidate for the quantum bit material. It shows, thatin the quantum dots being created by the gates where ∆(or voltage between layers V , if we think about doublelayer graphene system ) is on the same scale as thebounding potential U , huge part of the wave-function isbeing localized outside of the dot. As we can see, applyingperpendicular magnetic field (considering right amplitude)can lead to better localization of the particle inside thedot itself and therefore to more stable states with higherdecoherence times. FIG. 5. Wave function of the twisted bilayer quantum dot asa function of the increasing ratio of R and magnetic length l B : ρ = R/ √ l B . Parameters related to the considered model: m = 0, τ = 1, U = 2 δ , ∆ = 2 δ , R = 25 nm , θ = 30 ◦ . C. Effect of the twisting angle on the states andenergy levels
In the Fig. 6 we plot the energy levels as a function ofthe twisting angle θ . Notice, that qualitatively we obtainexpected behavior of decreasing energy scale towardssmaller angles (bending bands). One can easily notice,that the decreasing trend can not be explained purelyby the decrease of the effective velocity v F ( θ ) from theEq. (1) in the form described by the Eq. (9). However, FIG. 6. Evolution of the states in the conduction band asa function of the twisting angle θ considering both valleys τ = ± each of the curves plotted in the Fig. 6 can be in facteasily fitted by an empirical formula: ε (cid:48) = ε (cid:18) − γ csc (cid:18) θ (cid:19)(cid:19) , (15)considering fitting parameters ε and γ . This formulashould not be in fact surprising, since it reminds us thecombination of the mentioned formulae Eq. (1), (2) (whichare valid up to second order of the perturbation theory),together with the free binding parameter between thestates in the separate layers. Since the overlap betweendifferent states and Hamiltonian correction will be differ-ent, it makes sense that γ will differ. Notice also that thetwisting in our considered system respects the symmetryof the system, therefore we have elegant tuning knobfitted to the geometry of the dot.
1. Wave-function
On the other side, let us focus on the qualitative discus-sion of the effect of the decreasing velocity on the spatialdistribution of the wave function. Realizing Eq. (5), wecan argue that lowering twisting angle renormalizes theeffective velocity towards lower values and the state be-comes better localized and also slightly shifted towardsthe boundary of the dot.
FIG. 7. Wave function of the twisted bilayer quantum dot asa function of the twisting angle θ . Parameters related to theconsidered model: m = 0, τ = 1, U = ∆ = 2 δ , R = 25 nm , ρ = 1, and the twisting angle is changing from θ = 4 ◦ (topcurve) to θ = 20 ◦ (bottom curve). In the Fig. 7, we plot the probability distribution inthe radial direction, considering finite magnetic field( B = 2 . T ), in order to have better localized state. Thevalue of the magnetic field B is chosen in the way that wehave exactly one allowed state in the dot with the radius R , since ρ = 1. Notice, that slightly increased localizabil-ity of the realised state towards the boundary of the dot asa function of the decreasing θ is also in qualitative agree-ment with the numerical analysis published in the Ref. 10.Let us emphasize at the end of this subsection, that ournumerically obtained results, focusing on the behaviorof the energy states as well as envelope wave-function inthe BGM-QD system under the large twisting angles θ nicely match with the expectations based on the approachdiscussed in the Subsections. II C and II D.
2. Comment on the small twisting angles
From what we know so far, we can qualitativelyaddress the question of the small twisting angles θ , or i.e.situations when R ≈ D . In the cases, when R (cid:29) D thedot system is fully described either by AA or AB stackingtogether with their states. In the cases, when R starts tobe comparable with D , based on our calculations we knowthat the wave function can be localized in the regionclose to the boundary of the dot. On the other side,described region also corresponds to the huge mixture ofcommensurate and incommensurate parts of double layeralignment (partially, we can observe such a situationin the subfigures a) and d) of the Fig. 1). Therefore,the solution of the problem will contain signatures ofsolutions of the exact AA or AB stacked limit cases and also signature of our obtained solution considering limit R (cid:29) D .No wonder, that in this intermediate regime, the evo-lution of the energy levels as a function of θ will dependon the exact (from the direction of the limit values of thetwisting angles θ = { ◦ , ◦ } ) AA or AB stacked state, aswell as on the exact geometry of the twisted bilayer (fromthe limit cases of the solution realising the BGM-QD).Since the geometry on the boundary will be changingabruptly with small changes in angle θ , we should expectlarge changes in the evolution of the energy levels. No-tice, that this kind of behavior is visible in the alreadymentioned tight binding study in Ref. 11 (however notrecognized or highlighted by the original authors). Noticealso, that considering the angles with diminishing Fermivelocity, assumption of the linearized dispersion (whichwe stand on) fails completely. V. SUMMARY
Since the numerical results from the solution of theboundary condition considering by a large angle twistedbilayer graphene quantum (Moire) dot problem werealready discussed in the previous section, let us focuson its synthesis aiming for a better understandingof this highly tunable mesoscopic system. The mainmessage of this letter is therefore following: twisting theBGM-QD system in the regime of large angles shows tobe suitable tool in order to get better localized states.Such a state also has lower energy due to reduced velocity.Added value of our work is also in the effort toshed more analytical light into the field of twistron-ics, which is currently dominated by the numericalstudies based on the first-principle calculations alreadydiscussed in the referenced material . As we cannotice, some of the behavior of the energy levelsas well as corresponding states can be examined andguessed by simple analytic analysis provided in the Sec. II.Use of the better localized states is natural in thetunneling applications (and references therein).Energy levels distance on the infrared scale rangingfrom δ ≈ meV (considering R = 100 nm ) already upto δ ≈ meV (considering R = 1 nm ) craves for theapplication in the infrared spectroscopy area. Alreadydiscussed energy scale of few meV also opens questionsin the direction of exploration of superconductivephenomena on quantum dot systems, since the systemitself is (with the considered length scale) able to host aCooper pair. In this way, it shows that the technologyof the gate bounded quantum dots enriches already vastarea of the quantum dot research .It is well known, that one of the problems of grapheneregarding commercial applications is production of thelarge sheets of this material. However, small flakessuitable for the quantum dots are not a problem at all.Let us just mention recent development in the field ofgraphene dot physics towards the identification of thesuitable two dimensional dot flakes and let us connectthis potential with the progress of the understanding ofthe measured states . We should also definitely notomit successful experimental effort in the fabrication ofthe large angle twisted samples motivated by the beautyof the dodecagonal quasicrystal pattern appearance at θ = 30 ◦ , and its richness on the Dirac cone replicasstructure.All of the mentioned ideas show promising researchfuture in the field of the twisted bilayer graphenequantum dot physics. Our modest contribution to thefast developing field lays in the direction of the verysimple description of the tuning by the large angletwisting of two graphene flakes placed above each other.It shows, that this effect leads to the decrease of theelectron velocities. Described effect causes severalfeatures in the energy and state description including e.g.higher localization of the envelope wave-function. ACKNOWLEDGEMENTS
We are grateful for financial support provided by theSlovak Research and Development Agency under Con-tract No. APVV-19-0371 and by the agency VEGA underContract No. 1/0640/20. We are also grateful for finan-cial support from the Swiss National Science Foundationthrough Division II (Grant No. 184739). We are alsograteful to Eliska Greplova, Peter Rickhaus, Fokko deVries and Frank Sch¨afer for many useful and interestingcomments and discussions.
Appendix A: Solution strategy
In this Appendix, we further simplify the set of thedifferential equations defined by the Eq. (10) usingassumption of the polar symmetry of the confiningpotential U ( r ), as well as further matrix structure of theconsidered equations. At the end of the Appendix, wealso provide full explanation of the boundary conditiondefined by the Eq. (13).Due to the rotational symmetry of our problem, letus formulate the kinetic part of H τ ( θ ), which applies toΨ ( r ), using polar coordinates: H ( θ ) = (cid:126) v F ( θ ) i √ l B (cid:32) ∂ ξ − m − / ξ − sξ∂ ξ + m − / ξ + sξ (cid:33) , (A1)where we introduced dimensionless units of ξ = r/ √ l B and magnetic length l B = (cid:112) (cid:126) /eB , which is at the value of B = 1 T on the scale of l B ≈ nm .Two first order equations for components of Ψ ( r ) haveknown solutions in terms of the iterative relations :( ∂ ξ − ( m − / /ξ − sξ ) φ sm − = a s φ sm , ( ∂ ξ + ( m − / /ξ + sξ ) φ sm = a s φ sm − . (A2)The functions φ sm as well as coefficients a si differ insideand outside the dot. φ sm are proportional to confluenthypergeometric (also known as Kummers) functions: φ sm + α ( ξ ) ≡ e − ξ / ξ b − M ( a, b, z ) / Γ( b ) , r ≤ Re − ξ / ξ b − U ( a, b, z ) , r > R (A3)in which we used standard notation for confluent hy-pergeometric functions of the first M ( a, b, z ) and second U ( a, b, z ) kind. In our case: a = | m + α | + 1 + s ( m − − α )2 + κ ,b = 1 + | m + α | ,z = ξ . (A4)As for the coefficients a si defined in the Ref. 28 inside thedot: a s = κ / , m ≥ , m = 02 , m ≤ − a s = , m ≥ κ / , m = 0 κ / , m ≤ − r > R we can observe that the coefficients are inde-pendent from angular momentum m : a s = − [( s + 1) + κ (1 − s ) / ,a s = − [(1 − s ) + κ (1 + s ) / . (A6)Next, we can further factorize two-site envelope wavefunction Ψ τ to introduce Ψ τ :Ψ τ = (cid:18) φ sm φ sm − (cid:19) Ψ τ . (A7)With this choice and using formulae (A2) – (A7) we canmanipulate the equations so that we can further simplifyand replace: H Ψ τ = H (cid:18) φ sm φ sm − (cid:19) Ψ τ = (cid:18) φ sm φ sm − (cid:19) (cid:126) v F ( θ ) i √ l B (cid:18) a s a s (cid:19) Ψ τ . After all these steps, we end up with the factorizationof the original two-site envelope wave-function:Ψ τ ( r, ϕ ) = e imϕ √ r (cid:18) e − iϕ (cid:19) (cid:18) φ sm ( r ) 00 φ sm − ( r ) (cid:19) Ψ τ . (A8)The radial dependence of the envelope function is incor-porated in φ si functions. On the other hand, Ψ τ is notcoordinate dependent. Moreover, interaction part of thehamiltonian commutes with both matrices in the aboveenvelope function expansion and thus we can move H τ inside the envelope function in front of Ψ τ . Then we areleft with the task: (cid:126) v F ( θ ) i √ l B (cid:18) a s a s (cid:19) Ψ τ + H τ Ψ τ = E Ψ τ , (A9)which is equivalent to solving the following homogeneoussystem: (cid:32) τ ∆ + U ( r ) − E a s (cid:126) v F ( θ ) i √ l B a s (cid:126) v F ( θ ) i √ l B − τ ∆ + U ( r ) − E (cid:33) Ψ τ = 0 . (A10)We have the condition for zero determinant in (A10) inorder to have nontrivial Ψ τ . Full envelope function canbe then easily recovered. However, we have two freeparameters to fix, namely κ (present in a si ) and E . Onecan be fixed by the condition on singularity of the abovematrix which yields the following relation between κ and E : κ <,> = 2 l B (cid:0) ∆ − ε <,> (cid:1) / ( (cid:126) v F ( θ )) . (A11)Throughout the task we set: U ( r ) = (cid:40) , r ≤ RU , r > R (A12) and therefore ε < = E , resp. ε > = E − U . The secondcondition to fix remaining degree of freedom comes fromthe continuity of two-site envelope wave-function at theboundary of the dot.In both cases (inside and outside) we have just one valueof κ . In this way we basically have solution inside Ψ τ ,< ( ξ )and outside Ψ τ ,> ( ξ ) of the dot up to a multiplicativeconstants (let us call them C < and C > ). We find theenergy of the respective state by matching these two at r = R . Condition on the nontrivial solution with nonzero C < and C > immediately leads to linear combination of therealized envelope wave-function. So, the solution whichdiscloses the allowed values of R (or B , or θ ) and E canbe found by solving for zero of the determinant: Det τ ( R, B, θ, E ) = | Ψ τ ,< ( R, B, θ, E ) , Ψ τ ,> ( R, B, θ, E ) | . (A13)Ratio of the considered constants can be easily found fromthe ratio of the equations in the boundary condition andthe last fix is coming from the normalization condition ofthe envelope wave-function. D. G. Papageorgiou, I. A. Kinloch, and R. J. Young,Progress in Materials Science , 75 (2017). A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov,and A. K. Geim, Review of Modern Physics , 109 (2017). H. Min, B. Sahu, S. K. Banerjee, and A. H. MacDonald,Physical Review B , 1 (2007). S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin, andE. Kaxiras, Physical Review B , 1 (2017). Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi,E. Kaxiras, and P. Jarillo-Herrero, nature , 43 (2018). P. Moon and M. Koshino, Physical Review B , 195458(2012). S. Shallcross, S. Sharma, E. Kandelaki, and O. A. Pankra-tov, Physical Review B , 165105 (2010). R. Bistritzer and A. H. MacDonald, PNAS , 12233 –12237 (2011). E. Greplova, C. Gold, B. Kratochwil, T. Davatz, R. Pisoni,A. Kurzmann, P. Rickhaus, M. H. Fischer, T. Ihn, andS. D. Huber, Phys. Rev. Applied , 064017 (2020). A. Tiutiunnyka, C. Duqueb, F. Caro-Loperac, M. Mora-Ramosa, and J. Correa, Physica E: Low-dimensional Sys-tems and Nanostructures , 36 (2019). M. Mirzakhani, F. M. Peeters, and M. Zarenia, PhysicalReview B , 075413 (2020). P. Recher, J. Nilsson, G. Burkard, and B. Trauzettel,Physical Review B , 085407 (2009). G. Catarina, B. Amorim, E. V. Castro, J. M. V. P. Lopes,and N. M. R. Peres, Handbook of Graphene, Vol. 3 (JohnWiley & Sons, 2019) Chap. 6, pp. 177 – 231. C. W. J. Beenakker, Rev. Mod. Phys. , 1337 (2008). G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly,and J. van den Brink, Physical Review B , 073103 (2007). S. Y. Zhou, G.-H. Gweon, A. V. Fedorov, P. N. First,W. A. D. Heer, D.-H. Lee, F. Guinea, A. H. C. Neto, andA. Lanzarai, nature materials , 770 (2007). D. Griffiths, Quantum Mechanics (New Jersey: Pearson,2005). D. Tong, The Quantum Hall Effect (Department of Ap-plied Mathematics and Theoretical Physics, Centre forMathematical Sciences, Wilberforce Road, Cambridge, CB3OBA, UK, 2016). M. Eich, R. Pisoni, H. Overweg, A. Kurzmann, Y. Lee,P. Rickhaus, T. Ihn, K. Ensslin, F. Herman, M. Sigrist,K. Watanabe, and T. Taniguchi, Physical Review X ,031023 (2018). T. Wenz, J. Klochan, F. Hohls, T. Gerster, V. Kashcheyevs,and H. W. Schumacher, Physical Review B , 201409(R)(2019). N. Ubbelohde, F. Hohls, V. Kashcheyevs, T. Wagner,L. Fricke, B. K¨astner, K. Pierz, H. W. Schumacher, andR. J. Haug, Nature Nanotechnology , 46 – (2015). D. Bera, L. Qian, T.-K. Tseng, and P. H. Hol-loway, Quantum Dots and Their Multimodal Applications(Materials, 2010) pp. 2260 – 2345. C. Volk, S. Fringes, B. Terr´es, J. Dauber, S. Engels, S. Trel-lenkamp, and C. Stampfer, Nano Letters , 3581 – (2011). A. Kurzmann, M. Eich, H. Overweg, M. Mangold, F. Her-man, P. Rickhaus, R. Pisoni, Y. Lee, R. Garreis, C. Tong,K. Watanabe, T. Taniguchi, K. Ensslin, and T. Ihn, Phys-ical Review Letters , 1 (2019). J. Bucko, F. Herman, F. Sch¨afer, A. Kurzmann, C. Tong,T. Ihn, M. Sigrist, S. D. Huber, and E. Greplova, To bepublished (2020). S. J. Ahn, P. Moon, T.-H. Kim, H.-W. Kim, H.-C. Shin,E. H. Kim, H. W. Cha, S.-J. Kahng, P. Kim, M. Koshino,Y.-W. Son, C.-W. Yang, and J. R. Ahn, Science , 782(2018). S. Pezzini, V. Miˇseikis, G. Piccinini, S. Forti, S. Pace,R. Engelke, F. Rossella, K. Watanabe, T. Taniguchi, P. Kim,and C. Coletti, Nano Letters , 3313 (2020).28