# Impurity effects on Dirac modes in graphene armchair nanoribbons

Yuriy G. Pogorelov<EMAIL_ADDRESS>IFIMUP-IN, Departamento de Física e
Astronomía, Universidade do Porto, Porto, Portugal,  Vadim M. Loktev
<EMAIL_ADDRESS>N. N. Bogolyubov Institute of Theoretical Physics, NAS of
Ukraine, Kyiv, Ukraine,
&
Igor Sikorsky Kyiv Polytechnic Institute, Kyiv, Ukraine,

###### Abstract

We consider finite ribbons of graphene with armchair orientation of their
edges to study in detail impurity effects on specific Dirac-like modes. In the
framework of Anderson hybrid model of impurity perturbation, a possibility for
Mott localization and for opening of a mobility gap under local impurity
perturbations is found and analyzed in function of this model parameters: the
impurity energy level, its hybridization with the host Dirac modes, and the
impurity concentration. Possible electronic phase states in such disordered
system and subsequent phase transitions between them are discussed.

graphene nanoribbons, armchair edge states, impurity disorder, Anderson model

##  I Introduction

Between electronic properties of two-dimensional (2D) graphene layer, the
presence of linear gapless quasiparticle modes, or 2D Dirac modes, is
especially notable by defining unusual physical effects in graphene Geim2004 ;
Geim2005 ; Geim2009 . These modes are also a source for even finer, 1D Dirac
modes, in graphene nanoribbons Yang ; Yang2008 ; Ruffieux  with special
orientation of their edges and special adjustment of their atomic width
Wakabayashi1996 .

A broad family of Dirac semi-metals is of great interest for modern
electronics, in particular, their behavior under doping by different
impurities and the resulting restructuring of quasiparticle spectrum.
Comparing with the known impurity effects in common semiconductors and in 2D
graphene, the doped graphene nanoribbons can be expected to permit even higher
sensitivity to various external controls and their study may deepen our
general knowledge of disorder physics.

This work continues the recent study of impurity effects in graphene
nanoribbons Pogorelov2021 , focusing on their armchair edge orientation and
their width adjusted for presence of Dirac-like modes in the electronic
spectrum. In this course, we study various regimes of spectrum restructuring
under impurity perturbation in function of perturbation parameters and compare
the obtained results with the known such effects in other electronic
materials.

The following consideration begins from the description of an armchair
graphene nanoribbon (AGNR) and its spectral structure in terms of the second
quantization Hamiltonian (Sec. II) and the related Green functions (GFs, Sec.
III). The perturbation of Hamiltonian by impurity adatoms within the Anderson
hybrid model is introduced in Sec. IV, giving the solutions for perturbed GFs
in the T-matrix form and checking for Mott localization of perturbed
quasiparticles. The analysis of possible electronic phase states in doped AGNR
in function of perturbation parameters is developed in Secs. V, VI and
compared with the known behaviors of analogous electronic materials. The final
discussion of the obtained results with suggestions for their possible
practical use is presented in Sec. VII. An additional check of the solutions
for 1D modes, beyond the T-matrix framework, is done in Appendix.

It is our honour to dedicate this paper to 95th birthday of a prominent solid
state theorist Emmanuel Iosifovich Rashba. His brilliant scientific career
began in his native city Kyiv, Ukraine’s capital, where he studied, worked for
a long time and formed as a physicist under an outstanding supervise by S.I.
Pecar. Rashba’s results in non-ideal molecular crystal optics, in
semiconductors theory, and, especially, in spintronics the very emergence of
which was mostly due to his discovery of a new type spin-orbit coupling,
essentially promoted the understanding of processes and phenomena in
electronic systems of various dimensionalities.

##  II Hamiltonian and Green functions

Graphene armchair nanoribbon can be seen as a periodic sequence of $N$
segments where each segment is a slant stack of $M$ layers (collinear between
the segments) and each layer consists of two atomic sites (of graphene $a$\-
and $b$-types, see Fig. 1).

Figure 1: Graphene nanoribbon with armchair orientation of its edges. Blue
dashed lines delimit slant segments $n=1,\dots,N$. Each segment extends along
the graphene elementary translation vector ${\bf a}_{1}$ and includes
$m=1,\dots,M$ layers (red dashed lines) with $a$-type (white) and $b$-type
(black) atomic sites. The sequence of segments has its longitudinal period
$a=|{\bf a}_{1}+{\bf a}_{2}|$. The carbon dangling bonds at the edges are
passivated by hydrogens (small circles).

The respective local electronic states in $m$th layer of $n$th segment are
generated by the local operators $a_{n,m}^{\dagger}$ and $b_{n,m}^{\dagger}$.
Longitudinal translational invariance is imposed through the Born-von Karman
closure of the $N$th to the 1st segment. For an AGNR with $M$ layers
($M$-AGNR) and the nearest neighbor hopping $t$, the related tight-binding
Hamiltonian reads:

$\displaystyle H_{tb}$ $\displaystyle=$ $\displaystyle
t\left\\{\sum_{n=1}^{N}\left[\,\sum_{m=2}^{M}a_{n,m}^{\dagger}\left(b_{n,m}+b_{n-1,m+1}\right.\right.\right.$
(1)
$\displaystyle\qquad\left.\left.+\,\,b_{n,m-1}\right)+a_{n,1}^{\dagger}\left(b_{n,1}+b_{n-1,2}\right)\right]$
$\displaystyle\qquad\qquad+\left.\sum_{m=1}^{M-1}a_{1,m}^{\dagger}b_{N,m+1}+h.c.\right\\}.$

The last sum in Eq. 1 just generates the longitudinal translation invariance
and suggests the Fourier-transform to 1D-wave operators. The longitudinal
coordinates of $a$\- and $b$-sites in units of the longitudinal period
$a=|\mathbf{a}_{1}+\mathbf{a}_{2}|$ define this transform as:

$\displaystyle\alpha_{k,m}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}\sum_{n=1}^{N}{\rm e}^{i(2\pi
k/N)(\xi_{n,m}-1/6)}a_{n,m},$ $\displaystyle\beta_{k,m}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N}}\sum_{n=1}^{N}{\rm e}^{i(2\pi
k/N)(\xi_{n,m}+1/6)}b_{n,m},$ (2)

where $\xi_{n,m}=n+(m+1)/2$ is the longitudinal coordinate of the center of
$m$th layer from $n$th segment (see an example in Fig. 1). This readily
diagonalizes the Hamiltonian, Eq. 1, in the $k$ numbers. If the AGNR is
macroscopically long, $N\to\infty$, one can pass to a quasi-continuous
momentum variable: $2\pi k/N\to k$ (measured in $a^{-1}$ units). Also, for
simplicity, the energy $\varepsilon$ will be measured in units of $t$.

Then the system dynamics in the transversal $m$-index can be considered at a
fixed longitudinal $k$-momentum, and the overall spectrum structure results
qualitatively defined by the AGNR width $M$. In the known analytic approach by
Wakabayashi et al Wakabayashi1996 ; Wakabayashi2009 ; Wakabayashi2010 , $2M$
eigen-states at given $k$ are taken as running $k$-waves superposed by
standing waves in the transversal $q$-momentum, subject to the open edge
condition (reaching a node when continued by a half-period beyond an AGNR
edge). Namely, they are pairs of standing waves with discrete momentum values:

$q_{j}=\frac{\pi j}{M+1},\quad j=1,\dots,M,$ (3)

being just the combinations (symmetric and antisymmetric in $a$\- and
$b$-sites) of 1D-projected graphene states.

The related eigen-energies are simple uniform 1D-projections of the 2D
graphene spectrum for transversal momentum values $q_{j}$ by Eq. 3:

$\varepsilon_{j,k}=\sqrt{1+4\cos\tfrac{k}{2}\cos q_{j}+4\cos^{2}q_{j}},$ (4)

for conduction sub-bands (and $-\varepsilon_{j,k}$ for valence sub-bands). The
1D Brillouin zone (BZ) for all the $2M$ sub-bands is defined within the $0\leq
k\leq 2\pi$ range and the respective secular determinant reads:

$\det(\varepsilon-\hat{H})=\prod_{j=1}^{M}\left(\varepsilon^{2}-{\varepsilon_{j,k}}^{2}\right).$
(5)

The eigen-state associated to the $(j,k)$-mode is a combination of the running
$k$-wave and the standing $q_{j}$-wave Wakabayashi2010  with its amplitudes on
$a$\- and $b$-sites in $m$-layer:

$\displaystyle A^{(j,k)}_{m}=-\frac{{\rm e}^{-i\varphi_{j,k}}}{\sqrt{M+1}}\sin
mq_{j},$ $\displaystyle\qquad\qquad B^{(j,k)}_{m}=\frac{{\rm
e}^{i\varphi_{j,k}}}{\sqrt{M+1}}\sin mq_{j},$ (6)

where the phase is defined by the relation:

$\varphi_{j,k}=\frac{1}{2}\arctan\frac{\sin\tfrac{k}{2}}{\cos\tfrac{k}{2}+2\cos
q_{j}}+\frac{k}{6}.$ (7)

The standing waves are orthonormalized through the relations:

$\displaystyle\sum_{j=1}^{M}\sin mq_{j}\sin m^{\prime}q_{j}$ $\displaystyle=$
$\displaystyle\frac{M+1}{2}\delta_{m,m^{\prime}},$
$\displaystyle\sum_{m=1}^{M}\sin mq_{j}\sin mq_{j^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{M+1}{2}\delta_{j,j^{\prime}}.$ (8)

Then we construct the eigen-mode operators $\psi_{\pm j,k}$ from the wave
operators $\alpha_{m,k}$ and $\beta_{m,k}$ by Eq. 2 in order to reproduce the
mode amplitudes by Eq. 6:

$\displaystyle\psi_{\pm j,k}=\frac{1}{\sqrt{M+1}}\sum_{m=1}^{M}\sin
mq_{j}\left({\rm e}^{i\varphi_{j,k}}\beta_{m,k}\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\mp\left.{\rm
e}^{-i\varphi_{j,k}}\alpha_{m,k}\right).$ (9)

In their basis, the Hamiltonian, Eq. 1, turns fully diagonal:

$H_{tb}=\sum_{j,k}\varepsilon_{j,k}\left(\psi_{j,k}^{\dagger}\psi_{j,k}-\psi_{-j,k}^{\dagger}\psi_{-j,k}\right).$
(10)

By inversion of Eqs. 2, 9, the local operators are expanded in the eigen-mode
operators:

$\displaystyle a_{n,m}=\frac{1}{\sqrt{(M+1)N}}\sum_{j,k}{\rm
e}^{i(k\xi_{n,m}-\varphi_{j,k})}\sin mq_{j}$
$\displaystyle\qquad\qquad\qquad\qquad\times\,\left(\psi_{-j,k}-\psi_{j,k}\right),$
$\displaystyle b_{n,m}=\frac{1}{\sqrt{(M+1)N}}\sum_{j,k}{\rm
e}^{i(k\xi_{n,m}+\varphi_{j,k})}\sin mq_{j}$
$\displaystyle\qquad\qquad\qquad\qquad\times\,\left(\psi_{-j,k}+\psi_{j,k}\right),$
(11)

which is helpful for the next treatment of AGNR perturbations by local
impurity centers.

The notable feature of the spectrum by Eq. 4 is that it contains gapless modes
if the AGNR width satisfies a special condition Wakabayashi1996 :

$M+1=3\nu,\quad\nu=1,2,\dots.$ (12)

For such $M=3\nu-1$, the mode with $j=2\nu$ reaches zero energy at the BZ edge
$k=0$ as:

$\varepsilon_{2\nu,k}=2\left|\sin\frac{k}{4}\right|\approx\frac{|k|}{2}$ (13)

(see Fig. 2), and the mode with $j=\nu$ reaches zero energy at the opposite BZ
edge $k=2\pi$ as:

$\varepsilon_{\nu,k}=2\left|\cos\frac{k}{4}\right|\approx\frac{|k-2\pi|}{2}.$
(14)

Figure 2: Energy bands dispersion in 5-AGNR ($\nu=2$ by Eq. 12), showing the
Dirac-like modes with $j=2,4$.

The dispersion laws by Eqs. 13, 14 formally coincide with the standard linear
dispersion near the Dirac points of 2D graphene, hence they can be seen as
definitions of effective 1D Dirac points in ($3\nu-1$)-AGNR spectra. All other
modes there (with $j\neq\nu,2\nu$) have finite energy gaps.

##  III Green functions and observables

The following consideration goes in the framework of two-time GFs Zubarev1960
; Economou1979  defined by their Fourier transforms:

$\langle\langle
A|B\rangle\rangle_{\varepsilon}=\frac{i}{\pi}\int_{0}^{\infty}dt{\rm
e}^{it(\varepsilon+i0)}\left\langle\left\\{A(t),B(0)\right\\}\right\rangle$
(15)

where $A(t)={\rm e}^{iHt}A{\rm e}^{-iHt}$ is a Heisenberg picture operator for
the system Hamiltonian $H$, $\\{.,.\\}$ is the anti-commutator and
$\langle\dots\rangle$ is the quantum-statistical average. In what follows, the
GF’s energy subindex is mostly omitted.

Practical calculation of GFs is done through the general equation of motion:

$\varepsilon\langle\langle
A|B\rangle\rangle=\langle\left\\{A(0),B(0)\right\\}\rangle+\langle\langle\left[A,H\right]|B\rangle\rangle,$
(16)

involving the commutator $\left[.,.\right]$. So found GFs generate the
physical observable quantities (the averages of operator products) through the
spectral relation:

$\langle AB\rangle=\frac{1}{\pi}{\rm Im}\int_{0}^{\infty}\langle\langle
B|A\rangle\rangle_{\varepsilon}d\varepsilon.$ (17)

In the present case, the system electronic properties can be obtained from the
2M$\times$2M GF matrix $\hat{G}(k,k^{\prime})$ with its matrix elements
$G_{j,j^{\prime}}(k,k^{\prime})\equiv\langle\langle\psi_{j,k}|\psi_{j^{\prime},k^{\prime}}^{\dagger}\rangle\rangle$
built from the eigen-mode operators by Eq. 9 where the $j$-indices count the
transversal momenta as by Eq. 3 and also their opposites $-j$ ($2M$
altogether).

For the unperturbed AGNR system with its diagonal Hamiltonian, Eq. 10, the
above defined GF matrix results also diagonal:
$G_{j,j^{\prime}}^{(0)}(k,k^{\prime})=\delta_{j,j^{\prime}}\delta_{k,k^{\prime}}G_{j,k}^{(0)}(\varepsilon)$,
with its diagonal elements called propagators:

$G_{j,k}^{(0)}(\varepsilon)=\frac{1}{\varepsilon-\varepsilon_{j,k}}.$ (18)

They define an important observable, the density of quasiparticle states
(DOS), as a sum $\rho(\varepsilon)=\sum_{j=1}^{M}\rho_{j}(\varepsilon)$ where
a partial contribution from $(j,-j)$ subbands is:

$\rho_{j}(\varepsilon)=\frac{2}{\pi}\sum_{k}{\rm
Im}\,\left[G_{j,k}^{(0)}(\varepsilon)+G_{-j,k}^{(0)}(\varepsilon)\right]$ (19)

(including the implicit factor 2 for electron spins). Using Eqs. 4, 18 and
passing from sum in $k$ to integral:

$\frac{1}{N}\sum_{k}f(k)=\frac{1}{2\pi}\int_{0}^{2\pi}f(k)dk,$ (20)

gives this contribution as:

$\rho_{j}(\varepsilon)=\frac{8|\varepsilon|}{\pi(M+1)\sqrt{(\varepsilon^{2}-\varepsilon_{-,j}^{2})(\varepsilon_{+,j}^{2}-\varepsilon^{2})}},$
(21)

where $\varepsilon_{\pm,j}=1\pm 2\cos q_{j}$ are the $j$th subband energy
limits.

Figure 3: Density of states in 5-AGNR. Almost constant value $\approx
4/(3\pi)$ within the low energy range, $|\varepsilon|\leq\sqrt{3}-1$, comes
from the Dirac-like ($j=2,4$) subbands.

In particular, for the most relevant here Dirac-like modes with
$\varepsilon_{+,\nu}=\varepsilon_{-,2\nu}=2,\,\varepsilon_{+,2\nu}=\varepsilon_{-,\nu}=0$,
we have:

$\rho_{\nu}(\varepsilon)=\rho_{2\nu}(\varepsilon)=\frac{4}{3\pi\nu\sqrt{1-(\varepsilon/2)^{2}}},$
(22)

so DOS results almost constant at low energies:
$\rho(\varepsilon)\approx\rho(0)\equiv\rho_{0}=8/(3\pi\nu)$ at
$|\varepsilon|\ll 1$, as seen in the case of 5-AGNR presented in Fig. 3.

##  IV Impurity perturbations

The local impurity perturbations in graphene nanoribbons of both zigzag and
armchair types were recently considered within two basic impurity models: the
one-parameter Lifshitz model (LM) Lifshitz_1963 , more adequate for
substitutional impurities, and the two-parameter Anderson hybrid model (AM)
Anderson , for impurity adatoms Pogorelov2021 . For zigzag structures, the
overall conclusion was about their eigen-modes stability (topological
protection) against quasiparticles localization by the impurity disorder, both
in LM and AM models.

However such localization was found in AGNRs with LM impurities, though
reduced in that case to a narrow vicinity of the Dirac point (zero energy).
But this already opens a possibility for Mott’s metal/insulator phase
transitions in a nanosystem and generates the next interest for studying AGNR
behavior under more diversified AM perturbations. In the latter case, a more
complicated intermittence of conducting and localized states in other ranges
of energy spectrum and a broader variety of related phase states for this 1D
system can be expected. Then it would be also of interest to compare such
effects with the known analogs for 3D and 2D electronic systems under impurity
disorder.

The following consideration is focused on special ($3\nu-1$)-AGNRs and
restricted to only their Dirac-like modes. Since these modes with $j=\nu$ and
$2\nu$ give identical and independent contributions to the spectrum, one can
next focus on a single Dirac-like mode, say $j=2\nu$, then denoting
$\psi_{2\nu,k}\equiv\psi_{+,k}$ and $\psi_{-2\nu,k}\equiv\psi_{-,k}$.
Consequently, the above introduced $\hat{G}(k,k^{\prime})$ matrix gets reduced
to the 2$\times$2 form in the basis of $\psi_{\pm,k}$ operators. In
particular, the non-perturbed solution, Eq. 18, presents here as
$\hat{G}^{(0)}(k,k^{\prime})=\delta_{k,k^{\prime}}(\varepsilon-\varepsilon_{k}\hat{\sigma}_{3})^{-1}$
with $\varepsilon_{k}\equiv\varepsilon_{2\nu,k}$ by Eq. 13 and the Pauli
matrix $\hat{\sigma}_{3}$ in $\pm$ indices.

Figure 4: Position $p$ of a top impurity adatom over an $a$-site in
$m_{p}$-layer from $n_{p}$-segment of 5-AGNR.

Then we consider the impurity adatoms location restricted to the simplest top-
positions: $p_{a}$ over an $a$-type host atom in $m_{p}$ layer from $n_{p}$
segment (see Fig. 4) or $p_{b}$ over a $b$-type host atom.

In these notations, the AM perturbation Hamiltonian reads:

$\displaystyle H_{AM}$ $\displaystyle=$
$\displaystyle\sum_{p}\left\\{\varepsilon_{res}c_{p}^{\dagger}c_{p}+\frac{\omega}{\sqrt{3\nu
N}}\sum_{k}\left[\sin\frac{\pi m_{p}}{3}\right.\right.$ (23)
$\displaystyle\times$ $\displaystyle\left.\left.{\rm
e}^{i(k\xi_{p}\mp\phi_{k})}c_{p}^{\dagger}(\psi_{-,k}\mp\psi_{+,k})+{\rm
h.c.}\right]\right\\},$

where a local impurity operator $c_{p}$ with its resonance level
$\varepsilon_{res}$ is coupled to the Dirac-like modes $\pm\varepsilon_{k}$
through the hybridization $\omega$ with the neighbor host atom at the
longitudinal position $\xi_{n_{p},m_{p}}-1/6$ for its $a$-type as in Fig. 4
(or $\xi_{n_{p},m_{p}}+1/6$ for its $b$-type).

Then the complete Hamiltonian $H_{tb}+H_{AM}$ generates a perturbation of the
GF matrix: $\hat{G}^{(0)}\to\hat{G}$. In its simplest form, this is given by
the T-matrix approximation:

$\hat{G}=\left(\varepsilon-c\hat{T}(\varepsilon)-\varepsilon_{k}\hat{\sigma}_{3}\right)^{-1},$
(24)

where $c=(2MN)^{-1}\sum_{p}1$ is the impurity concentration and the T-matrix
in this case results diagonal:

$\hat{T}(\varepsilon)\equiv
T(\varepsilon)=\frac{\omega^{2}}{2}\left(\varepsilon-\varepsilon_{res}-\frac{i\omega^{2}}{4f\sqrt{1-(\varepsilon/2)^{2}}}\right)^{-1}.$
(25)

Figure 5: Dispersion law for the Dirac-like modes in 5-AGNR having AM
impurities with their parameters $\varepsilon_{res}=0.03$, $\omega=0.3$, and
concentration $c=0.02$. For comparison, the unperturbed Dirac-like modes are
shown by the red dashed lines.

Then the modified dispersion law $\tilde{\varepsilon}_{k}$ follows from the
standard GF secular equation ${\rm Re}\,[\det(\hat{G})^{-1}]=0$ Bonch  as:

$\tilde{\varepsilon}_{k}=\sqrt{\varepsilon_{k}^{2}+(c\,{\rm
Im}\,T(\varepsilon))^{2}}+c\,{\rm Re}\,T(\varepsilon),$ (26)

and its solution for 5-AGNR at the choice of AM impurity parameters
$\varepsilon_{0}=0.03$ and $\omega=0.3$, corresponding to Cu adatoms in top
position Pogorelov2020 , is presented in Fig. 5. Here the characteristic
impurity effects are seen as:

1) the shift of the Dirac point from its initial zero energy position down to

$\varepsilon_{\rm
D}\equiv\tilde{\varepsilon}_{k=0}\approx\frac{\varepsilon_{res}-\sqrt{\varepsilon_{res}^{2}+2c\omega^{2}}}{2},$
(27)

2) the resonance splitting between the initial linear $\varepsilon_{k}$ law
and the impurity resonance level $\varepsilon_{res}$ until its vicinity of
width

$\gamma_{res}\approx\omega\sqrt{\frac{c-c_{0}}{2}}.$ (28)

where $c_{0}\sim\omega^{2}/(8\nu^{2})$ is the critical concentration value for
this splitting to appear.

Also an anomalous negative dispersion formally appears inside this vicinity,
at $|\varepsilon-\varepsilon_{res}|\lesssim\gamma_{res}$, but this range
occurs unphysical when validity of the modified dispersion law is checked with
the Ioffe-Regel-Mott (IRM) criterion for conducting states IoffeRegel ; Mott :

$k\frac{\partial\tilde{\varepsilon}_{k}}{\partial k}\gtrsim c\,{\rm
Im}\,T(\tilde{\varepsilon}_{k}).$ (29)

This simply means that the quasiparticle lifetime (inverse of the r.h.s) is
longer than its oscillation period (inverse of the l.h.s) and such
quasiparticles are conductive indeed, otherwise they are localized near
impurity sites. The explicit form of the IRM criterion for the dispersion law
by Eq. 26 and the T-matrix by Eq. 25 reads:

$\frac{2}{\partial\ln{\rm Re}\,\left[\varepsilon-
cT(\varepsilon)\right]^{2}/\partial\varepsilon}\gtrsim c\,{\rm Im}\,T,$ (30)

All the energy ranges where this inequality does not hold are attributed to
localized states so the dispersion law by Eq. 26 for conducting states does
not apply there. The mobility edges between conducting and localized states
can be qualitatively estimated as the $\varepsilon$ values that make the
relation of Eq. 30 an equality. The results of such numerical estimation at
the choice of AM parameters as for Cu top impurities are shown in Fig. 6a.
They illustrate formation of two mobility gaps (ranges of localization), one
around $\varepsilon_{res}$ and another around $\varepsilon_{D}$. Their width
grows with the impurity concentration $c$: the first as $\sim
2\gamma_{res}(c)$ (by Eq. 28) and the second as $\gamma_{D}\approx
c\omega^{4}/[8\nu(\varepsilon_{D}-\varepsilon_{res})^{2}]$.

Figure 6: a) Mobility gaps (in meV scale) in function of impurity
concentration $c$ at the choice of AM parameters as in Fig. 5. Special
concentration values refer respectively: $c_{0}$ to opening of mobility gap
near the resonance level $\varepsilon_{res}$, $c_{1}$ to metal/insulator, and
$c_{2}$ to insulator/metal transitions (see below). b) Analogous developments
for the same impurities in 2D graphene, note here opening of a spectrum quasi-
gap inside the mobility gap at $c_{1}$, absent in AGNR.

Validity of the above T-matrix description is confirmed with an additional
check beyond the frame of this simplest approximation (see in Appendix).

An important issue for such disordered AGNRs is how the positions of mobility
edges compare with that of the Fermi level vs impurity concentration,
$\varepsilon_{\rm F}(c)$. The latter results from an extra amount of $c$
charge carriers by impurities (per host atom) being filled into the relevant
subbands:

$\int_{0}^{\varepsilon_{\rm F}(c)}\rho(\varepsilon)d\varepsilon=c.$ (31)

Notably, the above indicated low energy DOS,
$\rho(\varepsilon)\approx\rho_{0}$, holds its constancy even under impurity
disorder, as seen from stability of Eq. 22 at the passing
$\varepsilon\to\varepsilon-cT(\varepsilon)$ since $|\varepsilon-
cT(\varepsilon)|\ll 1$ for all $|\varepsilon|\ll 1$, one of the main specifics
of this 1D-like system.

Then from Eq. 31 we come to a simple linear relation $\varepsilon_{\rm
F}(c)\approx c/\rho_{0}$. Such behavior superposed onto the diagram of
mobility gaps in Fig. 6a shows a possibility for $\varepsilon_{\rm F}$ to
cross the mobility edges near $\varepsilon_{res}$, both into and out of this
mobility gap. Then an intermittency of the related metal/insulator transitions
(MIT) can be expected.

It is of interest to compare these AGNR results with their analogs for the
same impurities in 2D graphene (Fig. 6b). Besides a general similarity of two
pictures, they also present substantial differences. First of all, the
estimated critical concentration for localization on Cu impurities in 5-AGNR,
$c_{0}\approx 7\cdot 10^{-3}$, is more than an order of magnitude higher than
its analog for the 2D case, $c_{0}^{(2D)}\approx 4\cdot 10^{-4}$ Pogorelov2020
. This can be explained by the higher and almost constant low-energy DOS for
the 1D Dirac-like modes: $\rho_{0}\sim 1$, compared to its linear in
$\varepsilon$ smallness for the 2D Dirac modes: $\rho_{2D}(\varepsilon)\approx
4\varepsilon/(\pi\sqrt{3})\ll 1$.

Another difference is in behaviors of $\varepsilon_{\rm F}(c)$ for each
system, also caused by that of DOS. From Eq. 31, its almost linear growth in
AGNR: $\varepsilon_{\rm F}(c)\approx c/\rho_{0}$ is much slower than the very
fast, square root initial growth of $\varepsilon_{\rm
F}^{(2D)}(c)\approx\sqrt{c\pi\sqrt{3}/2}$ in 2D graphene. This defines a much
higher threshold in $c$ for occurrence of MIT in AGNR:
$c_{1}\sim\varepsilon_{res}\rho_{0}$, than in 2D graphene: $\sim
c_{0}^{(2D)}$. But the constancy of AGNR DOS, even in presence of impurities,
permits the linear $\varepsilon_{\rm F}(c)$ growth to persist also for
$c>c_{1}$ and then an inverse MIT to occur at its emergence from the mobility
gap at $c=c_{2}\sim\varepsilon_{res}\rho_{0}+(\omega\rho_{0})^{2}$. In
contrary, the 2D graphene DOS presents a sharp resonance peak near
$\varepsilon_{res}$ and the Fermi level, when approaching this peak at $c\to
c_{0}^{(2D)}$, gets fixed here, defining insulating phase for all
$c>c_{0}^{(2D)}$.

At least, the critical concentration $c_{0}$ by Eq. 28 that defines the onset
of localization near $\varepsilon_{res}$, indicates it to occur earlier for
weaker impurity-host coupling $\omega$.

##  V Electronic phase states and their tuning

An important physical issue is to determine the system electronic phase
states. For the considered nanoribbons, this refers first of all to their
electric conductivity.

In the limit of zero temperature, it is fully defined by the Fermi level
position with respect to the spectrum mobility edges: implying the metallic
state for $\varepsilon_{\rm F}$ out of the mobility gaps and the insulating
state for $\varepsilon_{\rm F}$ inside them. Thus, from the diagram in Fig. 6,
the system of 5-AGNR with given concentration $c$ of Cu impurities is expected
to be metallic if $0<c<c_{1}$ or $c>c_{2}$ and insulating if $c_{1}<c<c_{2}$.

But a practical interest arises in an effective tuning of possible MIT’s at a
given impurity composition (in analogy with the common gate controls in doped
semiconductors). First of all, the initial composition can be chosen to set
the Fermi level close enough to a mobility edge, for instance,
$\varepsilon_{\rm F}<\varepsilon_{g}<\varepsilon_{res}$ and $\varepsilon_{\rm
F}-\varepsilon_{g}\ll\varepsilon_{res}$. Then several factors can be
considered for MIT tuning: 1) temperature (in terms of its inverse $\beta$),
2) magnetic field, $\bf h$, and 3) electric field, $\bf E$.

The temperature control will result from the interplay between the metallic
conductivity $\sigma_{met}(\beta)$, due to extended states, and the hopping
conductivity $\sigma_{hop}(\beta)$, due to localized states.

The first type refers to the Kubo-Greenwood formula written here as:

$\sigma_{met}(\beta)\approx\rho_{0}\int_{0}^{\varepsilon_{g}}[1-(\varepsilon/2)^{2}]\tau(\varepsilon)\frac{\partial
n(\varepsilon,\beta)}{\partial\varepsilon}d\varepsilon,$ (32)

where $n(\varepsilon,\beta)=[{\rm e}^{\beta(\varepsilon-\varepsilon_{\rm
F})}+1]^{-1}$ is the standard Fermi function and the lifetime
$\tau(\varepsilon)=[\tau_{imp}^{-1}(\varepsilon)+\tau_{ph}^{-1}(\varepsilon)]^{-1}$
involves:

$\tau_{imp}^{-1}(\varepsilon)=c|{\rm Im}\,T(\varepsilon)|$
($\beta$-independent) due to impurities, and

$\tau_{ph}^{-1}(\varepsilon)\sim 1/\beta\Theta_{\rm D}$, due to 1D phonons
(with the Debye temperature $\Theta_{\rm D}$).

The second type refers to the Mott formula written for an 1D system:

$\sigma_{hop}(\beta)\propto{\rm e}^{-\sqrt{\beta T_{0}}},$ (33)

where $T_{0}\sim\tau^{-1}(\varepsilon_{\rm F})/\rho(\varepsilon_{\rm F})$.
Interplay of this decreasing $\sigma_{hop}(\beta)$ and the growing
$\sigma_{met}(\beta)$ by Eq. 32 results in an overall conductivity maximum at
$\beta\sim 1/|\varepsilon_{\rm F}-\varepsilon_{g}|$, but having a comparable
temperature width. So, this crossover between the types of conductivity is not
yet a canonical phase transition.

Figure 7: Realization of MIT in AGNR with Cu top impurities by tuning of the
composed mobility edges: a) with an applied magnetic field or b) with an
applied electric field (here for $M=5$).

But a true electronic phase transition at zero temperature can be realized,
for instance, applying a uniform static magnetic field $\bf h$ to AGNR. This
will produce a spin splitting of the Dirac-like subbands defined in Sec. III
and also of the impurity levels, implying respective splitting of IRM critical
points for spin subbands.

At the same time, the position of overall Fermi level for a given impurity
concentration $c$ will stay the same as it was for $h=0$, due to the
persisting constancy of the overall DOS. Then, in the situation of several
overlapping subbands, the overall mobility edges are determined by the Mott
principle: if, at a given energy, there is at least one state extended, all
other states at this energy are also extended. Then the overall mobility gap
is formed by the intersection of partial (formal) gaps for each spin
projection and it gets reduced with growing splitting $\mu_{\rm B}h$. In this
way, the overall mobility edges are tuned by the applied field and MIT is
realized at its critical value $h_{cr}\sim|\varepsilon_{\rm
F}-\varepsilon_{g}|/\mu_{\rm B}$ (see Fig. 7a). But for the relevant energy
scales of several meV, this may require high enough $h$ values of several tens
Tesla.

An alternative way may be sought in applying a static electric field
$E_{y}\equiv E$ across the nanoribbon (along the $y$-axis in Figs. 1, 4) to
produce linearly growing local potentials $V_{m}=meE$ on $m$-layers. This can
be shown not to influence the relevant Dirac-like subbands, however to produce
an $M$-fold splitting between the local energy levels for impurities on
different $m$-layers and so between the respective mobility gaps, with a
subsequent decrease of the resulting mobility gap (as in Fig. 7b). The
critical fields to achieve MIT in this case, $E_{cr}\sim|\varepsilon_{\rm
F}-\varepsilon_{g}|/(eM)$ of some mV/nm, could be reached with no experimental
difficulties.

##  VI Spin-orbit effects on electronic phase states

Yet one more tuning mechanism can result from the spin-orbit (SO) effect,
including the Rashba spin-orbit coupling Bychkov . The latter is generally
known to lift the spin degeneracy in the systems with broken mirror symmetry,
for instance, under electric field $E_{z}$, applied normally to the crystal
surface Vas'ko  or to the 2D graphene plane (along the $z$-axis in Fig. 4)
Kane ; Kuemmeth . A similar effect can be achieved in graphene nanotubes
Huertas-Hernando ; Rudner  and also in AGNRs Lenz .

Figure 8: Coupled pair of $m=1,2$ layers in AGNR.

For the AGNR case, we note that the relevant Dirac-like $\psi_{\pm,k}$ modes
(from Eq. 23) have their amplitudes $\sin m\tfrac{\pi}{3}$ equal zero at each
third $m$-layer, so the interlayer couplings through these modes are
restricted to the $m$-pairs: $(1,2),(4,5),\dots(M-1,M)$, and their overall
effect on AGNR can be represented by a single pair, for instance, with $m=1,2$
(see Fig. 8). All the couplings in such a pair of layers are suitably
presented in terms of local operators, now equipped with explicit
$\uparrow\downarrow$ spin indices and then composed into 4-spinors:

$f_{n,m}=\left(\begin{array}[]{c}a_{n,m,\uparrow}\\\ a_{n,m,\downarrow}\\\
b_{n,m,\uparrow}\\\ b_{n,m,\downarrow}\\\ \end{array}\right).$ (34)

In the basis of these local spinors, the SO Hamiltonian reads:

$\displaystyle H_{SO}$ $\displaystyle=$
$\displaystyle\sum_{n}\left[f_{n,1}^{\dagger}\hat{H}_{SO}\left(f_{n,1}+f_{n,2}+f_{n-1,2}\right)\right.$
(35) $\displaystyle+$
$\displaystyle\left.f_{n,2}^{\dagger}\hat{H}_{SO}\left(f_{n,2}+f_{n,1}+f_{n+1,1}\right)\right],$

where the 4$\times$4 matrix:

$\hat{H}_{SO}=\Delta\hat{\sigma}_{z}+\lambda\left(\hat{\sigma}_{x}\hat{\tau}_{y}-\hat{\sigma}_{y}\hat{\tau}_{x}\right).$

includes the Pauli matrices $\hat{\sigma}_{j}$ in spin
$\uparrow,\downarrow$-indices and $\hat{\tau}_{j}$ in sublattice $a,b$-indices
and also the parameters $\Delta$ for standard SO and $\lambda$ for Rashba SO
(the latter being $E_{z}$-dependent). The estimates for these local SO
couplings in 2D graphene (also plausible for nanoribbons) show the standard
$\Delta\sim 10^{-4}$ Kane , fixed and much smaller of the relevant energy
scales for MIT crossing. Otherwise, the Rashba $\lambda$ can be strongly
enhanced Marchenko  and yet tunable Kane , so it is taken as an effective SO
variable below.

Next, we pass to the basis of chain-wave 4-spinors:

$\psi_{k}=\left(\begin{array}[]{c}\psi_{+,k,\uparrow}\\\
\psi_{+,k,\downarrow}\\\ \psi_{-,k,\uparrow}\\\ \psi_{-,k\downarrow}\\\
\end{array}\right),$ (36)

which are related to the local spinors by Eq. 34 (with $m=1,2$) through a
$\hat{\tau}$-rotation:

$f_{n,m}=\frac{(-1)^{m-1}}{2\sqrt{\nu N}}\sum_{k}{\rm
e}^{ik\xi_{n,m}}\hat{U}_{k}\psi_{k}.$ (37)

Here the rotation matrix:

$\hat{U}_{k}=\cos\phi_{k}\left(\hat{\tau}_{x}-\hat{\tau}_{z}\right)+\sin\phi_{k}\left(\hat{\tau}_{y}+i\hat{\tau}_{0}\right),$
(38)

results from Eq. 11 (restricted to single $j=2\nu$) for the components of
$f_{n,m}$. At relevant $|k|\ll 1$, the phase $\phi_{k}\equiv\phi_{2f,k}$ is
approximated as $\phi_{k}\approx\phi_{0}=-\pi/4$. Then the Dirac-like part of
SO Hamiltonian reads:

$H_{SO}=\sum_{k}\psi_{k}^{\dagger}\left[-\Delta\hat{\sigma}_{z}\hat{\tau}_{x}+\lambda\left(\hat{\sigma}_{x}\hat{\tau}_{z}+\hat{\sigma}_{y}\hat{\tau}_{y}\right)\right]\psi_{k},$
(39)

and, together with the Dirac-like part of $H_{tb}$ by Eq. 10, it defines the
SO-split dispersion laws (in neglect of impurity disorder):

$\varepsilon_{\pm,k}=\sqrt{\varepsilon_{k}^{2}+\Delta_{\pm}^{2}}.$ (40)

Figure 9: SO splitting of the Dirac-like mode $\varepsilon_{k}$ (dashed lines)
into $\varepsilon_{\pm,k}$ at the choice of $\lambda=5\cdot 10^{-3}$ and
$\Delta=2.5\cdot 10^{-4}$.

Here the non-zero bandgaps:

$\Delta_{\pm}=\sqrt{\Delta^{2}+\lambda^{2}\left(2\pm\sqrt{3}\right)},$ (41)

are due to both SO types but their splitting is only due to the Rashba SO as
shown in Fig. 9 for the choice of SO parameters $\Delta=2.5\cdot 10^{-4}$ Kane
and $\lambda=5\cdot 10^{-3}$ Marchenko , and they create the low-energy DOS
singularities:

Figure 10: Singularities of the low-energy DOS of 5-AGNR with split subbands
at the same choice of SO parameters as in Fig. 9.

$\rho_{so}(\varepsilon)=\frac{\varepsilon\rho_{0}}{2}\left[\frac{\theta(\varepsilon-\Delta_{+})}{\sqrt{\varepsilon^{2}-\Delta_{+}^{2}}}+\frac{\theta(\varepsilon-\Delta_{-})}{\sqrt{\varepsilon^{2}-\Delta_{-}^{2}}}\right],$
(42)

as shown in Fig. 10.

Such DOS behavior, instead of its almost constancy at no account of SO (by Eq.
22), when used in Eq. 31 leads to the $\lambda$-dependence of the Fermi level
given by the equation:

$\sqrt{\varepsilon_{\rm F}^{2}-\Delta_{+}^{2}}+\sqrt{\varepsilon_{\rm
F}^{2}-\Delta_{-}^{2}}\approx 2c/\rho_{0}.$ (43)

Its numerical solution defines MIT to occur when $\lambda$ reaches a certain
critical value $\lambda_{cr}$ as shown in Fig. 11 .

Noting that for all relevant impurity concentrations $c>c_{0}$ we have
$c/\rho_{0}\gg\Delta$, the approximate solution of Eq. 43 for $\lambda\ll
c/\rho(0)$ reads:

$\varepsilon_{\rm
F}(\lambda)\approx\frac{c}{\rho_{0}}+\frac{\rho_{0}}{c}\left(\Delta^{2}+2\lambda^{2}\right).$
(44)

Figure 11: Dependence of the Fermi level on Rashba SO coupling for 5-AGNR with
Cu top impurities concentration fixed at $c=0.011$, its crossing with the
mobility edge $\varepsilon_{g}$ at $\lambda=\lambda_{cr}\approx 0.0078$
indicates an SO-tuned MIT to occur.

This quadratic $\lambda$-dependence relates to filling of both
$\varepsilon_{\pm,k}$ subbands by $c$ charge carriers, it follows the lowest
part of the numerical solution in Fig. 11. With further growing $\lambda$, a
break of $\varepsilon_{\rm F}(\lambda)$ occurs when it gets intercepted by the
faster growing upper bandgap $\varepsilon_{+}$ at the break value
$\lambda_{br}=c/(\sqrt{2\sqrt{3}}\rho_{0})$. The next slower $\varepsilon_{\rm
F}(\lambda)$ growth at $\lambda>\lambda_{br}$ relates to filling of only the
lower subband which is expressed as:

$\varepsilon_{\rm
F}(\lambda)=\sqrt{(c/\rho_{0})^{2}+\Delta^{2}+(2-\sqrt{3})\lambda^{2}}.$ (45)

Notably, for the considered case of Cu top impurities, both the impurity
resonance $\varepsilon_{res}$ and the mobility edge $\varepsilon_{g}$ lie in
the energy range $\varepsilon\gg\varepsilon_{\pm}$ where
$\rho_{SO}(\varepsilon)$ already reaches its asymptote $\approx\rho(0)$, so
the mobility gap structure stays practically insensitive to Rashba SO. Hence a
possibility arises here for MIT to be realized by SO tuning of
$\varepsilon_{\rm F}$ at fixed $\varepsilon_{g}$, unlike the above considered
regimes with tuning of mobility edges at fixed Fermi level.

Notably, this tuning process can be realized in a combined way: a rough
"tuning" of $\varepsilon_{\rm F}$ closeness to $\varepsilon_{g}$ by a proper
choice of impurity parameters $\varepsilon_{res},\,\omega,\,c$ and also by a
strong structural contribution to the Rashba parameter $\lambda$, say, from a
gold substrate atomic field Marchenko , and then its fine tuning by an applied
external field. Evidently, the expected MIT at such sub-meV energy scales
would require a range of liquid He temperatures for its sufficient resolution.

##  VII Discussion

The obtained results demonstrate how the difference of electronic states in
graphene nanoribbons defined by their edge orientations is reflected in their
stability against impurity disorder. Physically, this opens the possibility
for specific electronic phase transitions and for their controls by combining
the disorder and external bias effects.

It is of interest to compare the above picture of 1D spectrum reconstruction
under AM impurity disorder with the known such effects in 3D and 2D systems.
This comparison can be done between the corresponding correlators, defined by
the system dimensionality and low energy quasiparticles dispersion. For
instance, such a correlator in 3D semiconductors (quadratic dispersion)
Ivanov1979 , decays with distance $r$ as:
$\propto(\sin\sqrt{\varepsilon}r)/(\sqrt{\varepsilon}r)$, a similar behavior
for 3D acoustic phonons (linear dispersion) Ivanov1970  is found as:
$\propto(\sin\varepsilon r)/(\varepsilon r)$, and it is modified to
$\propto(\sin\varepsilon r)/\sqrt{\varepsilon r}$ in 2D graphene with linear
Dirac dispersion Pogorelov2020 . All these systems admit both extended and
localized quasiparticles based on impurity states ILP . But, unlike all those,
a constancy or only a weak exponential decay of correlators are found for the
present AGNR case (see Eqs. 47, 51 below) which defines the complete
localization of all the states near $\varepsilon_{res}$ (as in Fig. 6a, unlike
the 2D graphene case in Fig. 6b). Another AGNR specifics, the low-energy DOS
constancy, defines a higher sensitivity of Fermi level to doping and so more
possibilities for tuning of the system electronic properties.

The present study was limited to the simplest framework of tight-binding model
for pure nanoribbons and simplest models for impurity perturbations on them.
In principle, it can be extended to account for many other physical factors,
as electron-electron Hubbard correlations, spin-ordering effects, phonon and
spinon excitations, etc.

Also, the effects from passivating hydrogens, known to be commonly present at
the edges of experimental nanoribbon samples Wang2016 ; Ruffieux2016 ;
Wang2017 ; El Abbassi , may influence the dynamics of host nanoribbon
carriers. This factor can be naturally included into the above developed
Hamiltonians and resulting GFs, to be possibly an object of future study. At
least, an experimental check for the suggested effects, for instance, on
carriers mobility and its collapse under definite external factors should be
of considerable interest.

## Acknowledgments

We are grateful to Larissa S. Brizhik, Mark I. Dykman and Aleksandr A. Eremko
for useful discussions of this work. The work by V.M.L. was partially
supported by Grant Nos. 0117U000236 and 0117U000240 from the Department of
Physics and Astronomy of the National Academy of Sciences of Ukraine.

## References

  * (1) K.S. Novoselov, A.K. Geim, S.V. Morozov, S.V. Dubonos, Y. Zhang, D. Jiang, arXiv:cond-mat/0410631 (2004).
  * (2) S.V. Morozov, K.S. Novoselov, F. Schedin, D. Jiang, A.A. Firsov, A.K. Geim, Phys. Rev. B, 72, 201401 (2005).
  * (3) A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. Novoselov, A.K. Geim, Rev. Mod. Phys., 109, 109 (2009).
  * (4) L. Yang, C.-H. Park, Y.-W. Son, M.L. Cohen, S.G. Louie, Phys. Rev. Lett., 99, 186801 (2007).
  * (5) X. Yang, X. Dou, A. Rouhanipour, L. Zhi, H.J. Räder, K. Müllen, J. Amer. Chem. Soc., 130, 4216(2008).
  * (6) P. Ruffieux, J. Cai, N.C. Plumb, L. Pattney, D. Prezzi, A. Ferretti, E. Molinari, X. Feng, K. Müllen, C.A. Pignedoli, R. Fasel, ACS Nano, 6, 6930 (2012).
  * (7) M. Fujita, K. Wakabayashi, K. Nakada, K. Kusakabe, J. Phys. Soc. Jpn, 65, 1920 (1996).
  * (8) Y.G. Pogorelov, D. Kochan, V.M. Loktev, Low Temp. Phys., 47, 754 (2021).
  * (9) K. Wakabayashi, Y. Takane, M. Yamamoto, M. Sigrist, New Journal of Physics, 11, 095016 (2009).
  * (10) K. Wakabayashi, K. Sasaki, T. Nakanishi, T. Enoki, Sci. Technol. Adv. Mater., 11, 054504 (2010).
  * (11) D. Zubarev, Sov. Phys. Uspekhi, 3, 320 (1960).
  * (12) E.N. Economou, Green’s Functions in Quantum Physics, Springer (1979).
  * (13) I. M. Lifshitz, Soviet Physics JETP, 17, 1159 (1963).
  * (14) P. W. Anderson, Phys. Rev., 124, 41 (1961).
  * (15) V.L. Bonch-Bruevich and S.V. Tyablikov, The Green Function Method in Statistical Mechanics, North Holland Publishing House (1962).
  * (16) Y.G. Pogorelov, V.M. Loktev, D. Kochan, Phys. Rev. B, 102, 155414 (2020).
  * (17) A.F. Ioffe and A.R Regel, Progress Semicond., 2, 237 (1960).
  * (18) N.F. Mott, Phil. Mag., 102, 835 (1962).
  * (19) Yu. A. Bychkov and E. I. Rashba, Sov. Phys. - JETP Lett., 39, 78 (1984).
  * (20) F. T. Vas’ko, Sov. Phys. - JETP Lett., 30, 541 (1979).
  * (21) C.L. Kane and E.J. Mele, Phys. Rev. Lett., 95, 226801 (2004).
  * (22) Ferdinand Kuemmeth and Emmanuel I. Rashba, Phys. Rev. B, 80, 241409 (2009).
  * (23) D. Huertas-Hernando, F. Guinea and A. Braatas, Phys. Rev. B, 74, 155426 (2006).
  * (24) Mark S. Rudner and Emmanuel I. Rashba, Phys. Rev. B, 81, 125426 (2010).
  * (25) Lucia Lenz, Daniel F. Urban and Dario Bercioux, Eur. Phys. J B, 86, 502 (2013).
  * (26) D. Marchenko, A. Varykhalov, M.R. Scholz, G. Bihlmayer, E.I. Rashba, A.M. Shikin and O. Rader, Nat. Commun, 3, 1 (2012).
  * (27) S. Wang, L. Talirz, C.A. Pignedoli, X. Feng, K. Müllen, V. Meunier, R. Fasel and P. Ruffieux, Nat. Commun., 7, 11507 (2016).
  * (28) P. Ruffieux, S. Wang, B. Yang, C. Sanchez-Sánchez, J. Liu, T. Dienel, L. Talirz, P. Shinde, C.A. Pignedoli, D. Passerone, T. Dumslaff, X. Feng, K. Müllen and R. Fasel, Nature, 531, 489 (2016).
  * (29) S. Wang, N. Kharche, E. Costa Girão, X. Feng, K. Müllen, V. Meunier, R. Fasel, P. Ruffieux, Nano Lett., 17, 4277 (2017).
  * (30) M. El Abbassi, M.L. Perrin, G. Borin Barin, S. Sangtarash, J. Overbeck, O. Braun, C.J. Lambert, Q. Sun, T. Prechtl, A. Narita, K. Müllen and P. Ruffieux, H. Sadeghi, R. Fasel and M. Calame, ACS Nano, 14, 5754 (2020).
  * (31) I. M. Lifshits, S. A. Gredeskul, L. A. Pastur, Introduction to the theory of disordered systems, Wiley, New York (1988).
  * (32) M.A. Ivanov , V.M. Loktev and Yu.G. Pogorelov, Physics Reports, 153, 209 (1987).
  * (33) R. J. Elliott, J. A. Krumhansl and P. L. Leath, Rev. Mod. Phys., 46, 465 (1974).
  * (34) M.A. Ivanov and Yu.G. Pogorelov, Sov. Phys. JETP, 79, 1010 (1979).
  * (35) M.A. Ivanov, Sov. Phys. Solid State, 12, 1508 (1970).

##  Appendix A Beyond the T-matrix

Validity of the above T-matrix solution should be yet verified in view of the
effective 1D character of the relevant Dirac-like quasiparticles. It is known
that generic 1D systems are unstable against any disorder, up to a full
localization of all their eigen-states Lifshitz , and the IRM check with use
of the simplest (single-impurity) T-matrix is not sufficient to detect this.
Therefore the IRM results from Sec. III need a support by some T-matrix
extensions known for many other disordered systems. There are two such
possible extensions:

1) group expansions (GEs) in clusters of correlated impurity centers ILP  and

2) self-consistent T-matrix approximation Elliott .

For GEs, their basic elements are the correlators, defined for different GE
forms. The simplest form is that of non-renormalized GE, known to better apply
for the energy ranges of localized states. Here the correlator for the
considered Dirac-like quasiparticles is written as:

$A_{r}(\varepsilon)=\frac{2T(\varepsilon)}{3\nu N}\sum_{k}{\rm
e}^{ikr}\left(\frac{1}{\varepsilon-\varepsilon_{k}}+\frac{1}{\varepsilon+\varepsilon_{k}}\right)$
(46)

(taking into account equal contributions from $j=\nu$ and $j=2\nu$ modes).
Then, after integration in $k$ by Eq. 20, the related integral has its long
distance asymptotics at $r\gg 1$ as:

$A_{r}(\varepsilon)=\frac{2T(\varepsilon)\varepsilon}{3\pi\nu}\int_{-\pi}^{\pi}\frac{{\rm
e}^{ikr}dk}{\varepsilon^{2}-4\sin^{2}k/2}\approx\frac{2T(\varepsilon)\varepsilon}{3\nu}\sin\varepsilon
r,$ (47)

that is non-decaying. This contrasts with the decaying correlators in 3D and
2D systems and makes all the GE terms for the 1D-like system formally
divergent. To avoid that problem, some alternative, renormalized GE forms
(more adequate for conductive states) could be employed.

For instance, the first order renormalization for GE is obtained with the
simple change, $\varepsilon\to\varepsilon-cT(\varepsilon)$ in the denominators
of Eq. 46. Then the renormalized correlator
$\tilde{A}_{r}(\varepsilon)=2T(\varepsilon)[\varepsilon-
cT(\varepsilon)]\tilde{I}_{r}(\varepsilon)/(3\nu)$, involves the integral:

$\tilde{I}_{r}(\varepsilon)=\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{{\rm
e}^{ikr}dk}{[\varepsilon-cT(\varepsilon)]^{2}-4\sin^{2}k/2}.$ (48)

Figure 12: Integration contours for calculation of the renormalized correlator
$\tilde{A}_{r}(\varepsilon)$, Eq. 48: $C_{+}$ (blue lines) for $r>0$ and
$C_{-}$ (red lines) for $r<0$.

This can be found analytically, passing to complex momentum:
$k\to\zeta=k+ik^{\prime}$ and extending integration to one of the closed
contours shown in Fig. 8, depending on the 1D correlator direction. The
forward direction, $r>0$, relates to $C_{+}$ (blue lines) with the pole
$\zeta_{\varepsilon}=2\arcsin[\varepsilon-cT(\varepsilon)]/2$, and the
backward direction, $r<0$, does to $C_{-}$ (red lines) with the pole
$-\zeta_{\varepsilon}$.

The contour integral for the forward case:

$\frac{1}{2\pi}\oint_{C_{+}}\frac{{\rm e}^{i\zeta r}}{\varepsilon-
cT(\varepsilon)-2\sin\zeta/2}d\xi,$ (49)

presents a zero sum of three terms:
$0=I_{r}(\varepsilon)+R_{r}(\varepsilon)+V_{r}(\varepsilon)$. Here the residue
term, $R_{r}(\varepsilon)={\rm
e}^{i\zeta_{\varepsilon}r}/\cos(\zeta_{\varepsilon}/2)$, and the term from the
semi-infinite vertical segments, $V_{r}(\varepsilon)\approx(2i\sin\pi
r)/[(1-\varepsilon)r]$, define the sought correlator as:

$\displaystyle\tilde{A}_{r}(\varepsilon)=\frac{2T(\varepsilon)[\varepsilon-
cT(\varepsilon)]}{3f}$ $\displaystyle\qquad\qquad\times\,\left[\frac{{\rm
e}^{i\zeta_{\varepsilon}r}}{\cos(\zeta_{\varepsilon}/2)}+2i\frac{\sin\pi
r}{(1-\varepsilon)r}\right).$ (50)

For the low energy range, $|\varepsilon|,c|T(\varepsilon)|\ll 1$ and
$\zeta_{\varepsilon}\approx\varepsilon-cT(\varepsilon)$, Eq. 50 simplifies to:

$\displaystyle\tilde{A}_{r}(\varepsilon)\approx\frac{2T(\varepsilon)[\varepsilon-
cT(\varepsilon)]}{3\nu}$ $\displaystyle\times\,\left[{\rm
e}^{i[\varepsilon-c{\rm Re}\,T(\varepsilon)]r}{\rm e}^{-c|{\rm
Im}\,T(\varepsilon)r|}+2i\frac{\sin\pi r}{r}\right].$ (51)

Here, unlike Eq. 47 both terms in the brackets are already decaying with $r$.
For the backward case, integration over $C_{-}$ gives the same result.

The relevant criterion for GE convergence, the smallness of the dominating
contribution by impurity pair clusters into the quasiparticle self-energy
compared to that by single impurities ILP , is presented here in the form:

$B_{2}(\varepsilon)\approx
c\left|\int_{0}^{\infty}\tilde{A}_{-r}(\varepsilon)\tilde{A}_{r}(\varepsilon)^{2}{\rm
e}^{-i(\varepsilon-c{\rm Re}\,T(\varepsilon))r}dr\right|\ll 1.$ (52)

Figure 13: Smallness of the relative contribution to GE by impurity pairs,
$B_{2}(\varepsilon)$ by Eq. 53 (for Cu impurities with concentration $c=0.05$
in 5-AGNR), assuring GE convergence for this system.

Then the residue term in Eq. 51 with its slower exponential decay $\propto{\rm
e}^{-c|{\rm Im}\,T(\varepsilon)r|}$ dominates in the $r$-integral and converts
the above criterion into:

$B_{2}(\varepsilon)\approx\frac{8}{81\nu^{3}}\left|\frac{T^{3}(\varepsilon)[\varepsilon-
cT(\varepsilon)]^{3}}{{\rm Im}\,T(\varepsilon)}\right|\ll 1.$ (53)

The straightforward numerical check shows this criterion to surely hold for
all the above considered impurity parameters (see an example in Fig. 13).

Also, the self-consistent extension of the T-matrix function:

$T_{sc}(\varepsilon)=\frac{\omega^{2}}{2}\left\\{\varepsilon-\varepsilon_{res}-\frac{i\omega^{2}}{4\nu\sqrt{1-\left[\frac{\varepsilon-
cT_{sc}(\varepsilon)}{2}\right]^{2}}}\right\\},$ (54)

practically coincides with its non-renormalized version $T(\varepsilon)$ by
Eq. 25, due to the above noted smallness of $|\varepsilon-cT(\varepsilon)|\ll
1$, assured for all $|\varepsilon|\ll 1$.

Hence, the discussed T-matrix results for the quasiparticle spectra in
disordered AGNRs and the related estimates for their mobility edges can be
considered valid. We conclude that the inverse lifetime by the non-
renormalized T-matrix in r.h.s. of Eq. 30 is the main factor for quasiparticle
localization in AGNRs. So the results of Sec. III correctly determine the
system observable characteristics.