Analytical and numerical study of uncorrelated disorder on a honeycomb lattice
Kean Loon Lee, Benoît Grémaud, Christian Miniatura, Dominique Delande
aa r X i v : . [ c ond - m a t . d i s - nn ] N ov Analytical and numerical study of uncorrelated disorder on a honeycomb lattice
Kean Loon Lee, Benoˆıt Gr´emaud,
2, 1, 3
Christian Miniatura,
4, 1, 3 and Dominique Delande Centre for Quantum Technologies, National University of Singapore, 3 Science Drive 2, Singapore 117543, Singapore Laboratoire Kastler-Brossel, UPMC-Paris 6, ENS, CNRS; 4 Place Jussieu, F-75005 Paris, France Department of Physics, National University of Singapore, 2 Science Drive 3, Singapore 117542, Singapore INLN, Universit´e de Nice-Sophia Antipolis, CNRS; 1361 route des Lucioles, 06560 Valbonne, France
We consider a tight-binding model on the regular honeycomb lattice with uncorrelated on-sitedisorder. We use two independent methods (recursive Green’s function and self-consistent Bornapproximation) to extract the scattering mean free path, the scattering mean free time, the densityof states and the localization length as a function of the disorder strength. The two methodsgive excellent quantitative agreement for these single-particle properties. Furthermore, a finite-sizescaling analysis reveals that all localization lengths for different lattice sizes and different energies(including the energy at the Dirac points) collapse onto a single curve, in agreement with the one-parameter scaling theory of localization. The predictions of the self-consistent theory of localizationhowever fail to quantitatively reproduce these numerically-extracted localization lengths.
I. INTRODUCTION
Anderson localization (AL) of waves [1] in disor-dered media is a ubiquitous phenomenon which has beenobserved both for classical and quantum waves, e.g.light [2, 3], acoustics [4, 5], water waves [6], ultracoldatoms [7, 8], polaritons [9] and quantum Hall system [10].The scaling theory of localization [11] predicts that athree-dimensional (3D) system exhibits a metal-insulatortransition while 1D and 2D systems always display local-ization at any finite disorder strength. Approximate an-alytical expressions for the localization length in termsof the transport mean free path can be derived withinthe framework of the self-consistent theory of localiza-tion [12, 13]. Dimension two is in fact the critical di-mension for AL and symmetry considerations can playan important role. Indeed, while localization is expectedto take place for spinless time-reversal invariant systems(albeit with an exponentially large localization length),perturbative renormalization group studies on non-linear σ -models suggest that a metal-insulator transition mayoccur in 2D if chirality is present [14, 15]. Such a dis-ordered system with different chiral classes could be re-alized with the honeycomb lattice. The successful iso-lation of graphene flakes in 2004 [16], and the discoverythat graphene samples exhibit a finite electronic conduc-tivity at half-filling although the density of states (DoS)vanishes [17, 18], has thus spurred interest in studyingelectronic transport in graphene in the presence of disor-der, see [19] and references therein.However, even though graphene is a readily-availablephysical realization of a honeycomb lattice, its propertiesare invariably affected by the combined effects of interac-tion, disorder and phonons. The controlled study of dis-order alone in graphene sheets is thus difficult, notwith-standing the fact that engineering disorder with givenstatistical properties seems out of reach. In that respect,ultracold atoms loaded on a graphene-like optical lat-tice [20–22] offer an alternative route and have alreadyproven their key impact in weak and strong localization studies [7, 23–27]. Furthermore, key transport quanti-ties, like the scattering and transport mean free paths orthe localization length, have already been analyzed forspeckle optical potentials in the Born approximation inRef. [13] while engineering disorder with different corre-lation properties is possible.The aim of this paper is to study transport in a dis-ordered honeycomb lattice. We first stick to the simplercase of the tight-binding model for the regular honey-comb lattice in the presence of uncorrelated on-site disor-der characterized by a symmetric box distribution. Themore interesting (but more complicated) case of corre-lated on-site disorder will be the scope of a forthcomingpublication. The novelty of our work lies in the gen-eralization of two known methods, the recursive Green’sfunction method and the self-consistent Born approxima-tion, (i) to extract single-particle properties, such as thescattering mean free time τ , the scattering mean free path ℓ , the density of state (DoS) ν , and the localization length ξ , and (ii) to relate them to experimentally-controlledparameters, such as the tunneling amplitude J and thedisorder strength W . The numerical data obtained fromthese two methods show remarkable agreement and givean accurate estimation of these single-particle properties.Our results further confirm the one-parameter scalinghypothesis[11] for localization but also reveal a quantita-tive disparity with the predictions of the self-consistenttheory of localization [12, 13]. This disparity is yet to beunderstood.Currently, there are three pieces of numerical worksthat are technically relevant to ours. Using a transfer ma-trix technique, Schreiber and Ottomeier [28] have shownthat the localization lengths for various lattices (square,honeycomb and triangular) at the energy band centre E = 0 – thus including the Dirac point of the honeycomblattice – and for various disorder strengths obey the samescaling laws. Using the same method, Xiong and Xiongconcluded that all states are localized but found that thescaling behavior at the charge neutrality point is differentfrom the one at different energies [29]. On the other hand,Lherbier et al. considered the time dynamics of a ran-dom phase wave packet using a real space order- N Kubomethod [30]. They subsequently extracted the diffusionconstant, and hence the scattering mean free path (whichcoincides with the transport mean free path as scatteringby our δ -correlated potential is isotropic), from the timeevolution of the spatial spread of the wave packet. Thisextracted scattering mean free path was then used todeduce the semi-classical conductivity through Einstein-Kubo relation [31] and the localization length throughthe self-consistent theory of localization, a procedure thatour numerical data question.The paper is organized as follows. Sec. II gives theessential ingredients of our model and the eigenstruc-ture of the disorder-free honeycomb lattice in the tight-binding regime. In Sec. III, we introduce the self-energyand detail the Born and self-consistent Born approxima-tions (SCBA). We further derive analytical expressionsfor the self-energy at some particular energies in theweak-disorder limit. In Sec. IV, we introduce the recur-sive Green’s function (RGF) method which is designedto compute exact matrix elements of the Green’s func-tion for large system sizes, with the caveat that actualcomputations can take months on a computer cluster. Afaster but more restricted variant is the recursive trans-fer matrix method. Together with a finite-size analysis,these methods allow us to extract the localization lengthof the disordered honeycomb lattice at any given energy.In Sec. V, we perform a finite-size analysis of the localiza-tion lengths. We show that they can be simultaneouslyscaled for all energies, including the charge neutralitypoint. Furthermore, our results indicate that this univer-sal curve is valid for all lattice types, all energies withinthe energy band and possibly for all types of uncorrelateddisorders, which is not surprising from the viewpoint ofthe one-parameter scaling hypothesis. In addition to thelocalization length, we also extract the scattering meanfree path and the DoS from the recursive Green’s functionmethod evaluated at complex energies. The extractedquantities show remarkable agreement with our resultsfrom the self-consistent Born approximation. The com-parison of the numerically computed localization lengthto the prediction of the self-consistent theory of local-ization shows a fair qualitative agreement, but markedquantitative differences. We finally conclude in Sec. VI.Additional details are given in the Appendices. II. MODEL
We consider here a tight-binding Hamilton operatoracting on a regular honeycomb lattice with on-site disor-der H = H + V = − J X h i,j i (cid:0) | i ih j | + | j ih i | (cid:1) + X i ε i | i ih i | , (1)where | i i refers to a Wannier state localized on site i and h i, j i denotes a sum over nearest neighbors. The hoppingparameter J is usually positive and it fixes the energy B .......................................................................................................................................................... ........................................................................................................................................................................................................................................................... ........................................................ A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... ........................................................ A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... ........................................................ A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... ........................................................ A ................................................................................................................................................................................................................................................................................................................................................................................................................. A ......................................................................................................................................................... B ......................................................................................................................................................... ................................................................................................................................... ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... ................................................................................................................................... ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... ................................................................................................................................... ........................................................................................................................................................................................................................................................................... . ..................................................................................................................................... B .......................................................................................................................................................... ........................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... A ....................................................................................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... ......................................................... ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... ......................................................... ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... ......................................................... ........................................................................................................................................................................................................................................................................... . ................................................................................................................................................................................................................................................................................................ ................................................................................................................................................................................................................................................................................................................................................................................................................................................................................. a c c c ✻❄ ✻ ✻ ✻ FIG. 1. A honeycomb lattice with lattice constant a andits diamond-shaped two-point basis cell (dashed line). Thevectors c l ( l = 1 , ,
3) connect an a -site to its three b -sitenearest-neighbors. scale. Throughout the paper, we assume the diagonal el-ements ε i to be independent random variables character-ized by the same symmetric box probability distribution P ( ε i ) = (cid:26) W for | ε i | ≤ W , . (2)where W is the disorder strength. The disorder has thuszero mean average ε i = 0 and its two-point correlator isthen C ij = ε i ε j = W δ ij , where the bar denotes averag-ing over disorder configurations. The disorder being spa-tially δ -correlated, scattering is isotropic; the scatteringand transport mean free paths are consequently equal.We now briefly review the eigenstructure of thedisorder-free Hamiltonian H . We refer to Ref. [20] formore details. The regular honeycomb lattice being a tri-angular Bravais lattice with a two-site basis cell, it can bepictured as two shifted triangular sublattices denoted by a and b , see Fig. 1. As a consequence, the coordinationnumber of the honeycomb lattice is three. We denote by c l ( l = 1 , ,
3) the link vectors connecting a site i ∈ a toits three nearest-neighbors j l ∈ B ( | c l | = a , a being thelattice constant). We next define the structure factor ofthe honeycomb lattice for nearest-neighbor hopping as f ( k ) = | f ( k ) | e i ϕ ( k ) = X l =1 e i k · c l . (3)For the honeycomb lattice depicted in Fig. 1, where c points along the y -axis, we have | f ( k ) | = 1+4 cos √ k x a ! +4 cos √ k x a ! cos (cid:18) k y a (cid:19) . (4) π/ (3 √
3) 4 π/ √ k x a − π/ π/ k y a K ′ K FIG. 2. The first Brillouin zone B (in grey) of the honey-comb lattice shown in Fig. 1 and the associated reciprocal lat-tice vectors (arrows). The two empty circles mark the Diracpoints, around which the energy dispersion is linear. By usingconvenient reciprocal lattice translations, the first Brillouinzone can be mapped onto the shown rectangle which is thenused as the integration space in all integrals in the paper. As shown in Ref. [20], the eigenstructure of H is de-fined by H = X k ∈B ,s = ± ε k s | k s ih k s | , (5) ε k s = s ε k = s J | f ( k ) | , (6) | k s i = 1 √ N c (cid:2) X i ∈ A e i k · r i | i i − s e − i ϕ ( k ) X j ∈ B e i k · r j | j i (cid:3) (7)where B is the first Brillouin zone of the honeycomb lat-tice (see Fig. 2), s is the band index ( s = +1 for theupper, or conduction, band; s = − N c is the number of Bravais cells of thelattice.The full spectrum span the energy interval [ − J, J ]and band crossing can only occur at zero energy since ε + ( k ) = ε − ( k ) implies that f ( k ) = 0. This happens atthe Dirac points K and K ′ = − K where K = (cid:18) π √ a , (cid:19) (8)for the honeycomb lattice in Fig. 1. In the solid-statecommunity, the energy at the Dirac points is usually re-ferred to as the charge neutrality (energy) point becausethe number of energy states above and below this pointare equal. As a consequence, when the gate voltage isfixed at ε ± ( K ) = 0 in a graphene sample, the graphenesheet is charge neutral since the particle and hole statesare balanced.In the rest of the paper, we will be mainly interestedin the different localization properties of our model forenergies E lying near the band edges, E ≈ ± J , andnear the band centre, E ≈
0. In the first case, k liesnear the centre of the Brillouin zone ( ka ≪
1, where k = | k | ), while in the second case k lies near the Diracpoints ( qa ≪
1, where q = k − K , and similarly around K ′ ).Near the band edges E = ± J , the dispersion relation is quadratic, thus representing free massive particles ε k s ≈ Js (cid:18) − a k (cid:19) ka ≪ , (9)the effective mass m s = − s ~ / (3 Ja ) being negativefor the upper band and positive for the lower band.Near the charge neutrality point E = 0, the disper-sion relation is linear, thus representing the celebrated“relativistic” massless Dirac particles propagating witha velocity c playing the role of an effective speed of light, ε k s ≈ s J aq = s ~ cq ( qa ≪ , (10) c = 3 aJ ~ . (11) III. SELF-ENERGY
An initial Bloch state | k s i propagating in the latticewill suffer scattering by the disorder fluctuations and thuswill be depleted, decaying exponentially over time witha time constant which is the scattering mean free time τ k s . This coherent propagation and decay are describedby the disorder-averaged Green’s function G which obeythe Dyson’s equation G ( z ) = G ( z ) + G ( z )Σ( z ) G ( z ) , (12) z being a point in the complex energy plane. Within ourmodel, the disorder-free Green’s function is given by G ( z ) = 1 z − H = X k s | k s ih k s | z − ε k s . (13)The Dyson equation features the self-energy Σ which isa central quantity given by a perturbative sum of irre-ducible diagrams [32]. As the disorder average restoresthe lattice translation symmetries and characteristics,one has h k ′ s ′ | G ( z ) | k s i = G k s ( z ) δ kk ′ δ ss ′ , (14) h k ′ s ′ | Σ( z ) | k s i = Σ k s ( z ) δ kk ′ δ ss ′ , (15)where the dependence of Σ k s ( z ) on k is usually smooth.Due to particle-hole symmetry, the spectrum of H issymmetric with respect to E = 0. Since the on-site energies are themselves independent symmetrically-distributed random variables, it is easy to show that h k s | G ( z ) | k s i ∗ = −h k , − s | G ( − z ∗ ) | k , − s i , where the starstands for complex conjugation. In turn, the self-energysatisfies Σ k ,s ( z ) = Σ ∗ k , − s ( − z ∗ ) . (16) These identities no longer hold for a speckle potential since itbreaks the V → − V symmetry. Furthermore, because of time-reversal invariance, theself-energy is the same at ± k . This means that it is suf-ficient to study the negative energy sector and forwardpropagation ( k x ≥ ~ τ k s ( E ) = − ImΣ k s ( E ) , (17)is independent of the band index. A. Born approximation
An analytical expression for Σ k s ( z ) is generally notavailable and one has to resort to approximations to findthe self-energy. For weak disorder, the simplest approxi-mation is the Born approximation which consists in dis-carding all terms of the full diagrammatic perturbativeexpansion in Fig. 3 except the first one. Its lattice matrixelements are h i | Σ Born ( z ) | j i = C ij h i | G ( z ) | j i = W I ( z ) δ ij (18)with I ( z ) = h i | G ( z ) | i i = Z B d k Ω zz − J | f ( k ) | , (19)where Ω = 8 π / (3 √ a ) is the area of the first Brillouinzone. For uncorrelated on-site disorder, the self-energyat the Born approximation is a scalar:Σ k s ( z ) ≈ Σ Born ( z ) = W I ( z ) . (20)The average Green’s function then reads G Born ( z ) = 1 z − H − Σ Born ( z ) = G ( z − Σ Born ( z )) . (21)The expression of I ( z ) in terms of elliptic integrals [33,34] is given in Appendix A. B. Self-consistent Born approximation (SCBA)
This approximation scheme builds on the Born approx-imation by replacing G by G in Eq. (18). It is morepowerful as it amounts to sum the infinite subclass of“rainbow” diagrams given in Fig. 3. It gives the follow-ing self-consistent equation:Σ SCBA ( z ) = W h i | G ( z ) | i i = W I ( z − Σ SCBA ( z )) , (22)which is easy to solve numerically.In the following we will use the parametrizationΣ SCBA ( z ) = γJe − i θ with γ positive. For the scatter-ing time to be positive, we must have ImΣ k s <
0, which = + + + ... + + + ...
FIG. 3. The “rainbow” subclass of diagrams retained tocompute the self-energy in the self-consistent Born approx-imation. The double line with arrow denotes an averagedGreen’s function G while a single line with arrow denotes adisorder-free lattice Green’s function G . Two vertices (soliddots) connected by a dashed line represent a 2-point correlator C ij = ε i ε j . The Born approximation consists in computingthe self-energy with the first “rainbow” diagram only. For un-correlated disorder, the connected vertices correspond thus tothe same site. In this case the self-energy is a scalar operatordepending only on the energy and the disorder strength. enforces 0 ≤ θ ≤ π . Eq. (16) then translate into γ ( − E, W ) = γ ( E, W ) (23) θ ( − E, W ) = π − θ ( E, W ) . (24)This implies θ (0 , W ) = π/ W .Figs. 4 and 5 show the E -dependence of θ and γ in theSCBA for some particular disorder strengths W , whileFigs. 6 and 7 show the W -dependence of θ and γ in theSCBA at some particular energies E . In the followingsubsection, we investigate SCBA analytically in the weakdisorder regime. C. Weak disorder limit
In the weak disorder regime W ≪ J , one expects γ ≪
1. Several analytical results can then be derived in thislimit. Some details are exposed in Appendix B.
1. Charge neutrality point
At the charge neutrality point, we have z = E = 0and the function I in (22) is thus evaluated at the di-mensionless complex number Z = − Σ SCBA (0 , W ) /J .Furthermore, because of Eq. (16), θ (0 , W ) = π/
2, andΣ
SCBA (0 , W ) = − i γ (0 , W ) J is purely imaginary. Equa-tions (22) and (A3) have two solutions. One is the trivialsolution γ (0 , W ) = 0 while the nontrivial solution solves Z B d k Ω 1 | f ( k ) | + γ (0 , W ) = 12 J W . (25) . . . . . . . E/J . . . . . . . . . θ / π BA W = 0 . JW = JW = 5 J FIG. 4. (Color online) Using the parametrizationΣ( z ) = γJ e − i θ ( γ ≥
0, 0 ≤ θ ≤ π ), the figure shows theangle θ as a function of the energy E in the Born approxima-tion (black open circles) and in the SCBA for various disorderstrengths. In the Born approximation, θ is independent of W .At W = 0 . J , SCBA and the Born approximation are essen-tially identical. . . . . . . . E/J . . . . . . . . . γ γ, W = 0 . J γ, W = Jγ, W = 5 J FIG. 5. (Color online) Using the parametrization Σ( z ) = γJ e − i θ ( γ ≥
0, 0 ≤ θ ≤ π ), the figure shows the amplitude γ as a function of the energy E in the SCBA for various disorderstrengths. As W →
0, the peak at E = J develops into thevan Hove singularity. An expansion of the elliptic integrals appearing in thefirst line of Eq. (A3) for | Z | ≪ γ (0 , W ) ≈ (cid:0) − π √ J W (cid:1) . (26)Shon et al. found essentially the same type of result γJ = ε c exp( − A /W ), see Eq. (3.21) in [35]. The differ-ence between their parameter values and ours arises fromtheir introduction of the cut-off energy ε c . In our treat-ment, we use the exact dispersion relation, beyond thelinear approximation around the Dirac points, making W/J . . . . . . . . . . θ / π E = 0 . JE = JE = 3 J FIG. 6. (Color online) Using the parametrization Σ( z ) = γJ e − i θ ( γ ≥
0, 0 ≤ θ ≤ π ), the figure shows the angle θ asa function of the disorder strength W in the SCBA near thecharge neutrality point, at the van Hove singularity and atthe band edge. W/J . . . . . . γ E = 0 . JE = JE = 3 J FIG. 7. (Color online) Using the parametrization Σ( z ) = γJ e − i θ ( γ ≥
0, 0 ≤ θ ≤ π ), the figure shows the amplitude γ as a function of the disorder strength W in the SCBA nearthe charge neutrality point, at the van Hove singularity andat the band edge. the introduction of an artificial cut-off unnecessary.
2. Van Hove singularities
At the van Hove singularities, the disorder-free DoSdiverges. For our honeycomb tight-binding model, thisoccurs at z = E = ± J [20] and the function I in (22)is now evaluated at the dimensionless complex number Z + = 1 − Σ SCBA ( J, W ) /J . A small-parameter expansionof equations (22) and (A3) around Z = 1 now leads to γ ( J, W ) ≈ W πJ (cid:18) ln (cid:16) πJ W (cid:17) + ln ln (cid:16) πJ W (cid:17)(cid:19) , (27) θ ( J, W ) ≈ π (cid:0) ln( πJ W ) − ln ln( πJ W ) (cid:1) ! . (28)One sees that, at lowest order, the self-energy at thevan Hove singularities is purely imaginary too, θ ≈ π/
3. Band edges
At the disorder-free band edges, z = E = ± J . By thesame token, we evaluate the function I in (22) at the di-mensionless complex number Z + = 3 − Σ SCBA (3 J, W ) /J .A small-parameter expansion around Z = 3 now leads toΣ SCBA (3 J, W ) ≈ √ W πJ (cid:16) ln (cid:16) √ πJ W (cid:17) − i π − ln(ln (cid:16) √ πJ W (cid:17) − i π ) (cid:17) . (29)Anticipating results displayed and discussed in Para-graph V C, within the SCBA scheme, the DoS van-ishes outside a finite energy band, with a square-rootbehavior near the band edge. This SCBA band edgeis approximately given by the solution of the equation E − ReΣ
SCBA ( E ) = 3 J and does not coincide with theexact band edge 3 J + W/ J and 3 J + W/ IV. RECURSIVE GREEN’S FUNCTIONMETHOD
The RGF method is based upon the division of the sys-tem in smaller sections for which the Green’s functions can be calculated more easily [36]. These sections arethen “glued together” one after one and Dyson’s equa-tion [37] is repeatedly used to derive the full Green’s func-tion in terms of the Green’s functions of the smaller sec-tions. To study localization, we will use a generalizedversion [38, 39] of the RGF method. This generalizedversion enables us to extract any lattice matrix elementof the Green’s function conveniently and with high nu-merical stability.Applying the RGF scheme to our case amounts to con-sider a finite quasi-1D lattice strip and to divide it into N vertical slices, the two open ends being along the horizon-tal direction, see Fig. 8. Denoting by H N − the nearest-neighbor tight-binding Hamilton operator for a strip with( N −
1) slices, the Hamilton operator H N obtained by glu-ing an additional slice N can be split into three terms, H N = H N − + H slice N + H hop N − ,N + H hop N,N − , (30)where H hop N − ,N is the nearest-neighbor hop operator con-necting sites within slice ( N −
1) to sites within slice N (and vice-versa for H hop N,N − ). Since no external gaugefields are present in our model, we safely consider the hopoperators to be real from now on. H slice N is the nearest-neighbor tight-binding Hamilton operator for the isolatedslice N before it is stacked to the others. It thus includesthe on-site disorder diagonal term V in Eq. (1).Using Dyson’s equation, the submatrix G ( N ) l,n of thefull retarded Green’s function G ( N ) = ( E + i0 + − H N ) − coupling slice l to slice n at energy E can be obtainedthrough the following recursion relations, G ( N ) l,n = G ( N − l,n + G ( N − l,N − H hop N − ,N G ( N ) N,n , (31a) G ( N ) l,N = G ( N − l,N − H hop N − ,N G ( N ) N,N , (31b) G ( N ) N,n = G ( N ) N,N H hop N,N − G ( N − N − ,n , (31c) G ( N ) N,N = (cid:18) E + i0 + − H slice N − H hop N,N − G ( N − N − ,N − H hop N − ,N (cid:19) − (31d)with 1 ≤ l ≤ n ≤ ( N − FIG. 8. Schematic illustration of the RGF method with astrip of square lattice of length L N and width L M . The systemunder consideration is constructed by repeated stacking ofvertical slices of length L M . The Green’s function at eachstep of the stacking process is calculated recursively usingDyson’s equation. Knowing the Green’s function G ( N − =( E +i0 + − H N − ) − for the system with ( N −
1) slices, the hopoperator H hop N = H hop N − ,N + H hop N,N − coupling slices ( N − N , and the Hamilton operator H slice N of the isolated slice N , one can exactly compute the Green’s function G ( N ) forthe whole system of N slices, see Eqs. (31). chair geometries. Once the stacking is properly defined,each lattice site i ≡ ( n, m ) will then be labelled by twointegers. The first one 1 ≤ n ≤ N corresponds to theslice it belongs to from left to right, N being the totalnumber of such slices. The second one 1 ≤ m ≤ M corre-sponds to its position along the slice from bottom to top, M being the total number of sites per slice. For bothgeometries we have checked that the matrix elements ofthe Green’s function obtained through the recursive algo-rithm agree well with those computed by direct inversionof the operator ( E + i0 + − H N ). A. Zigzag configuration
We first consider the zigzag (ZZ) configuration de-picted in Fig. 9, where L vertical zigzag chains of thehoneycomb lattice are stacked along the horizontal di-rection. In the following we will be essentially interestedin the case of periodic boundary conditions in the verticaldirection and open boundary conditions in the horizon-tal one. This imposes the number of sites in each zigzagchain to be an even integer 2 M . For each zigzag chain,one can define a -type (resp. b -type) vertical slices con-taining only a -sites (resp. b -sites). Each of these verticalslices contain M sites and there are N = 2 L such slices, a -type slices alternating with b -type slices. Lattice sitesare thus parametrized by ( n, m ) with 1 ≤ n ≤ N = 2 L and 1 ≤ m ≤ M , see Fig. 9. The width and length of thishoneycomb ZZ strip are L M = √ M a and L N ≈ N a/ N ≫ H slice n asso-ciated to the isolated slice n is simply a M × M diagonal A ......................................................................................................................................................... B ......................................................................................................................................................... . .................................................................................................................................. ...................................................................................................................................... ..................................................................................................................................... ...................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... . .................................................................................................................................. ...................................................................................................................................... ..................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... . .................................................................................................................................. ...................................................................................................................................... ..................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... . .................................................................................................................................. ...................................................................................................................................... ..................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... . .................................................................................................................................. ...................................................................................................................................... ..................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... . .................................................................................................................................. ...................................................................................................................................... ..................................................................................................................................... ...................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... . .................................................................................................................................. ...................................................................................................................................... ..................................................................................................................................... ...................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... . .................................................................................................................................. ...................................................................................................................................... ..................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... . .................................................................................................................................. ...................................................................................................................................... ..................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... . .................................................................................................................................. ...................................................................................................................................... ..................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... . .................................................................................................................................. ...................................................................................................................................... ..................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... . .................................................................................................................................. ...................................................................................................................................... ..................................................................................................................................... ............................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................ ...................................................................................................................................................................................................................................................................................................................................................................................................................................................................................... a ✲✛ FIG. 9. In the honeycomb zigzag configuration with periodicboundary conditions along the vertical direction, L zigzagvertical chains, each containing 2 M sites, are stacked alongthe horizontal direction ( L = 4 and M = 3 in the figure).Each zigzag chain defines two vertical slices (dashed verticallines). These slices are labeled from left to right by the integer1 ≤ n ≤ N = 2 L . Each slice contains M sites labeled frombottom to top by the integer 1 ≤ m ≤ M . The slice genderalternates between a -type and b -type. Each site in the latticeis uniquely parametrized by the couple of integers ( n, m ). matrix with entries given by the M on-site disorder ele-ments { ε n,m } . Furthermore, one can easily see that H hop n,n +1 = − J if n ≡ . H hop n,n +1 = − J D if n ≡ . H hop n,n +1 = − J if n ≡ . H hop n,n +1 = − J D T if n ≡ . H hop n +1 ,n = ( H hop n,n +1 ) T (32)where the T-superscript means matrix transposition and D is the M × M matrix given by D = · · · p · · · · · · · · · . (33)When periodic boundary conditions in the vertical di-rection are used (which is our case), p = 1. For openboundary conditions in the vertical direction, one wouldhave p = 0.An astute reader may have noticed that the determi-nant of D vanishes for M even and periodic boundaryconditions. In this case the transfer matrix method [36]cannot be implemented as it would require the inversionof D or D T . To avoid this pitfall, one can neverthe-less always choose M odd. Note however that the RGFscheme is perfectly immune to this breakdown and itsimplementation does not suffer any flaw as we have dulychecked. B. Armchair configuration A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... ........................................................ A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... ........................................................ A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... ........................................................ A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... ........................................................ A ......................................................................................................................................................... B ......................................................................................................................................................... ................................................................................................................................... ..................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... ................................................................................................................................... ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... ................................................................................................................................... ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... ................................................................................................................................... ........................................................................................................................................................................................................................................................................... . ..................................................................................................................................... A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ........................................................................................................................................................................................................................................................................................... . ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... ......................................................... ..................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... ......................................................... ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... ......................................................... ........................................................................................................................................................................................................................................................................... A ......................................................................................................................................................... B ......................................................................................................................................................... ......................................................... ........................................................................................................................................................................................................................................................................... . ............................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................... .......................................................................................................................................................................................................................................................................................................................................................................................................................................................................... a ✻❄ FIG. 10. In the armchair configuration (AC), M horizontalzigzag chains, each containing N sites, are stacked vertically( M = 4 and N = 8 in the figure). There are thus N verticalslices labeled from left to right by the integer 1 ≤ n ≤ N (dashed lines). Each vertical slice contains M sites labeledfrom bottom to top by the integer 1 ≤ m ≤ M . If periodicboundary conditions are imposed along the vertical direction, M must be even. The site gender within a slice alternatesbetween a -type and b -type. Each site in the lattice is uniquelyparametrized by the couple of integers ( n, m ). We now turn to the honeycomb armchair configura-tion where M horizontal zigzag chains, each containing N sites, are stacked along the vertical direction, see Fig. 10.When periodic boundary conditions are imposed alongthe vertical direction, M must be even. The left ver-tical boundary of the lattice is now reminiscent of theshape of an armchair. Using the same recipe, we slicethe lattice with vertical lines. There are now N suchvertical slices, each containing M sites. The width andlength of this honeycomb AC strip are L M = 3 M a/ L N ≈ √ N a/ N ≫ M × M hop matrices satisfy H hop n,n +1 = H hop n +1 ,n = − J , while H slice n = − J X n + ε n where ε n is a M × M diagonal matrix with entries equal to the M on-sitedisorder { ε n,m } in slice n . X n is a M × M sparse matrixthat couples each site to its nearest neighbors within slice n , namely X n = · · · · · · · · · · · · · · · · · · if n is odd, · · · p · · · · · · · · · · · · · · · · · · p · · · if n is even,(34)where p = 1 for periodic boundary conditions along thevertical direction (and M even) and p = 0 for open ones. V. NUMERICAL RESULTSA. Localization length ξ Localized wave functions are expected to decrease ex-ponentially at large distances. We thus compute the lo-calization length along the horizontal direction as λ M = − lim N →∞ L N ln Tr (cid:16) J G ( N )1 ,N [ G ( N )1 ,N ] † (cid:17) , (35)where the bar indicates the average over disorder config-urations. Note that it corresponds to the log-averagedtransmission in quasi-1D systems which has the niceproperty of being additive and self-averaging when newslices are added [40].This finite-size localization length λ M depends on theenergy E and the disorder strength W , but also on thelattice configuration (ZZ or AC) and width L M . In allour simulations, we have used a sufficiently large numberof randomly-generated disorder configurations for eachset of parameters ( E , W , M and lattice configuration)such that the estimated relative error in computing 1 /λ M is less than 0 . /λ M is always com-puted with lengths L N greater than 5 λ M . This meansthat larger number of samples are required as λ M in-creases. To avoid numerical underflow, a rescaling of G ( N )1 ,N is done through Eq. (31b) every 10 multiplications Some authors define the localization length λ M through ln | ψ | (as we do here) or through ln | ψ | . There is a factor of 2 betweenthese two possible definitions. approximately. RGF equations require one matrix in-version per slice. To speed up the computation for theAC lattice we switched to the recursive transfer matrixmethod [36] where a true matrix inversion is carried outonly after the transfer matrix equation is applied for 10slices.Based on ideas of the renormalization group [36] andthe scaling theory of localization [11], it was conjecturedthat all data points in a λ M /L M versus ξ/L M plot shouldcollapse onto the same universal curve, λ M L M = F (cid:18) ξL M (cid:19) (36)where the infinite-lattice localization length ξ dependsonly on energy E and disorder strength W . Our resultsin Fig. 11 fully supports this conjecture for sufficientlylarge L M .When L M ≫ ξ , the system becomes insensitive to thevertical boundary conditions and one expects λ M ≈ ξ .The scaling function should thus satisfy F ( x ) ≈ x forsmall x . We chose the honeycomb AC configuration at E = 3 J and W = 10 J as the reference data set forthis limiting behavior, therefore fixing ξ = (2 . ± . a in this case, see black data at the bottom left cornerin Fig. 11. At finite system width, we nevertheless ex-pect the computed localization length λ M to be slightlyshorter than the true infinite-lattice localization length ξ . This behavior is accounted for by a Taylor expansionof the scaling function F ( x ) = x − αx + O ( x ) . Fromour data we infer α = (1 . ± . x , the full scalingfunction F ( x ) is then built by finding, for each of our datasets, the corresponding ξ allowing to patch them on asingle smooth curve. This procedure has been done withour AC data sets obtained in the range E ∈ [0 . J, J ] and W ∈ [1 . J, J ]. We have also checked that AC data setsobtained for E slightly larger than 3 J and various W , aswell as ZZ data sets obtained at different energies ( E =0 . J, J, . J ) and W = 1 . J , all collapse onto this verysame curve too. In particular, we numerically confirmthat at E = 0 . J , E = 2 . J , and W = 1 . J , we get thesame ξ for the ZZ and AC configurations. This might beunderstood as a consequence of the isotropic dispersionrelation in the two energy ranges, see Eqs. (9) and (11).On the other hand, at E = J , the extracted ξ are slightlydifferent for the two honeycomb configurations.For the AC configuration, an important finding is thatdata sets obtained at the charge neutrality point ( E = 0)for W ≥ . J can all be scaled by F ( x ). Fig. 12 shows thecollapsed data for W = 2 . J , 4 J and 6 J ( ξ = 500 a, . a and 16 . a respectively). Performing similar calculationsfor different lattice models (honeycomb, square and tri-angular), Schreiber and Ottomeier also found that theirdata obtained at E = 0 and W ≥ J could be col-lapsed onto a single universal curve [28]. Since our dataat W = 4 J and 6 J with M = 46 agree with theirs at M = 40 (after proper unit conversion), we infer that weindeed found the same universal curve. While it is difficult to numerically confirm that data at E = 0 and W < . J can be scaled by F ( x ) (see analysisbelow for the reason of this difficulty), there is no appar-ent reason why they should not. Therefore, our results,combined with those in Ref. [28], imply that all curves λ M /L M as a function of 1 /L M can be collapsed onto thesame universal curve by an appropriate one-parameterscaling, independently of the lattice type and disorderstrength, at least for uncorrelated box-distributed on-sitedisorder and energies within (and slightly outside of) theenergy band of the clean system. This is in marked con-trast with the findings of Ref. [29], which claim a differentscaling behavior at E = 0 (see below for a possible ex-planation of this apparent contradiction).For the clean system near E = 0, the energy shellslices the band structure near the Dirac points, defin-ing two circles of allowed wavevectors centred on K and − K . The initial state | k s i belongs to one of these cir-cles, say the one around K . Introducing the deviationwavevector q = k − K , the circles have radius | q | . In thepresence of disorder, the two circles are broadened in theenergy shell and form rings of allowed wavevectors withmean radius | q | . Since disorder is uncorrelated in space,the two rings are coupled by scattering. However, thescattering processes being elastic, starting from the ini-tial state q , only these two rings around the Dirac pointswill be populated in the course of time. Loosely speak-ing, if scattering does change the direction of q , it cannotchange its modulus. Around E = 0, the dynamics is thuscharacterized by small wavevectors, i.e. by long distancesin real space. This is why we expect finite-size effects tobe larger around E = 0. This can be seen in Fig. 13where the dependence of λ M on L M shows a damped os-cillatory behavior when L M /a is not large enough. Asexplained below, the very existence of oscillations can betraced back to the opening of new scattering channelswhen L M increases whereas the damping can be tracedback to disorder broadening.Consider the finite-size AC configuration. For pe-riodic boundary conditions, the allowed values for thewavevector k along Ox and Oy are quantized accordingto k x = n x ∆ k N and k y = n y ∆ k M with ∆ k N = 2 π/L N ,∆ k M = 2 π/L M , n x ∈ { , · · · , N } and n y ∈ { , · · · , M } .In the vicinity of E = 0, the curve at constant E willenclose a small circular area around the Dirac point K in the Brillouin zone containing few such discrete points.When the system size, notably L M , is increased, morepoints enter this area, which means that more channelsopen for propagation. The period of the oscillations ob-served in Fig. 13 thus corresponds to the change in L M allowing one new channel to open, i.e. ∆ L M = 2 π/q, where q is the radius of the circle at energy E, . Us-ing Eq. (11), we estimate the period of the oscillationsto be roughly of the order of ∆ L M = 3 πaJ/ | E | . How-ever, each quantized k -mode corresponding to energy E is broadened by disorder, typically by ℓ − . As a conse-quence modes cannot be distinguished from each otherwhen ∆ k M ℓ ≈
1. Near the Dirac point, one has ℓ = cτ FIG. 11. (Color online) Single-parameter scaling law for the AC honeycomb lattice with uncorrelated on-site disorder. Wenumerically compute the finite-size localization length λ M at various energy 0 . J ≤ E ≤ J and disorder strength 1 . J ≤ W ≤ J . The fact that all data sets can be collapsed onto a single log-log curve where ξ is the infinite-lattice localization lengthconfirms the scaling theory of localization [11]. Each disjoint color represents a parameter pair ( E , W ) as shown in the inset. and we find that the oscillations due to the opening ofnew channels are washed out as soon as aL M . − E, W )3 πJ . (37)Using SCBA at W = 1 . J , this simple estimate gives a/L M ≈ × − at E = 0 . J and a/L M ≈ . × − at E = 0 . J , in perfect agreement with the data collectedin Fig. 13.To obtain a reliable estimate of ξ , it is necessary togo beyond the oscillatory region. This is numericallychallenging at E = 0 for weak disorder, see (26). Forexample, the SCBA estimate at W = 1 . J now gives a/L M ≈ × − , way out of the capabilities of the RGFscheme. We believe this is the reason why a differentscaling at E = 0 was claimed in [29]: the sample sizesused by the Authors were probably not large enough toreach the scaling region. As seen in Fig. 13, it is easyto miss the oscillations at E = 0 for small W. The curvethere is clearly different from the other ones and cannotbe collapsed by translation (in log scale) onto the uni-versal curve, F ( x ) . This might be the origin of the claim in Ref. [29] that a different scaling is observed at E = 0 . This is however only a finite-size effect.In Fig. 14, we show the variations of the localizationlength ξ as a function of the energy E for different dis-order strengths W . In the linear dispersion regime (butnot too close to E = 0), ξ appears to be a constantwhose magnitude increases as W decreases. When get-ting closer to E = 0, ξ decreases (see curve at W = 2 . J )and a small dip occurs. A similar feature is predictedin Ref. [30] when computing ξ from the numerically-evaluated semi-classical conductivity. Extrapolating ourresults, we find that ξ is finite at E = 0 but increasessharply when W decreases. Taking graphene as an ex-ample ( a = 0 . , J ≈ . ξ is of the order of 10 a ≈ . W = 1 . J ≈ . E = J . This is surpris-ing as the DoS at the van Hove singularity diverges inthe absence of disorder, which implies more propagationchannels. Finally we see that ξ decreases as E gets closer1 FIG. 12. (Color online) Collapsing AC data obtained at E = 0and disorder strengths W ≥ . J onto the scaling function F ( x ) shown in Fig. 11. For W = 2 . J (circles), 4 J (squares)and 6 J (diamonds), we have ξ = 500 a , 70 . a and 16 . a re-spectively. The largest transverse sizes are given by M = 300,70, and 70 for W = 2 . J , 4 J and 6 J respectively. Data for W = 2 . J display a small oscillation around the universalcurve.FIG. 13. (Color online) Plot of λ M /L M as a function of a/L M at weak disorder W = 1 . J . The oscillations observed atsmall L M results from the opening of new scattering chan-nels when L M is gradually increased. In the linear dispersionregime, the oscillation period in L M is crudely approximatedas inversely proportional to E , with a small correction due todisorder. to the band edge ( E = 3 J ) and slightly extends beyondit.To conclude this subsection, we compare in Fig. 14 ournumerical data near E = 0 and E = 3 J at W = 1 . J withthe analytical prediction of the self-consistent theory oflocalization [12] (see section V B 3). As can be seen, ifthe general trend is qualitatively satisfactory and evensemi-quantitatively correct (within a factor 10) near theband edge, the prediction is off by more than 3 orders ofmagnitude near the band centre E = 0. FIG. 14. (Color online) Infinite-lattice localization length ξ (in units of a ) versus energy E (in units of J ) for variousdisorder strengths W . For each set ( E, W ), ξ is extractedfrom our numerically-computed λ M , Eq. (35), with the helpof the scaling function F ( x ) shown in Fig. 11. The orangesolid line gives ξ = 2 ℓ p e πκℓ − W = 1 . J as inferredfrom the self-consistent theory of localization, see Eqs. (45).In this analytical prediction, the scattering mean free path ℓ is estimated using SCBA and κ is defined in Eq. (45). If thequalitative behavior is relatively satisfactory, the quantitativepredictions are off by several orders of magnitude in the lineardispersion regime and by about an order of magnitude in thequadratic dispersion regime. B. Scattering mean free path ℓ
1. RGF numerical estimate
The elastic scattering mean free path ℓ defines the dis-tance a particle travels on average without being scat-tered. For the negative energy sector (see the reason whyin the following subsection), we compute the 1D averagedretarded wave function asΨ N = JM M X m,m ′ =1 (cid:20) G ( N +2 n ) n,n + N (cid:21) m,m ′ , (38)where the retarded submatrix G ( N +2 n ) n,n + N connecting slices n and ( n + N ) is evaluated at complex energy E + i η ( E ≤
0) and disorder strength W , η being a small positivenumber. The indices m and m ′ label the sites within eachslice respectively. In the following, n is chosen such thatthe distance between slice n and the closest open-endboundary is large enough (see below).As shown in Fig. 15, | Ψ N | falls off exponentially witha decay constant ℓ ( η ). The scattering mean free path isfurther obtained as ℓ = lim η → ℓ ( η ). In our simulations,the number of disorder configurations has been chosensuch that we have 2 × samples near the band edgeat − J (quadratic dispersion regime) and 2 × sam-ples near the band centre (linear dispersion regime). We2 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6-0.4-0.200.20.40.60.81 I m Ψ N Re Ψ N ← − N E = − . J, η = 0 . JW = 1 . J, M = 200 , Maximum N = 30AC FIG. 15. Disorder-averaged wave function ψ N shown in thecomplex plane as the longitudinal size of the lattice, fixed bythe number of slices N , increases (see arrow for the flow).The plot has been made for the honeycomb AC configurationwith M = 200 sites per slice. The complex energy is E +i η = − . J + i0 . J and the disorder strength is W = 1 . J .The spiraling feature shows that Ψ N ∼ exp (cid:0) i k w ( η ) L N (cid:1) × exp (cid:16) − L N ℓ ( η ) (cid:17) , from which we can extract the decay constant ℓ ( η ). The scattering mean free path is further obtained aslim η → ℓ ( η ). have imposed a minimum distance of 5 ℓ ( η ) between slice n and its closest open-end boundary, thereby fixing thevalue of n in Eq. (38). This ensures that the exponen-tial decay is insensitive to the open-end boundaries. Ournumerical results are thus equivalent to those that wouldhave been obtained with an infinite tube. We would liketo mention here an important technical remark. Whensumming the submatrix entries in (38), serious numericalcancellation errors occur when E and J have same signs.In our case J >
E <
0. The reason behind this numerical instabil-ity is the relative sign between the quasi-Bloch sublatticecomponents of | k s i , see Eq. 7.The decay constant ℓ ( η ) is extracted by performinga linear fit on ln | Ψ N | as a function of L N . We see inFig. 16 that a/ℓ ( η ) displays oscillations similar to thoseobserved for ξ when the transverse size (fixed by the num-ber of sites M ) increases. Again, the oscillations can beaccounted for by the opening of new scattering channels.By further increasing M , the oscillations damp and fi-nally the fluctuations in a/ℓ ( η ) reach the level of the er-ror bars themselves. This happens when the width L M roughly exceeds 10 ℓ ( η ). Beyond this point, the systemdoes not feel anymore the transverse periodic boundarycondition and the corresponding ℓ ( η ) are reliable esti-mates of the infinite-lattice case. We show how to extract ℓ from these reliable estimates in Fig. 17. Note that adirect evaluation of ℓ at η = 0 is in principle possible, ↓ a / ℓ ( η ) M E = − . J, η = 0 . JW = 0 . J AC 10 ℓ ( η ) FIG. 16. Oscillation of 1 /ℓ ( η ) as a function of the trans-verse size M . The oscillation period can be accounted forby the presence of new scattering channels. The oscillationstops when the system does not sense the periodic boundarycondition, i.e. L M > ℓ ( η ). a / ℓ ( η ) η/J E = − . J, W = 1 . J FIG. 17. (Color online) Evaluation of the infinite-lattice scat-tering mean free path ℓ at a fixed energy E and disorderstrength W ( E = − . J and W = 1 . J in the figure). Asequence of reliable estimates a/ℓ ( η ) is first obtained by theRGF technique for smaller and smaller η (black open circles).A quadratic fit (solid black line) is then used to interpolatethe value at η = 0, from which ℓ is deduced. For comparison,we give the results obtained with the SCBA method (openred squares). For the chosen disorder strength the agreementis excellent. but is unfortunately affected by huge fluctuations as theGreen’s function is computed on the real axis where itspoles lie. This direct method proves thus much less effi-cient in practice.3
2. SCBA estimate
In the same spirit, an SCBA estimate of ℓ can be ob-tained by using G ( E − Σ( E )) in (38) and by consideringa finite horizontal lattice strip of length L with periodicboundary conditions in the transverse direction in thelimit L M , L → ∞ . In the limit of an infinite AC hon-eycomb lattice, a careful but straightforward calculationusing Eqs. (7) and (13) leads toΨ N = J Z π dq π e i Nq z + J (1 + 2 cos q ) , (39)where z = E − Σ( E ) and q = √ k x a/
2. The self-energyΣ( E ) is computed using Eqs. (22) and (A3). The distancebetween the two slices being L N = √ aN/ ℓ = − lim N →∞ ln | Ψ N | L N . (40)The exponential decay of Ψ N is given by the imaginarypart of the complex pole Q = Q r + i Q i solvingcos Q = − z/J Q i ≥
0. The scattering mean free path is then just a/ℓ = 4 Q i / √
3. Approximate solutions are given in theAppendix B.A word of caution is here necessary. Inspection of (39)in the absence of disorder (Σ = 0) shows that the allowedenergy range is restricted to − J ≤ E ≤ J . A part ofthe positive energy sector is thus missed and will keepmissed at weak enough disorder. As a consequence us-ing Eq. (39) is only well adapted to the negative energysector. The physical reason is that the honeycomb lat-tice has a two-point Bravais cell and the prescription (38)amounts to consider a symmetric combination of ampli-tudes associated to a -sites and b -sites within a Bravaiscell of the initial slice. Would one had chosen the anti-symmetric combination, then one would have gottenΨ N = J Z π dq π e i Nq z − J (1 + 2 cos q ) , (42)which is well adapted to the positive energy sector. Asa rule of thumb, we thus only use Eqs. (38) and (39) forthe negative energy sector E = −| E | and then resort to(16) whenever necessary.Once the scattering mean free path ℓ ( E ) and the scat-tering mean free time τ ( E ) have been calculated, one cancompute the ratio v ( E ) = ℓ ( E ) /τ ( E ). From (41), we get v ( E ) c = 2 √ Q r sinh Q i Q i , (43)where c is the Dirac fermions speed and where 0 ≤ Q r ≤ π/ Q i ≪
1, and thus v ( E ) /c ≈ (2 sin Q r ) / √
3, where Q r solves ε ( Q r ) = − J (1 + 2 cos Q r ) ≈ E − ReΣ( E ) . (44)In the weak disorder regime, the real part of self-energyis generally small compared to the value of E and can beusually discarded.Returning to fully dimensioned quantities, we thus findthe usual result that v ( E ) = | ∂ k x ε ( k x ) | / ~ is the groupvelocity (here along Ox ) when disorder is sufficientlyweak. Quantitative comparison between the SCBA andRGF results are shown in Figs. 18 and 19. The twomethods show remarkable agreement. We note howeverthat the SCBA underestimates ℓ in the quadratic disper-sion regime but overestimates ℓ in the linear dispersionregime. The deviation of the Born approximation fromthe true (RGF) ℓ might be understood from two perspec-tives: first as a consequence of the smoothing of the DoS ν ( E ) by disorder (see Fig. 20), second as a consequenceof the enhancement of the group velocity (see Fig. 21and 22). Indeed, from the relation ℓ = vτ ( v is the groupvelocity), we know that 1 /ℓ is directly proportional toImΣ( E ), which is in turn directly proportional to theDoS. As disorder is increased, the sharp variations of theDoS at the band edges ( E = ± J ), at the van Hove sin-gularities ( E = ± J ), and at the charge neutrality point( E = 0) are smoothed. Since the DoS area must be con-served (particle number is conserved), this smoothing isaccompanied by a redistribution of states over energiesand by an increase of the DoS at some energies. For ex-ample, since the van Hove singularity peaks in the DoSare decreased, there is a corresponding increase of theDoS in the linear dispersion regime. Similarly, as theband edges are smoothed, the DoS at | E | < J decreasesbut increases for | E | > J . Therefore, the departureof the actual 1 /ℓ from the Born value as disorder is in-creased, and whether the Born value underestimates oroverestimates the actual 1 /ℓ in Fig. 18 and Fig. 19, is asimple consequence of the smoothing effect of the DoS indifferent energy regimes.
3. Comparison with the self-consistent theory of localization
The self-consistent theory of localization (SCTL) [12]provides a self-consistent recipe which extends the weak-disorder diagrammatic approach to the regime of Ander-son localization. The diagrammatic approach describesperturbatively the weak localization corrections to clas-sical transport due to interference along closed loops. Asthis correction diverges for infinite system size, the self-consistent theory of localization aims at computing theminimum size of the system beyond which the correc-tion is strong enough to stop the diffusive transport andidentifies this length scale with the localization length.For isotropic scattering and an isotropic dispersion re-lation, SCTL establishes a simple link between the local-ization length ξ , the scattering mean free path ℓ and the4 FIG. 18. (Color online) Inverse of the scattering mean freepath ℓ (in units of the inverse of lattice constant a ) extractedfrom the recursive Green’s function (RGF) method and fromthe self-consistent Born approximation (SCBA) as a func-tion of the square of the disorder strength W in units ofthe hopping energy J . To compare results obtained for thehoneycomb and square lattices, ( W/J ) is renormalized by afactor proportional to the corresponding densities of states, ν ho = 0 .
14 (honeycomb) and ν sq = 0 .
081 (square), in the ab-sence of disorder as defined by Eq. (47). For both lattices, the(negative) energy E is chosen near the band edge where thedispersion is quadratic. The inset identifies the different lat-tices, configurations, energies and calculation methods. Theoverlap between the results of the armchair (AC) honeycomband square (SQ) lattices show that the average propagationnear the band edge is independent of the lattice type andsolely determined by the quadratic nature of the dispersionrelation. Comparison with the Born approximation is shown.We also plot the scattering mean free path corresponding tothe numerically observed localization length (calculated withthe Recursive Green Function method), assuming that thetwo quantities are connected by eq. (45) derived from theSelf-Consistent Theory of Localization. This shows that thepredictions of the Self-Consistent Theory of Localization arequalitatively correct, especially for weak disorder, but thatlarge quantitative deviations are observed at large disorder. wavevector κ , ξ = 2 ℓ p e πκℓ − ≈ ℓ e π κℓ , (45)where κ = | k | when the energy is chosen near the bandedges and κ = | q | = | k − K | when the energy is chosennear the charge neutrality point. SCLT being valid when κℓ ≫
1, we further get the approximation πκℓ ≈ ln( πκξ/ − ln (cid:0) ln( πκξ/ (cid:1) . (46)We show this SCTL prediction in Fig. 18 and Fig. 19, ξ being evaluated by the RGF method. Much to oursurprise, the agreement between our numerical data andthe SCTL prediction proves very poor. It should be re-minded however that SCLT is supposed to be valid when κℓ ≫ , that is in the weak disorder limit where a/ℓ is FIG. 19. (Color online) Inverse of the scattering mean freepath ℓ (in unit of the inverse of the lattice constant a ) as afunction of the square of the disorder strength W in unitsof the hopping energy J for the AC honeycomb lattice. The(negative) energy E = − . J is chosen in the linear disper-sion regime. Connected filled black circles: RGF results. Bluedashed line: SCBA results. Black solid line: Born approxi-mation (BA). Connected filled red squares: scattering meanfree path estimated from the localization length, assumingthat they are connected by eq. (45) derived from the Self-Consistent Theory of Localization. While the SCBA predic-tion agrees well with the RGF numerical computation of themean free path, the predictions of the Self-Consistent Theoryof Localization are quantitatively off. simply proportional to W . Fig. 18 (and to a lesser extentFig. 19 too) indeed shows that the SCLT estimate tendsto be a linear function of W . In fact all predictions,SCBA, RGF, SCTL and Born approximation agree wellwhen W ≪ J . It would be desirable to have numericalresults at smaller disorder strengths. Unfortunately, thelocalization length is too large to be measurable. In anycase, our numerical results show that higher orders in W are completely different for the true (RGF) a/ℓ and theSCLT prediction. Furthermore, the trend in the lineardispersion regime is completely different. These resultsare yet to be understood. C. Density of states
The disorder-averaged DoS per lattice site is definedby ν ( E/J, W/J ) = 12 N c X µ δ ( E/J − λ µ /J ) , (47)where N c is the total number of Bravais cells and where λ µ is the eigenvalue of the Hamilton operator H inEq. (1). This definition is related to the diagonal ele-ments of the Green’s function by ν ( E/J, W/J ) = − Jπ lim η → + Im h j | G ( E + i η ) | j i , (48)5 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 000.511.52 W = 0.5 JW = 1.0 JW = 1.4 J -3 -2.9 -2.80.250.5 -0.5 -0.4 -0.30.250.5 π ν ( E / J , W / J ) E/J M = 200 FIG. 20. (Color online) Density of states (DoS) ν ( E/J, W/J ),Eq. (48) as a function of (negative) energy E (in units of thehopping energy J ) for several disorder strengths W . Thenumber of sites in the transverse direction is M = 200. TheDoS is increasingly smoothed by disorder as W is increased.Because of particle number conservation, the area under thecurve is also a conserved quantity and states are simply re-distributed over a broader energy range. For example, as thesmoothed van Hove singularity peak at E = − J decreases,the corresponding states are redistributed in the wings as in-dicated by the arrows. This translates into an increased DoSnear the charge neutrality point E = 0. Similarly, states atthe band edge at E = − J are partly redistributed outsidethe energy band of the clean system ( E < − J ), resultingin a decrease of the DoS in the quadratic dispersion regime E & − J . The upper inset shows the variations of the DoS inthe quadratic dispersion regime while the lower inset showsthe variations of the DoS in the linear dispersion regime nearthe charge neutrality point. where j labels an arbitrary lattice site in the infinitely-large lattice. Like what we did to extract ℓ with the RGFmethod, we compute G ( E + i η ) for η/J ∈ [0 . , .
07] andthen perform a quadratic fit in η to extract the limit η → + . For each value of η , the longitudinal length L N and the transverse width L M are both chosen to begreater than 10 ℓ ( η ). The value of h j | G ( E + i η ) | j i is thenestimated by considering any site j at a minimum dis-tance of 5 ℓ ( η ) from the two ends of the tube.The comparison between the DoS calculated with theRGF and SCBA methods is shown in Fig. 23. The twomethods agree remarkably well except near the bandedges ( E = ± J ) and at the charge neutrality point( E = 0). The breakdown of SCBA is probably due tothe fact that qℓ ≪
1. For moderate disorder strength( W = 1 . J ), the RGF results display Lifshitz tails nearthe band edges while the SCBA results show a square-root cut-off.The DoS at E = 0 as a function of W is shown inFig. 24. The large deviation between the RGF and SCBAresults is expected since the latter is strictly not valid inthe range of disorder strengths shown. Note that Shon v / v ( W/J ) E = − . J FIG. 21. Variation of the group velocity v = ℓ/τ with respectto the square of the disorder strength W in units of of thehopping energy J within the SCBA. The (negative) energy E = − . J is chosen near the linear dispersion regime. Thegroup velocity of the clean system (chosen along Ox ) is v =1 . c , where c is the massless Dirac particles velocity. As onecan see, the group velocity near the charge neutrality pointis only slightly enhanced by disorder. The reason why SCBAoverestimates ℓ near the charge neutrality point compared tothe RGF calculation is thus essentially a consequence of thesmoothing of the density of states by disorder, see Fig. 20. v / v ( W/J ) E = − . J FIG. 22. Variation of the group velocity v = ℓ/τ with respectto the square of the disorder strength W in units of of thehopping energy J within the SCBA. The (negative) energy E = − . J is chosen in the quadratic dispersion regime. Thegroup velocity of the clean system (chosen along Ox ) is givenby v = 0 . c , where c is the massless Dirac particles veloc-ity. As one can see, the group velocity near the band edge issignificantly enhanced by disorder. This is the main reasonwhy SCBA underestimates ℓ near the band edge compared tothe RGF calculation. -4 -2 0 200.511.52 W =1.4 JW =0.1 J SCBA, W =1.4 J SCBA, W =0.1 J π ν ( E / J , W / J ) E/J M = 200 FIG. 23. (Color online) Density of states ν ( E/J, W/J ) com-puted by the recursive Green’s function method and bythe self-consistent Born approximation at two weak disorderstrengths W . The number of sites in the transverse directionis M = 200. Since ν ( E/J, W/J ) is an even function of E , wehave plotted the RGF (symbols) and SCBA (continuous line)predictions at W = 1 . J in the negative energy sector. Thecorresponding predictions at W = 0 . J have been plotted inthe positive energy sector. The agreement is excellent. RGFSCBA π ν ( E / J , W / J ) W/J E = 0 M = 400 FIG. 24. (Color online) Density of states ν ( E/J, W/J ) com-puted by the recursive Green’s function method (full blackcircles) and by the self-consistent Born approximation (redsquares) as a function of the disorder strength W (in units ofthe hopping energy J ) at the charge neutrality point E = 0.SCBA underestimates ν ( E/J, W/J ) by more than a factor 2when disorder is large enough. This is the main reason whySCBA overestimates ℓ (see text). and Ando [35] obtained slightly different SCBA resultsfor the DoS at the Dirac points probably because theyused a strictly linear dispersion relation, hence discard-ing trigonal warping [41]. As already witnessed by the re-sults on ℓ , the SCBA underestimates ν ( E ), and hence theself-energy Σ( E ), in the linear dispersion regime. How- ever, although quantitatively incorrect, the SCBA qual-itatively captures an essential property of the system,namely the very fast decrease of the DoS when W goesto zero, essentially like exp( − A/W ) . VI. CONCLUSION
In this paper, we have studied coherent transportin the honeycomb lattice subjected to the effect of aspatially-uncorrelated on-site disorder with symmetricbox-like distribution. We have used the recursive Green’sfunction method to reliably extract the scattering meanfree path ℓ and the density of states. We have comparedthese quantities to the analytic predictions of the self-consistent Born approximation and found good agree-ment at weak disorder. We have also used the recur-sive Green’s function method to extract the localizationlengths for different transverse sizes. We have shownthat all of these finite-size localization lengths can becollapsed onto a single curve. We have checked that thiscurve is universal as it applies equally well to the squareand honeycomb lattices at any energy. In particular, itapplies to the honeycomb lattice at the charge neutral-ity point, at the van Hove singularities, and at the bandedges. These findings validate the one-parameter scalinghypothesis which is thus not restricted to particles witheither quadratic dispersion or linear dispersion relations. ACKNOWLEDGMENTS
LKL acknowledges support from the French Merlion-PhD programme (CNOUS 20074539) and would like tothank M. Schreiber for sharing of data. BG and ChMacknowledge support from the LIA FSQL and from theFrance-Singapore Merlion programme (FermiCold grant2.01.09). The Centre for Quantum Technologies is a Re-search Centre of Excellence funded by the Ministry ofEducation and the National Research Foundation of Sin-gapore. ChM is a Fellow at the Institute of AdvancedStudies (NTU, Singapore).
Appendix A: Diagonal elements of the clean latticeGreen’s function
The diagonal elements of the disorder-free latticeGreen’s function G turn out to be independent of thesite iI ( z ) = h i | G ( z ) | i i = Z B d k Ω zz − J | f ( k ) | , (A1)with Ω = 8 π / (3 √ a ) the area of the Brillouin zone.We use the rescaling √ k x a/ → α and 3 k y a/ → β .Defining Z = z/J and I ( z ) = H ( Z ) /J , little algebra7 . . . . . Re Z . . . . . . . . I m Z Im µ < FIG. 25. The boundary Im µ = 0 (open circles) divides thefirst quadrant of the complex- Z plane in two regions. The leftregion satisfies Re Z < µ <
0. Because of the analyt-ical continuation of the elliptic integral across the boundaryIm µ = 0, different analytic expressions have to be used tocompute the rescaled diagonal element of the Green’s function H ( Z ) in these two regions, see Eq. (A3). gives H ( Z ) = Z Z π dαdβπ ZZ − (1 + 4 cos α + 4 cos α cos β ) , (A2)where Z can be any point in the complex plane exceptthe real interval [ − , Z outside thisinterval and do an analytic continuation in the complexplane. In fact, direct inspection shows that Re H is evenin Re Z and in Im Z , while Im H is odd in Re Z and inIm Z . This means that it is sufficient to compute H ( Z )is the first quadrant (Re Z ≥
0, Im Z ≥
0) of the complexplane and use these parity properties to infer H ( Z ) in theother quadrants. By demanding I ( Z ) to vary smoothlyas Z varies in the complex plane and noticing how µ (see below) crosses the branch cut of the elliptic integral,we get [33, 34] H ( Z ) = iΓ g (Γ) K ( µ ) for Z = iΓ , Γ > Z g ( Z ) (cid:16) K ( µ ) + 2i K (1 − µ ) (cid:17) for Im µ ≤ < Re Z ≤ ,Z g ( Z ) K ( µ ) otherwise , (A3) where g ( Z ) = 2 π ( Z − ( Z + 3) , (A4a) µ = 16 Z ( Z − ( Z + 3) , (A4b) g (Γ) = − π (Γ + 1) (Γ + 9) , (A4c) µ = 16 − (cid:16)p (Γ + 9)(Γ + 1) − (Γ + 1) (cid:17) + 1) (Γ + 9) (A4d)and K ( ρ ) = Z π dθ p − ρ sin θ (A5)is the complete elliptic integral of the first kind [42]. Fig-ure 25 gives the boundary Im µ = 0 in the first quadrantof the complex- Z plane. Appendix B: Self-consistent Born approximation1. Van Hove singularities
With the parametrization Σ = γJ e − i θ , one expects γ ≪ θ ≈ π/ γ ln (cid:18) γ (cid:19) = 64 πJ W (B1)Thus, γ ≈ W πJ Ω (cid:18) πJ W (cid:19) . (B2)in terms of the Lambert W -function Ω( α ) defined for anycomplex number α through the identity α = Ω( α )e Ω( α ) . (B3)At weak disorder, the asymptotic form for | α | ≫ α ) = ln α − ln ln α + O (cid:18) ln ln α ln α (cid:19) (B4)yields Eq. (27).The self-energy in fact obeys the more general relationΣ SCBA = 4 Jρ + i u Ω (cid:0) e − i2 π/ ( ρ + i u ) (cid:1) , (B5)where ρ = (1 − i4 π/
3) and u = 64 πJ /W . Assuming | ρ | ≪ u , we get from the latter equation θ = π/ δθ with δθ ≈ π (cid:0) ln( πJ W ) − ln ln( πJ W ) (cid:1) . (B6)8
2. Band edges
We consider the dimensionless “detuning” δ = 3 − E/J from the band edge at 3 J . We consider Z = δ + Σ /J ≪ Z approximately obeysthe equation Z = β (cid:16) ln(12 /Z ) − i π (cid:17) + δ, (B7)where β = √ W πJ . Using the Lambert W -function, wefind Σ ≈ βJ Ω (cid:0) − δ/β β (cid:1) − δJ. (B8)The large-argument expansion of the Lambert W -function yieldsΣ ≈ βJ (cid:18) ln( 12 β ) − i π − ln (cid:20) ln( 12 β ) − i π + δ/β (cid:21)(cid:19) . (B9)
3. Scattering mean free path
For
E <
0, we solve Eq. (41) to obtain the followingapproximate solutions: aℓ = γ ( E ) √ − Re( Z ) , | E | ≪ J, √ γ ( E ) , E = − J, √ δ γ ( E ) , δ = 3 + E/J ≪ . (B10)We note that the factor q − Re( Z ) is needed whenthe energy is sufficiently far from the charge neutralitypoint (where deviations from the linear dispersion regimecome into play) but not too near the van Hove singu-larity (where a different approximation applies). Indeedthe real part of the self-energy is then no longer small compared to the energy E itself, at least for the disor-der strengths considered. This is the case at E/J = − . Appendix C: Derivation of the generalized version ofthe recursive Green’s function method
Starting with Eq. (30), the Born series for the Green’sfunction G N reads G N = G N − + G N − ( H hop N + H slice N ) G N , (C1a) G N = G N − + G N ( H hop N + H slice N ) G N − , (C1b)where H hop N = H hop N − ,N + H hop N,N − . In Eq. (C1a), theterm G N − H slice N G N vanishes since H slice N couples onlysites within slice N and one immediately gets Eq. (31a)after sandwiching with bras and kets of sites not exceed-ing the ( N − G N − does not couple to sites in slice N , we immedi-ately recover Eq. (31b) by setting n = N in Eq. (31a).Eq. (31c) can be obtained from Eq. (C1b) through a sim-ilar procedure.To obtain Eq. (31d), we need to go back to the defi-nition of the Green’s function ( E − H N ) G N = , whichimplies ( E − H N − − H hop N − H slice N ) G N = . (C2)Sandwiching Eq. (C2) with bras and kets of all siteswithin slice N , we get EG ( N ) N,N − H hop N,N − G ( N ) N − ,N − H slice N G ( N ) N,N = . (C3)As Eq. (31b) yields G ( N ) N − ,N = G ( N − N − ,N − H hop N − ,N G ( N ) N,N , (C4)substitution into Eq. (C3) immediately gives Eq. (31d). [1] P. W. Anderson, Phys. Rev. , 1492 (1958).[2] M. P. V. Albada and A. Lagendijk,Phys. Rev. Lett. , 2692 (1985).[3] Y. Bidel, B. Klappauf, J. C. Bernard, D. Delande,G. Labeyrie, C. Miniatura, D. Wilkowski, and R. Kaiser,Phys. Rev. Lett. , 203902 (2002).[4] H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, andB. A. van Tiggelen, Nat. Phys. , 945 (2008).[5] S. Faez, A. Strybulevych, J. H. Page,A. Lagendijk, and B. A. van Tiggelen,Phys. Rev. Lett. , 155703 (2009).[6] M. Belzons, E. Guazzelli, and O. Parodi, J. Fluid Mech. , 539 (1988).[7] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht,P. Lugan, D. Cl´ement, L. Sanchez-Palencia, P. Bouyer,and A. Aspect, Nature , 891 (2008). [8] J. Chab´e, G. Lemari´e, B. Gr´emaud, D. De-lande, P. Szriftgiser, and J. C. Garreau,Phys. Rev. Lett. , 255702 (2008).[9] Z. Cheng and S.-W. Gu, Phys. Rev. B , 3128 (1990).[10] S. Ilani, J. Martin, E. Tetelbaum, J. Smet, V. Umansky,D. Mahalu, and A. Yacoby, Nature , 328 (2004).[11] E. Abrahams, P. W. Anderson, D. C. Licciardello, andT. V. Ramakrishnan, Phys. Rev. Lett. , 673 (1979).[12] D. Vollhardt and P. W¨olfle,Phys. Rev. B , 4666 (1980); P. W¨olfle and D. Voll-hardt, Int. J. Mod. Phys. B , 1526 (2010).[13] R. C. Kuhn, O. Sigwarth, C. Miniatura, D. Delande, andC. A. M¨uller, New J. Phys. , 161 (2007).[14] P. M. Ostrovsky, I. V. Gornyi, and A. D. Mirlin,Phys. Rev. B , 235443 (2006).[15] A. D. Mirlin, F. Evers, I. V. Gornyi, and P. M. Ostro- vsky, Int. J. Mod. Phys. B , 1577 (2010).[16] 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).[17] Y. Zhang, Y. Tan, H. Stormer, and P. Kim, Nature ,201 (2005).[18] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, , andA. A. Firsov, Nature , 197 (2005).[19] E. R. Mucciolo and C. H. Lewenkopf, J. Phys.: Condens.Matter , 273201 (2010).[20] K. L. Lee, B. Gr´emaud, R. Han, B.-G. Englert, andC. Miniatura, Phys. Rev. A , 043411 (2009).[21] P. Soltan-Panahi, J. Struck, P. Hauke, A. Bick,W. Plenkers, G. Meineke, C. Becker, P. Windpassinger,M. Lewenstein, and K. Sengstock, Nat. Phys. , 434440(2011).[22] L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu, andT. Esslinger, Nature , 302305 (2012).[23] G. Labeyrie, D. Delande, C. A. M¨uller, C. Miniatura,and R. Kaiser, Europhys. Lett. , 327 (2003);Phys. Rev. A , 033814 (2003).[24] D. Cl´ement, A. F. Var´on, M. Hugbart, J. A.Retter, P. Bouyer, L. Sanchez-Palencia, D. M.Gangardt, G. V. Shlyapnikov, and A. Aspect,Phys. Rev. Lett. , 170409 (2005).[25] J. E. Lye, L. Fallani, M. Modugno, D. S.Wiersma, C. Fort, and M. Inguscio,Phys. Rev. Lett. , 070401 (2005).[26] F. Jendrzejewski, A. Bernard, K. M¨uller, P. Cheinet,V. Josse, M. Piraud, L. Pezz´e, L. Sanchez-Palencia,A. Aspect, and P. Bouyer, Nat. Phys. , 398403 (2012).[27] M. Greiner, O. Mandel, T. Esslinger, T. W. H¨ansch, and I. Bloch, Nature , 39 (2002).[28] M. Schreiber and M. Ottomeier, J. Phys.: Condens. Mat-ter , 1959 (1992).[29] S.-J. Xiong and Y. Xiong,Phys. Rev. B , 214204 (2007).[30] A. Lherbier, B. Biel, Y.-M. Niquet, and S. Roche,Phys. Rev. Lett. , 036803 (2008).[31] R. Kubo, J. Phys. Soc. Japan , 570 (1957).[32] E. Akkermans and G. Montambaux, Mesoscopic Physicsof Electrons and Photons , 1st ed. (Cambridge UniversityPress, 2007).[33] T. Morita and T. Horiguchi, J. Math. Phys. , 981(1971).[34] T. Horiguchi, J. Math. Phys. , 1411 (1972).[35] N. H. Shon and T. Ando, J. Phys. Soc. Japan , 2421(1998).[36] A. MacKinnon and B. Kramer, Z. Phys. B , 1 (1983).[37] G. D. Mahan, Many Particle Physics , 3rd ed. (Springer,2000).[38] A. MacKinnon, Z. Phys. B , 385 (1985).[39] D. Wurtz and B. M. Pohlmann, J. Phys. C , 5631(1988).[40] C. A. M¨uller and D. Delande, in Les Houches 2009: Ul-tracold Gases and Quantum Information (Oxford Uni-versity Press, 2011) pp. 442–533.[41] A. H. Castro Neto, F. Guinea, N. M. R.Peres, K. S. Novoselov, and A. K. Geim,Rev. Mod. Phys. , 109 (2009).[42] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals,Series, and Products , 7th ed., edited by A. Jeffrey andD. Zwillinger (Academic Press, 2007).[43] R. Corless, G. Gonnet, D. Hare, D. Jeffrey, andD. Knuth, Adv. Comput. Math.5