Dislocation and node states in bilayer graphene systems

TL;DR

This paper studies the electronic properties of partial dislocation networks in bilayer graphene, revealing charge accumulation, current flow patterns, and localization behaviors influenced by dislocation configurations.

Contribution

It provides new insights into how dislocation networks affect electronic states and current flow in bilayer graphene, highlighting the role of partials and intersections.

Findings

Charge accumulation at partials depends on Fermi energy and Burgers vector.

Current flows spiral along dislocation lines with interlayer components.

Localization shifts from partials to intersections near the Dirac point.

Abstract

We investigate the electronic structure of realistic partial dislocation networks in bilayer graphene that feature annihilating, wandering, and intersecting partial lines. We find charge accumulation states at partials that are sensitive to Fermi energy and partial Burgers vector but not to the screw versus edge character of the partial. These states are shown to be current carrying, with the current density executing a spiral motion along the dislocation line with a strong interlayer component to the current. Close to the Dirac point localization on partials switches to localization on intersections of partials, with a corresponding complex current flow around the nodes.

Figures

Tables

Table 1. Table 1: M 𝑀 M matrices for the C 3 subscript 𝐶 3 C_{3} star of Bernal stacked AB and AC stacked bilayer graphene.

Stacking	$M_1$	$M_2$	$M_3$
AB
AC

Equations

$$H(\mathbf{r},\mathbf{p})=\left(\begin{array}{cc}{v}_F\mathbf{\sigma }.\mathbf{p}& S(\mathbf{r})\\ S(\mathbf{r}{)}^{†}& v_F{\mathbf{\sigma }}^{\ast }.\mathbf{p}\end{array}\right)$$

$$S(\mathbf{r})=\frac{1}{V_{UC}}\sum _i{M}_i{t}_{\perp }(K_i^2)e^{i{\mathbf{G}}_i.\Delta \mathbf{u}(\mathbf{r})},$$

$$H(\mathbf{p},{\mathbf{u}}_1)=\left(\begin{array}{cc}{v}_F\mathbf{\sigma }.\mathbf{p}& S({\mathbf{u}}_1)\\ S({\mathbf{u}}_1{)}^{†}& v_F{\mathbf{\sigma }}^{\ast }.\mathbf{p}\end{array}\right),$$

$${\mathbf{u}}_1=\frac{{\mathbf{B}}_i}{2}\left(1+\frac{\tanh \gamma x}{\tanh \gamma L_0/2}\right)$$

$$S_{AB}$$

$$S_{AC}$$

$$S_{AA}$$

$$\mathbf{j}(\mathbf{r})=\Im \Psi (\mathbf{r}{)}^{†}\left[{\nabla }_{\mathbf{p}}H(\mathbf{r},\mathbf{p})\Psi (\mathbf{r})\right]$$

$${\mathsf{j}}_{\mathsf{l}}^{(\mathsf{z})}({\mathbf{R}}_{\mathsf{i}})=\frac{1}{{\mathsf{V}}_{\mathsf{UC}}}\sum _{{\alpha }^{\prime }\in \mathsf{l}}\left[\mathsf{n}({\mathbf{R}}_{\mathsf{i}}+{\mathbf{\nu }}_{{\alpha }^{\prime }}){\mathsf{v}}_{\mathsf{z}}+{\mathsf{v}}_{\mathsf{z}}\mathsf{n}({\mathbf{R}}_{\mathsf{i}}+{\mathbf{\nu }}_{{\alpha }^{\prime }})\right]$$

$$t_{\alpha \beta }^{HS}({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha },{\mathbf{R}}_j+{\mathbf{\nu }}_{\beta })+\delta t_{\alpha \beta }^{LS}({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha },{\mathbf{R}}_j+{\mathbf{\nu }}_{\beta }),$$

$$\begin{array}{rl}{H}_{TB}& =\sum _{\alpha {\mathbf{R}}_i\beta {\mathbf{R}}_j}{t}_{\alpha \beta }({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha },{\mathbf{R}}_j+{\mathbf{\nu }}_{\beta })\\ & \times c_{{\mathbf{R}}_j+{\mathbf{\nu }}_{\beta }\beta }^{†}{c}_{{\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha }\alpha }.\end{array}$$

$$|{\Phi }_{{\mathbf{k}}_1\alpha }⟩=\frac{1}{\sqrt{N}}\sum _{{\mathbf{R}}_i}{e}^{i{\mathbf{k}}_1.({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha })}|{\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha }⟩$$

$$|{\Psi }_l⟩=\sum _{{\mathbf{k}}_1\alpha \in l}{c}_{{\mathbf{k}}_1\alpha }|{\Phi }_{{\mathbf{k}}_1\alpha }⟩.$$

$$⟨{\Psi }_l\left|j_l^{(z)}({\mathbf{R}}_i)\right|{\Psi }_{l^{\prime }}⟩=\frac{1}{V_{UC}}$$

$$\times \sum _{{\mathbf{k}}_1,{\mathbf{k}}_2,\alpha \in l,\beta \in l^{\prime }}\left[⟨{\Phi }_{{\mathbf{k}}_1\alpha }\left|n({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha })v_z\right|{\Phi }_{{\mathbf{k}}_2\beta }⟩+\text{H.c.}\right].$$

$$t_{\alpha \beta }(\mathbf{r},\mathbf{\delta })=\int d{\mathbf{q}}^{\prime }d\mathbf{q}{e}^{-i{\mathbf{q}}^{\prime }.\mathbf{r}}{e}^{-i\mathbf{q}.\mathbf{\delta }}{t}_{\alpha \beta }({\mathbf{q}}^{\prime },\mathbf{q}),$$

$$\frac{(z_l-z_{l^{\prime }})}{V_{UC}Ni\hslash (2\pi {)}^{2d}}\int d{\mathbf{q}}^{\prime }d\mathbf{q}{t}_{\alpha \beta }({\mathbf{q}}^{\prime },\mathbf{q})e^{-i({\mathbf{k}}_1+{\mathbf{q}}^{\prime }-\mathbf{q}).({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha })}$$

$$\times e^{i({\mathbf{k}}_2-\mathbf{q}).{\mathbf{\nu }}_{\beta }}\sum _{{\mathbf{R}}_j}{e}^{i({\mathbf{k}}_2-\mathbf{q}).{\mathbf{R}}_j}$$

$$\frac{1}{\sqrt{V}}{e}^{-{\mathbf{k}}_1.({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha })}$$

$$\times \left[\frac{\text{sign}(z_l-z_l^{\prime })}{i\hslash V_{UC}^{\perp }}\sum _{{\mathbf{G}}_j}{t}_{\alpha \beta }({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha },{\mathbf{k}}_2+{\mathbf{G}}_j){\left[M\right]}_{\alpha \beta }\right]$$

$$\times \frac{1}{\sqrt{V}}{e}^{{\mathbf{k}}_2.({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha })}$$

$${\varphi }_{\mathbf{k}\alpha }(\mathbf{r})=\frac{1}{\sqrt{V}}{e}^{i\mathbf{k}.\mathbf{r}}|1_{n\alpha }⟩,$$

$$\frac{\text{sign}(z_l-z_l^{\prime })}{i\hslash }{\varphi }_{{\mathbf{k}}_1\alpha }^{†}(\mathbf{r})H({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha },\mathbf{p}){\varphi }_{{\mathbf{k}}_2\beta }(\mathbf{r})$$

$$H(\mathbf{r},\mathbf{p})=\frac{1}{V_{UC}^{\perp }}\sum _{{\mathbf{G}}_j}{\left[M\right]}_{\alpha \beta }{t}_{\alpha \beta }(\mathbf{r},{\mathbf{K}}_j+\mathbf{p})$$

$$j_1^{(z)}(\mathbf{r})=\frac{1}{i\hslash }\left[-{\Psi }_1^{†}(H_{12}{\Psi }_2)+(H_{12}{\Psi }_2{)}^{†}{\Psi }_1\right]$$

$$j_1^{(z)}(\mathbf{r})=-j_2^{(z)}(\mathbf{r}).$$

$${\partial }_t{n}_1(\mathbf{r})$$

$$\left(\begin{array}{cc}{H}_{11}& H_{12}\\ H_{21}& H_{22}\end{array}\right)\left(\begin{array}{c}{\Psi }_1\\ {\Psi }_2\end{array}\right)=i\hslash {\partial }_t\left(\begin{array}{c}{\Psi }_1\\ {\Psi }_2\end{array}\right)$$

$${\partial }_t{n}_1(\mathbf{r})+\nabla .\mathbf{j}(\mathbf{r})+j_1^{(z)}(\mathbf{r})=0$$

$$\int d\mathbf{r}{j}_l^{(z)}(\mathbf{r})=0$$

Full text

I Introduction

Bernal stacked bilayer graphene is an imperfect material. Transmission electron microscopy (TEM) experiments have revealed the existence of a network of partial dislocations extending throughout the materialald13 ; butz14 , which has been shown to impact electronic properties dramatically in further investigationskiss15 ; 2 ; 3 ; 8 ; 12 ; 4 ; 6 ; 9 ; shall17 . In fact, even a single partial dislocation can qualitatively change transport at the Dirac point, driving a nominally minimally metallic sample to an insulating stateshall17 . Alternatively, in a semi-conducting bilayer, partial dislocations provide a network of conducting channels through topologically protected bound states4 ; 6 , creating a metal out of a semi-conductor. Such extended defects arising from the stacking degree of freedom are likely to be found in all van der Walls bound multilayers – a broad class of materials. Bilayer graphene, probably the simplest such material, thus offers an ideal test ground for elucidating the impact of this important defect on low-dimensional systems.

However, even for bilayer graphene the electronic structure of partial dislocations is not fully understood, a situation attributable to the formidable computation challenge posed by simulating a realistic dislocation network. For example, in graphene grown epitaxially on SiC the domain size is of the order of $50\,$nm and the partial width of the order of $5\,$nm. A network of $\sim 20$ domains then requires of the order of a square micrometer, equivalent to $\sim 10^{8}$ carbon atoms. This lies well beyond the scope of any atomistic approach and taxes even continuum methodologies. To date, therefore, most theoretical studies have considered idealized partial dislocations that are (i) straight and (ii) do not intersect with other partial dislocations. The purpose of the present paper is to address the electronic structure of realistic partial dislocation networks in which partials both wander and intersect.

We consider two networks, one derived directly from experimental TEM images – of the order of 1$\mu m^{2}$ in area – and a model network of trigonally intersecting partials 100 times smaller in area. Our principal focus will be the question of localized states on partial dislocations and their intersections. While in a finite layer-perpendicular electric field bound states on partials are topologically protected, there is evidence that they exist even in the absence of such a field8 ; 12 . We find such localized states, and show that their existence is sensitive to the partial Burgers vector and energy, although not to the screw versus edge character of a wandering partial. At energies below the Dirac point so-called type 3 partials feature localized states, switching to type 2 partials above the Dirac point, with close to the Dirac point localization instead occurring on intersections of the network. This coupling of energy and Burgers vector is a robust feature found in both the experimental as well as model network.

We further analyze the current density of partial dislocation networks, finding intense currents both along partials and at nodes, with a strong interlayer component. For example, the localized states on partials of type 2 and 3 are shown to possess a current density spiraling around the dislocation line, with the spiral chirality opposite in opposite K valleys. In general, we find that a current density impinging upon a partial is driven from one layer onto the other and then back again; a truncated version of the “snake states” found in twist bilayer graphenemel10 ; shall10 ; san12 ; shall13 ; shall16 . On the other hand, for the node states found at energies close to the Dirac point we find strong current densities around the vertices of intersecting partials, with a complex pattern of in-plane and interlayer currents.

II Theory

II.1 Continuum description of a

partial network

A partial dislocation in bilayer graphene represents a continuous change in stacking order from AB to AC stacking. Their electronic description therefore requires a stacking potential dependent on the local displacement of the layers at each point of the lattice. Taking AB stacking as the reference state this potential will, in the low-energy Dirac-Weyl sector, be represented by a $2\times 2$ matrix $S({\bf r})$ in a Hamiltonian of the form:

[TABLE]

such that $S({\bf r})$ couples the two Dirac-Weyl operators of each layer. Such an interlayer stacking potential has been derived in Ref. M, and previously employed to treat dislocations in Refs. kiss15, and shall17, . It takes the form

[TABLE]

where $\Delta{\bf u}({\bf r})={\bf u}_{2}({\bf r})-{\bf u}_{1}({\bf r})$ is the local displacement of the two layers. In this expression the $M_{i}$ matrices are given in Table 1, the sum is over the translation group of the high-symmetry K point ${\bf K}_{1}=\frac{4\pi}{3a}(1,0)$ with $a$ the lattice constant of graphene (i.e. ${\bf K}_{i}={\bf K}_{1}+{\bf G}_{i}$ with ${\bf G}_{i}$ a reciprocal vector), and $t_{\perp}(q^{2})$ is the Fourier transform of the interlayer hopping $t_{\perp}(\delta^{2})$. From Eq. (1) we see that the potential $S({\bf r})$ is a matrix whose elements describe the scattering amplitude from sublattice $\alpha$ of layer 1 to sublattice $\beta$ of layer 2; in the Dirac-Weyl language $S({\bf r})$ couples the pseudospin channels of each layer.

Note that the in-plane deformation accompanying a local shift of the two layers is not included in Eq. (1). The impact of such in-plane deformation upon the electronic structure of partial dislocations has been explored in Ref. shall17, and found to be negligible, and so we do not consider the resulting gauge, scalar, and Fermi velocity correction terms in the layer-diagonal blocks of our effective Hamiltoniangupta18 .

In order to establish that this approach can adequately describe the electronic properties of arbitrarily stacked bilayer graphene, we first analyze a very simple test case: a constant global shift of one layer. For a constant shift $\mathbf{u}_{1}$ applied to the first layer the stacking potential in Eq. (2) becomes a constant matrix, with the effective Hamiltonian now being

[TABLE]

which may be directly diagonalized. In Fig. 1 we show the low-energy spectrum of the bilayer for four shift vectors on a path that translates the bilayer from AB to AA stacking: ${\bf u}_{1}={\bf 0}$ (the AB stacked bilayer), ${\bf u}_{1}=1/10({\bf a}_{1}+{\bf a}_{2})$, ${\bf u}_{1}=1/5({\bf a}_{1}+{\bf a}_{2})$, and ${\bf u}_{1}=1/3({\bf a}_{1}+{\bf a}_{2})$ (AA stacked). Tight-binding results are shown by the full (black) lines with the results of the effective Hamiltonian presented with broken (red) lines. As may be seen, an excellent agreement exists between the two methods for (the relevant) energies of $\sim 1\,$eV around the Dirac point, and so Eq. (3) captures the low-energy electronic structure for all mutual translations of the bilayer.

The requirement of a continuous graphene membrane in each layer leads to the condition that separate AB and AC domains be connected by one of three possible partial Burgers vectors: ${\bf B}_{1}=a(1/2,1/(2\sqrt{3}))$, ${\bf B}_{2}=a(-1/2,1/(2\sqrt{3}))$, and ${\bf B}_{3}=a(0,-1/\sqrt{3})$. Traversing from one domain to the other then involves a local shift of one layer by one of these partial Burgers vectors. This transition occurs, according to experiment, over a width of $\sim 5\,$nm and the detailed atomic structure of this transition region has been carefully investigated via semi-empirical bond potential calculations in Ref. butz14, . From this data we are able to extract a model form of ${\bf u}_{1}({\bf r})$,

[TABLE]

with ${\bf B}_{i}$ the partial Burgers vector of the partial dislocation, $L_{0}$ the partial width, and $\gamma$ the partial stiffness, and the values of $L_{0}$ and $\gamma$ chosen to reproduce closely the data of Ref. butz14, . This is shown in Fig. 2 for the $AB\to AC$ transition mediated by ${\bf B}_{3}$ partial Burgers vector. Also shown is the interlayer field $S({\bf r})$, projected onto the 3 distinct high-symmetry stacking types that exist through this transition, $S_{AB}$, $S_{AC}$, and $S_{AA}$:

[TABLE]

(Deployed as a constant interlayer block these matrices will generate the standard AB/AC and AA stacked band structures.) The matrix function $S({\bf r})$ can be seen to transition between AC and AB stacking with the maximum of the AA component in the middle of the partial.

II.2 Current operator in multi-layer systems

In a recent article Gupta et al. gupta18a have presented a derivation of the current operator in the continuum theory we employ here, finding the intuitive form

[TABLE]

with $H({\bf r},{\bf p})$ the effective Hamiltonian. This form, however, only describes the diagonal blocks of Eq. (1), where it reproduces the well known $v_{F}\bm{\sigma}$ current operator of the Dirac-Weyl equation. We therefore require an extension of the concept of the current density to multilayer systems within the continuum approximation. This, however, requires some care. As each layer is represented by a strictly two-dimensional (in this case Dirac-Weyl) operator there is no natural definition of a 3-vector current density. The description of current flow in a multilayer therefore breaks down into an in-plane current density (a 2-vector) and an interlayer current density (a scalar), with the latter representing the net charge flow between the layers at any point ${\bf r}$. In what follows we will derive a general form of this scalar current density for a bilayer system, and show that it satisfies an appropriately generalized continuity equation.

We begin with the tight-binding form of the interlayer current operator

[TABLE]

where the density operator is defined as $\hat{n}({\bf R}_{i}+\bm{\nu}_{\alpha^{\prime}})=\left|{\bf R}_{i}+\bm{\nu}_{\alpha^{\prime}}\right\rangle\left\langle{\bf R}_{i}+\bm{\nu}_{\alpha^{\prime}}\right|$, the velocity operator is $v_{z}=1/(i\hbar)[z,H_{TB}]$, and $V_{UC}$ is the unit-cell volume of the multilayer system tod02 ; boy10 . The Hamiltonian $H_{TB}$ we take to be that of a standard two-center tight-binding approximation in a basis of Wannier orbitals, with an important caveat. We consider that the system of interest be obtained from a deformation applied to some high-symmetry reference state, for instance in this case the high-symmetry AB stacked bilayer. Note that we do not assume that this deformation is small, and thus this approach is not perturbative in nature. We implement this deformation simply by changing the values of the tight-binding integrals while keeping the labels of the local orbitals fixed, which, evidently, represents the choice of a local co-moving coordinate system with the deformation. The tight-binding overlap integrals $t_{\alpha\beta}({\bf R}_{i}+\bm{\nu}_{\alpha},{\bf R}_{j}+\bm{\nu}_{\beta})$ are then just

[TABLE]

where $\alpha$ and $\beta$ are sub-lattice labels, and ${\bf R}_{i}$ and ${\bf R}_{j}$ lattice vectors, of the high-symmetry system. Any spin or angular momentum indices are subsumed into the $\alpha$ and $\beta$. The tight-binding Hamiltonian is then given by

[TABLE]

The advantage of this seemingly baroque construction is that we can then naturally introduce the Bloch functions of the high-symmetry system as a complete basis for the full Hamiltonian. These are given by

[TABLE]

with ${\bf k}_{1}$ the crystal momentum and $\bm{\nu}_{\alpha}$ a basis vector.

We now consider a general state of the system on layer $l$ expressed in the Bloch-function basis

[TABLE]

Using this basis expansion the interlayer current density is given by

[TABLE]

We now manipulate the form of the matrix element $\left\langle\Phi_{{\bf k}_{1}\alpha}\left|n({\bf R}_{i}+\bm{\nu}_{\alpha})v_{z}\right|\Phi_{{\bf k}_{2}\beta}\right\rangle$ toward a continuum representation by introducing a hopping envelope function $t_{\alpha\beta}({\bf r},\bm{\delta})$ that takes values $t_{\alpha\beta}({\bf R}_{i}+\bm{\nu}_{\alpha},{\bf R}_{j}+\bm{\nu}_{\beta})$ for ${\bf r}={\bf R}_{i}+\bm{\nu}_{\alpha}$ and $\bm{\delta}={\bf R}_{j}+\bm{\nu}_{\beta}-{\bf R}_{i}-\bm{\nu}_{\alpha}$. Inserting this function as its Fourier transform

[TABLE]

the matrix element can be written as

[TABLE]

with the factor $(z_{l}-z_{l^{\prime}})$ taken outside the phase sum as all vectors in this expression are 2-vectors with no $z$-component. This may then be evaluated via the Poisson sum relation yielding

[TABLE]

where we have used $(z_{l}-z_{l^{\prime}})/V_{UC}=\text{sign}(z_{l}-z_{l}^{\prime})/V_{UC}^{\perp}$ to relate the volume of the three-dimensional unit cell $V_{UC}$ to the in-plane unit-cell area $V_{UC}^{\perp}$. The plane-wave normalization factor $1/V$ arises from the $(2\pi)^{2}/V_{UC}^{\perp}$ prefactor of the Poisson sum relation along with the $1/N$ factor due to Bloch-function normalization. Introducing spinor plane waves

[TABLE]

where $\left|1_{\alpha}\right\rangle$ is a unit ket in a space with dimensional equal the number of atomic degrees of freedom, i.e. $\left|1_{\alpha}\right\rangle=(0_{1},\ldots,1_{\alpha},\ldots)^{T}$, we then find for the tight-binding matrix element the final result

[TABLE]

where, by writing the Bloch-state momentum ${\bf k}_{2}$ relative to an expansion point ${\bf K}_{0}$ i.e. substituting ${\bf G}_{j}+{\bf k}_{2}={\bf G}_{j}+{\bf K}_{0}+{\bf p}_{2}$ in Eq. (17), and promoting ${\bf p}_{2}$ to an operator, we have identified

[TABLE]

as the effective continuum Hamiltonian described in Ref. M, from which Eq. (2) may be obtained.

Using Eq. (14) immediately yields for the $z$-component of current in layer 1

[TABLE]

where we take the origin of the $z$-coordinate to be within layer 1 pointing toward layer 2. This result is unsurprising, being identical in form to the standard “bond-current” expressiontod02 . However, rather than a current along a bond it represents a current between layers whose coupling is through the continuum field $H_{12}({\bf r})$. The use of such a local coordinate system ensures that, due to the factor $\text{sign}(z_{l}-z_{l^{\prime}})$, we have

[TABLE]

We now demonstrate that this form of the interlayer current satisfies the appropriate generalization of the continuity equation. To that end we consider the time dependence of the density within layer 1

[TABLE]

where we have used

[TABLE]

in order to replace time derivatives. The first term on the right hand side of Eq. (22) is just the negative of the divergence of intralayer current, with the second term just the negative of $j^{(z)}_{1}$, as may be seen from Eq. (20). We then immediately find a generalized continuity equation given by

[TABLE]

with a corresponding relation for layer 2. If $\Psi({\bf r})$ is a stationary state of the system then

[TABLE]

which follows from the Hermiticity of the layer diagonal blocks $H_{ll}$.

III Partial dislocation networks

III.1 Numerical details

We employ as a basis the Dirac-Weyl states of each layer, which has the advantage of compactness as the energy range required for convergence is $\sim 1.5$ the energy window of interest. The interlayer stacking potential $S({\bf r})$ is defined on a fine grid (we use $3000\times 3000$ for the experimental network and $800\times 800$ for the model network) and the matrix elements of $S({\bf r})$ in the Dirac-Weyl basis obtained via FFT. For the large networks that we consider, a basis of up to 25,000 Dirac-Weyl states from each layer is required, necessitating a distributed-memory parallelization approach to diagonalization. For the interlayer tight-binding hopping function we use a Gaussian form $t_{\perp}(\bm{\delta}^{2})=A\exp(-B\bm{\delta}^{2})$, with $A=0.4\,$eV and B=1\,$$\mathrm{\SIUnitSymbolAngstrom}-2 chosen to reproduce the $0.8\,$eV splitting of bonding and anti-bonding bands in the Bernal bilayer.

III.2 Charge pooling on the mosaic network and charge accumulation at partials

We will first consider the electronic structure of a partial dislocation network taken directly from transmission electron microscopy (TEM) imageskiss15 of bilayer graphene grown epitaxially on SiC. This network is illustrated in Fig. 3 with the various partial dislocations colored according to the their Burgers vector. Exploratory calculations of this system have already been presentedkiss15 , here we will investigate in greater detail the electronic structure. The area of the TEM image in experiment was 1$\mu m^{2}$ equivalent to $\approx\!10^{8}$ carbon atoms which, within the continuum approach described in the previous section, is numerically feasible. The interlayer stacking potential $S({\bf r})$, Eq. (2), that depends directly on the local translation $\Delta{\bf u}({\bf r})$, must now encode the mutual translation of the layers that occurs upon crossing any partial line in the network (within the domains of the mosaic structure the function $\Delta{\bf u}({\bf r})$ will obviously be constant). The partial dislocations in Fig. 3 evidently have a mixture of edge and screw dislocation character, corresponding to the cases in which the Burgers vector is, respectively, parallel or perpendicular to the tangent of the partial line. To model this we take the form of the interlayer displacement Eq. (4) shown in Fig. 2 and use this to describe the displacement along the line perpendicular to the partial tangent, forming a superposition of the resulting displacements at each point in space as this procedure is traced out along a partial line. In this way we find a smooth interlayer displacement field describing a network of wandering partials. Note that a local shift by a multiple of the lattice vectors ${\bf a}_{1,2}$ does not alter the interlayer potential, due to the presence of the reciprocal lattice vector in the exponent of Eq. (2).

Of principle interest is the form of the wave functions of the partial network, which are expected to be very different from the wave functions of the structurally perfect bilayer. In Fig. 4(a-d) we show the electron density integrated over a 13 meV window (of the order of the Fermi smearing at 150 K) with this window placed at four different energies: -20 meV, -5 meV, +5 meV, +90 meV. Even within this small energy window, of the order of $10^{3}$ individual eigenstates contribute to the probability density. Quite clearly, the mosaic structure of the bilayer has a dramatic impact on the wave functions. For the density integrated in the energy window situated at -20meV, see panel (a) of Fig. 4, charge accumulation is found on the type 3 partials (compare with Fig. 3) but not on the partials of type 1 and 2. Above the Dirac point this switches to a charge accumulation on type 2 partials, see panels (c) and (d). Close to the Dirac point substantial charge pooling is observed as some segments of the mosaic network have significantly higher density than others (see panels (b) and (c)), a point first noticed in Ref. kiss15, . Intense localization on intersections of the partial network can also be noted. For the energy window of +90meV, shown in panel (d), the pronounced charge pooling seen near the Dirac point is absent, although one still notes substantial density inhomogeneity created by the partial network.

Now turning to the features of the density of states (DOS), we note qualitative differences between pristine AB-stacked graphene bilayer and that of the bilayer with the experimentally observed partial dislocation network. At the Dirac point, the experimental system features a much more pronounced minimum of the DOS, with two side peaks found at approximately $\pm 0.14$ eV.

How much of this structure, and in particular the existence of “hot” partials on which charge density is accumulated, is related to the specific partial network shown in Fig. 4? In this respect it should be noted that charge accumulation on partials occurs on differently oriented segments of a wandering partial, and so the local orientation of the partial with respect to the underlying lattice is evidently only of secondary importance. This is particularly the case for the type 2 partial, as may be seen in panel (c) of Fig. 4.

To investigate further these features we consider a designed hexagonal network of partial dislocations shown in Fig. 3. In this network all partials are globally aligned along the zigzag direction of the underlying lattice, however we introduce some random disorder into the partials such that they are not perfectly straight, and so locally are not necessarily aligned with the zigzag direction.

As may be observed in Fig. 5, qualitatively similar features are seen to those noted in the experimental partial network - in particular the energy order in which the type 3 and 2 partials become ”hot” and accumulate charge is the same; the type 1 partial, as before, does not accumulate charge. The width of the energy window is again $13\,$meV and the position of the energy windows similar to those used in the analysis of the experimental network. This charge accumulation is therefore largely independent of the global details of partial network but sensitive to energy and Burgers vector type. A number of features are much more pronounced in this more ordered partial network, in particular the localization of charge on the nodes, see panel (d) of Fig. 5. A complete sweep of the electron density as a function of energy can be obtained from the electronic supplementary information.

III.3 Current density in the partial network

The strong imprint of the partial network upon electron density raises the natural question as to the corresponding current densities. Using the formalism outlined in Section II.2 we now examine this question. In Fig. 6 we present the current density for those states whose electron density is exhibited in Fig. 5. In each panel is shown the in-plane current for each layer and the interlayer current density. For clarity we also show the sign of the interlayer current. All current densities shown here are, of course, exactly compensated by a corresponding state in the K∗ valley with opposite current.

We first investigate the current density of the charge accumulation states shown in panels (a) and (c) of Fig. 5. As can be seen in the corresponding panels of Fig. 6 both these states are current carrying. This is particularly evident in panel (c) for the type 2 partial dislocation: the interlayer current exhibits an alternating sign along the partial dislocation line that, in conjunction with an in-plane current density alternating between the layers, signifies that the current spirals around the partial, from one layer to the other and back again. The type 2 charge accumulation state thus appears similar to the topological current-carrying state found in the gap generated by application of finite interlayer bias. The results here raise the intriguing possibility of current-carrying partial dislocation states in the absence of such an external field. This can also be seen in panel (a) for the type 3 partial, although the structure is somewhat obscured.

In general we find that the partial dislocations are associated with strong interlayer currents, as may be seen from panels (a)-(d) of Fig. 6. It is evident by comparison of the in-plane currents that the interlayer current is intense near partials. This can be seen for example in panel (a) upon examination of the charge repelling type 1 partial. One notes that flow of charge in layer 1 appears to be interrupted by the partial. However, comparison with the current density in layer 2 and the interlayer current shows that the charge simply tunnels underneath the dislocation. This finding of strong interlayer currents associated with partial dislocations is consonant with the situation found in the other common stacking fault of bilayer graphene, a mutual rotation of the layers. In that case snake states are found that oscillate between the two layers, and the spiral state found on the type 2 partial can be seen as the analogy of this. Evidently, in bilayer graphene stacking modulation generate strong interlayer current densities.

Finally, we consider the situation of the nodal states and the high-energy states that show no pronounced pattern of charge accumulation. In panels (b) and (d) of Fig. 6 we show the current density corresponding to the electron densities shown in panels (b) and (d) of Fig. 5. The nodal states are seen to exhibit intense current patterns in both the in-plane and interlayer currents. Evidently, nodal states are associated with interlayer current loops, quite tightly bound to the nodes. The current density does not, however, exhibit a six-fold pattern corresponding to the basic geometry of the intersection point – the low-symmetry environment of the node due to partial wandering precludes this. For high-energy states, that are neither nodal nor partial dislocation states, the current density is correspondingly featureless. One does however note that the strong interlayer current persists and again generally changes sign across a partial dislocation line.

IV Discussion

We have explored the electronic structure of massive realistic partial dislocation networks featuring wandering partials that both intersect with each other and annihilate in pairs at local defects. Strong charge accumulation at partials is found that is sensitive to energy and Burgers vector but less sensitive to the orientation of the partial line with respect to the underlying lattice or to the screw versus edge character of the partial. These charge accumulation states are all found within $30\,$meV of the Dirac point and are current-carrying with the current density exhibiting a spiral motion around the partial line. We find that partial dislocations are generally associated with strong interlayer current flow, a feature found not only for charge accumulation states but also in states for which intralayer current density is blocked by a partial but “short circuits” by flowing into the other layer. Very close to the Dirac point, charge strongly localizes on the nodes of the partial dislocation network, in particular for more ordered networks but less so for disordered partial networks typically found in experiment. Correspondingly, the current density exhibits a complex current flow around the nodes.

While current-carrying topological states are known to be found in the band gap of a bilayer subjected to interlayer bias, these states require no external bias for their existence. Partial dislocations in realistic networks are thus natural lines for the flow of current density, with potential implications for the quantum Hall effect in mosaic bilayer graphene.