Nonstandard Transition GUE‐GOE for Random Matrices and Spectral Statistics of Graphene Nanoflakes

Abstract

Spectral statistics of weakly disordered triangular graphene flakes with zigzag edges are revisited. Earlier, we have found numerically that such systems may show spectral fluctuations of Gaussian unitary ensemble (GUE), signaling the time‐reversal symmetry (TRS) breaking at zero magnetic field, accompanied by approximate twofold valley degeneracy of each energy level. Atomic‐scale disorder induces the scattering of charge carriers between the valleys and restores the spectral fluctuations of Gaussian orthogonal ensemble (GOE). A simplified description of such a nonstandard GUE‒GOE transition, employing the mixed ensemble of 4 × 4 real symmetric matrices was also proposed. Here, we complement our previous study by analyzing numerically the spectral fluctuations of large matrices belonging the same mixed ensemble. Resulting scaling laws relate the ensemble parameter to physical size and the number of atomic‐scale defects in graphene flake. A phase diagram, indicating the regions in which the signatures of GUE may by observable in the size‐doping parameter plane, is presented.

1. Introduction

The notion of emergent phenomena was coined out by Anderson in his milestone science paper of 1979 [1]. In brief, emergence occurs when a complex system shows qualitatively different properties then its building blocks. Numerous examples of emergent systems studied in contemporary condensed matter physics, including high‐temperature superconductors and heavy‐fermion compounds [2], are regarded as systems with spontaneous symmetry breaking [3]. A link between emergence and spontaneous symmetry breaking, however, does not seem to have a permanent character. In a wide class of electronic systems, such as semiconducting heterostructures containing a two‐dimensional electron gas (2DEG), physical properties of itinerant electrons are substantially different than properties of free electrons (or electrons in atoms composing the system), and are also highly‐tunable upon variation of external electromagnetic fields [4]. To give some illustration of this tunability, we only mention that electrons in GaAs heterostructures can be usually described by a standard Schrödinger equation of quantum mechanics with the effective mass $m_{\text{eff}}=0.067m_{\text{e}}$ (where $m_{\text{e}}$ is the free electron mass), whereas in extreme cases of quantum states formed in quantum Hall systems, effective quasiparticles may not even show the Fermi‒Dirac statistics [5, 6].

It is rather rarely noticed that graphene, a two‐dimensional form of carbon just one atom tick [7], also belongs to the second class of emergent systems (i.e., without an apparent spontaneous symmetry breaking) described briefly above. In a monolayer graphene, effective Hamiltonian for low‐energy excitations has a Dirac‒Weyl form, namely

${\mathscr{H}}_{\text{eff}}=v_{\text{F}}\left[p+eA(r,t)\right]\cdot \sigma +U\left(r,t\right),$E1

where $v_{\text{F}}={10}^6$ m/s is the energy‐independent Fermi velocity, $\sigma =({\sigma }_x,{\sigma }_y)$ with the Pauli matrices ${\sigma }_x$ and ${\sigma }_y$, $p=-iħ({\partial }_x,{\partial }_y)$ is the in‐plane momentum operator, the electron charge is $-e$, and the external electromagnetic field is defined via scalar and vector potentials, $U(r, t)$ and $A(r, t)$, with the in‐plane position $r=(x, y)$ and the time $t.$¹

Strictly speaking, ${\mathscr{H}}_{\text{eff}}$ Eq. (1) applies to quasiparticles near the $K$ valley in the dispersion relation. To obtain the effective Hamiltonian for other valley (${\mathrm{K\text{'}}}^{}$) it is sufficient to substitute ${\sigma }_y\to -{\sigma }_y$.

In other words, the system build of nonrelativistic elements (carbon atoms at normal conditions) turns out to host ultrarelativistic quasiparticles, providing a beautiful example of an emergent phenomenon, which binds together two rather distant areas of relativistic quantum mechanics and condensed matter physics [8]. This observation applies generically to bilayer or multilayer graphenes [9], as well as to HgTe/CdTe quantum wells [10], although microscopic models describing such other Dirac systems are slightly different. It is also worth to mention so‐called artificial graphenes, in which waves (of different kinds) obey their effective Dirac equations [11‒13].

A peculiar nature of Dirac fermions in graphene originates from the chiral structure of the Hamiltonian ${\mathscr{H}}_{\text{eff}}$, accompanied by the fact that coupling to the external electromagnetic field is described by additive terms, which are linear in both scalar and vector potentials. A remarkable consequence of these facts is the quantization of the visible light absorption [14], an unexpected macroscopic quantum effect recently found to have analogs in other Dirac systems [15, 16], and even in a familiar graphite [17]. Another intriguing effect of this kind appears for dc conductivity of ballistic graphene [18]. In the so‐called pseudodiffusive transport regime, the conductance of a rectangular sample (with the width $W$ and the length $L$) scales as $G={\sigma }_0\times W/L$ for $W\gg L$, where ${\sigma }_0=(4/\pi )e^2/h$ is the universal quantum value of the conductivity [19, 20], whereas the shot‐noise power and all the other charge‐transfer characteristics are indistinguishable from those of a classical diffusive conductor regardless the sample shape [21]. At high magnetic fields, the pseudo diffusive charge transport is predicted theoretically to reappear for resonances with Landau levels in both monolayer [22] and bilayer graphene [23]. In the presence of disorder, a fundamental property of the Hamiltonian—the time reversal symmetry (TRS)—starts to play a decisive role. In particular, effective TRS in a single valley may be broken even in the absence of magnetic fields, leading to observable (and having the universal character) consequences for the conductance and spectral fluctuations [24, 25], as well as for the peculiar scaling behavior predicted for the conductivity [26, 27].

Although the interest in graphene and other Dirac systems primarily focus on their potential applications [28, 29], quite often linked to the nonstandard quantum description [8], we believe that the fundamental perspective sketched in the above also deserves some attention. In the remaining part of this article, we first overview basic experimental, theoretical and numerical findings concerning signatures of quantum chaos in graphene and its nanostructures (Section 2). Next, we present our new numerical results concerning the additive random matrix model originally proposed in Ref. [25] to describe a nonstandard GUE‒GOE transition, accompanied by lifting out the valley degeneracy (Section 3). The consequences of these findings for prospective experiments on graphene nanoflakes, together with the phase diagram depicting the relevant matrix ensembles in the system size‐doping plane, are described in Section 4. The concluding remarks are given in Section 5.

2. Gauge fields, fluctuations and chaos in nanoscale graphene structures

Dirac fermions confined in graphene quantum dots [30] have provided yet another surprising situation, in which a piece of handbook knowledge needed a careful revision [31].

Quantum chaotic behavior appears generically for systems, whose classical dynamics are chaotic, and manifest itself via the fact that energy levels show statistical fluctuations following those of Gaussian ensembles of random matrices [32]. In particular, if such a system posses the time‐reversal symmetry (TRS), its spectral statistics follows the Gaussian orthogonal ensemble (GOE). A system with TRS and half‐integer spin has the symplectic symmetry and, in turn, shows spectral fluctuations of the Gaussian symplectic ensemble (GSE). If TRS is broken, as in the presence of nontrivial gauge fields, and the system has no other antiunitary symmetry [33], spectral statistics follow the Gaussian unitary ensemble (GUE). For a particular case of massless spin‐1/2 particles, it was pointed out by Berry and Mondragon [34], that the confinement may break TRS in a persistent manner (i.e., even in the absence of gauge fields), leading to the spectral fluctuations of GUE.

When applying the above symmetry classification to graphene nanosystems [24, 25], one needs, however, to take into account that Dirac fermions in graphene appear in the two valleys, $K$ and $K^{\text{'}}$, coupled by TRS. (In particular, real magnetic field breaks TRS and has the same sign in the two valleys, whereas the strain‐induced gauge field preserves TRS and has opposite signs in the two valleys.) If the valley pseudospin is conserved, a special (symplectic) time‐reversal symmetry (STRS) becomes relevant, playing a role of an effective TRS in a single valley [24]. Both real magnetic and strain‐induced gauge fields may break STRS leading to the spectral fluctuations of GUE [35]. As demonstrated numerically in Ref. [25], such fluctuations also appear for particular closed nanosystems in graphene in the presence of random scalar potentials slowly varying on the scale of atomic separation. Such nanosystems include equilateral triangles with zigzag or Klein edges, i.e., with terminal atoms belonging to one sublattice. Generic graphene nanoflakes with irregular edges show spectral fluctuations of GOE [24], as strong intervalley scattering restores TRS in the absence of gauge fields (see Figure 1). In contrast, the boundary effects are suppressed in open graphene systems, for which signatures of the symplectic symmetry class were reported [36].

Figure 1.

Transitions between symmetry classes and random matrix ensembles relevant for closed nanosystems in graphene characterized by the disorder strength, the intervalley scattering rate, and (optionally) placed in the weak magnetic field $B$. (Reprinted with permission from Ref. [25].)

It is worth mention here, that triangular graphene flakes, similar to studied theoretically in Ref. [25], have been recently fabricated [37, 38]. However, due to the hybridization with metallic substrates, quantum‐dot energy levels in such systems are significantly broaden, making it rather difficult to determine the symmetry class via spectral statistics.

3. Transition GUE‒GOE for real symmetric matrices

3.1. Additive random‐matrix models: brief overview

Additive random‐matrix models are capable of reproducing the evolutions of spectral statistics in many cases when a complex system undergoes transition to quantum chaos or transition between symmetry classes [32, 39]. The discussion usually focus on the auxiliary random Hamiltonian of the form

$H\left(\lambda \right)=\frac{H_0+\lambda V}{\sqrt{1+{\lambda }^2}},$E2

where $H_0={(H_0)}^{†}$ and $V=V^{†}$ are members of different Gaussian ensembles,²

To describe transition to quantum chaos rather then transition between symmetry classes in a chaotic system, one can choose $H_0$ to be a diagonal random matrix, elements of which follow a Gaussian distribution with zero mean and the variance ${⟨(H_0)}_{ij}^2⟩={\delta }_{i\mathrm{j.}}$

and the parameter $\lambda \in [0,\infty ]$.

For instance, if elements of $H_0$ are real numbers chosen to follow a Gaussian distribution with zero mean and the variance ${⟨(H_0)}_{ij}^2⟩=(1+{\delta }_{ij})/N$, where ${\delta }_{ij}$ is the Kronecker delta and $N$ is the matrix size, while elements of $V$ are complex numbers in which real and imaginary parts are generated independently according to Gaussian distribution with zero mean and the variance ${⟨(\text{Re}{V}_{ij})}^2⟩=(1+{\delta }_{ij})/2N,$${⟨(\text{Im}{V}_{ij})}^2⟩=(1-{\delta }_{ij})/2N$ (respectively), the Hamiltonian $H(\lambda )$ (2) refers to transition GOE‒GUE. For $N=2$, statistical distribution of the spacing between energy levels $S=|E_1-E_2|$ can be found analytically [40], and reads

$P_{\text{GOE}-\text{GUE}}\left(\lambda ;S\right)=\sqrt{\frac{2+{\lambda }^2}{2}}S{c}^2\left(\lambda \right)\exp \left[-\frac{S^2{c}^2\left(\lambda \right)}{2}\right]\text{erf}\left[\frac{Sc\left(\lambda \right)}{\lambda }\right],$E3

where $\text{erf}(x)$ is the error function, i.e., $\text{erf}(x)=(2/\sqrt{\pi })\underset{0}{\overset{x}{{\displaystyle \int }}}\exp (-t^2)dt$, and

$c\left(\lambda \right)=\sqrt{\frac{\pi \left(2\text{}+\text{}{\lambda }^2\right)}{4}}\left[1-\frac{2}{\pi }\left(\arctan \left(\frac{\lambda }{\sqrt{2}}\right)-\frac{\sqrt{2}\lambda }{2\text{}+\text{}{\lambda }^2}\right)\right].$E4

The above follows from the normalization condition

$⟨S⟩=\underset{0}{\overset{\infty }{{\displaystyle \int }}}S{P}_{\text{GOE}-\text{GUE}}\left(\lambda ;S\right)dS=1 \text{for} 0<\lambda <\infty .$E5

The limiting forms of the spacing distribution given by Eqs. (3) and (4) are

$P_{\text{GOE}-\text{GUE}}\left(\lambda \to 0;S\right)=\frac{\pi }{2}S\exp \left(-\frac{\pi S^2}{4}\right)\equiv P_{\text{GOE}}\left(S\right),$E6

$P_{\text{GOE}-\text{GUE}}\left(\lambda \to \infty ;S\right)=\frac{32}{{\pi }^2}{S}^2\exp \left(-\frac{4S^2}{\pi }\right)\equiv P_{\text{GUE}}\left(S\right),$E7

coinciding with well‐known Wigner surmises for GOE and GUE, respectively [32]. For $N\gg 1$, it was also shown that actual spacing distributions obtained numerically can be approximated (with an astonishing accuracy) by $P_{\text{GOE}-\text{GUE}}({\lambda }_{\text{fit}};S)$, where the empirical parameter ${\lambda }_{\text{fit}}\propto \lambda \sqrt{N}$ [39]. Similar scaling laws apply generically to all transitions between basic symmetry classes.

Relatively recently, spectra of models employing self‐dual random matrices have attracted some attention [41]. In such models, the matrix $H_0$ in Eq. (2) is equivalent (up to a unitary transformation) to the matrix having a block structure

${\tilde{H}}_0=\left(\begin{array}{cc}C& 0\\ 0& C^{\text{T}}\end{array}\right),$E8

where random matrix $C$ is an $N\times N$ member of one of Gaussian ensembles, $C^{\text{T}}$ denotes the transposition of $C$. The matrix $V$ in Eq. (2) is a generic $2N\times 2N$ member of the other ensemble (hereinafter, we redefine the $H(a)$ size as $2N$). In turn, for $\lambda =0$, each eigenvalue is doubly degenerate. For $\lambda \ne 0$, we have the degeneracy splitting accompanied by transition between selected symmetry classes. Even in the simplest case of $N=2$, closed‐form analyticexpressions for level‐spacing distributions corresponding to arbitrary $0<\lambda <\infty$ are missing. The approach presented in Ref. [41] employs the relevant expressions for joint probability densities for eigenvalues [42], allowing one to express level‐spacings distribution in terms of two‐dimensional integrals to be evaluated numerically.

In the remaining part of section, we focus on the transition between self‐dual GUE to GOE, show that the corresponding Hamiltonian $H(\lambda )$, and can be represented as real‐symmetric random matrix, and present our empirical expressions approximating spacing distributions obtained numerically.

3.2. Self‐dual GUE to GOE via $4\text{}\times \text{}4$ real‐symmetric matrices

We focus now on the situation, when the matrix $C$ in Eq. (8) is chosen to be an $N\times N$ member of GUE, whereas $V$ in Eq. (2) is a $2N\times 2N$ member of GOE.

For $N=2$, the matrix ${\tilde{H}}_0$ can be written as

${\tilde{H}}_0^{4\times 4}=\left(\begin{array}{cccc}a& c+id& 0& 0\\ c-id& b& 0& 0\\ 0& 0& a& c-id\\ 0& 0& c+id& b\end{array}\right),$E9

where $a$ and $b$ are real random numbers following Gaussian distribution with zero mean and the variance $⟨a^2⟩=⟨b^2⟩=1/2$, whereas $c$ and $d$ are real random numbers following Gaussian distribution with zero mean and the variance $⟨c^2⟩=⟨d^2⟩=1/4$. Exchanging the second row with the third row, as well as the second column with the third column, we find the matrix ${\tilde{H}}_0^{4\times 4}$ is equivalent, up to an orthogonal transformation, to

${\tilde{H}}_0^{4\times 4} \stackrel{O}{\leftrightarrow } \left(\begin{array}{cccc}a& 0& c+id& 0\\ 0& a& 0& c-id\\ c-id& 0& b& 0\\ 0& c+id& 0& b\end{array}\right).$E10

The matrix on the right‐hand side of Eq. (10) is self‐dual, and can be further transformed as

$U\left(\begin{array}{cccc}a& 0& c+id& 0\\ 0& a& 0& c-id\\ c-id& 0& b& 0\\ 0& c+id& 0& b\end{array}\right)U^{†}=\left(\begin{array}{cccc}a& 0& c& d\\ 0& a& -d& c\\ c& -d& b& 0\\ d& c& 0& b\end{array}\right),$E11

where

$U=\frac{1}{\sqrt{2}}{{I}}_{2\times 2}\otimes \left(\begin{array}{cc}1& 1\\ i& -i\end{array}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}{cccc}1& 1& 0& 0\\ i& -i& 0& 0\\ 0& 0& 1& 1\\ 0& 0& i& -i\end{array}\right).$E12

Exchanging the second with the third row and column in the rightmost matrix in Eq. (11) we arrive to

$H_0^{4\times 4}=\left(\begin{array}{cccc}a& c& 0& d\\ c& b& -d& 0\\ 0& -d& a& c\\ d& 0& c& b\end{array}\right)=\left(\begin{array}{cc}A& B\\ -B& A\end{array}\right),$E13

where the blocks $A$ and $B$ are real‐symmetric ($A^{\text{T}}=A$) and skew‐symmetric ($B^{\text{T}}=-B$) random matrices.

Spectral statistics of the Hamiltonian $H(\lambda )=\left(H_0^{4\times 4}+\lambda V^{4\times 4}\right)/\sqrt{1+{\lambda }^2}$, with $V^{4\times 4}$ being a $4\times 4$ GOE matrix, were thoroughly studied before [43]. Here, we revisit our findings, before discussing spectra of larger matrices in next subsection.

The nearest‐neighbor spacings distribution can be approximated by

$P\left(\alpha ,\kappa ;S\right)=\frac{\alpha P_{\text{GOE}}\left(\alpha S\right)+\beta \left(\alpha \right)P_{\text{GOE}-\text{GUE}}\left(\kappa ;\beta \left(\alpha \right)S\right)}{2},$E14

with $\beta (\alpha )=\alpha /(2\alpha -1)$, $P_{\text{GOE}}(S)$ given by Eq. (6), $P_{\text{GOE}-\text{GUE}}(\kappa ;S)$ given by Eq. (3), and the parameters $\alpha$ and $\kappa$ which can be approximated by empirical functions

$\alpha \approx {\overline{\alpha }}_{4\times 4}\left(\lambda \right)=1.118\times {\left[\sqrt[3]{1+{(0.60/\lambda )}^3}\right]}^{0.98},$E15

and

$\kappa \approx {\overline{\kappa }}_{4\times 4}\left(\lambda \right)=\sqrt{{\left(\frac{1+{\lambda }^{-2}}{1+{(0.33)}^{-2}}\right)}^{0.29}-1}.$E16

Eqs. (15) and (16) represent simplified versions of the corresponding formulas given in Ref. [43]. A comparison with the numerical will be given later in this section.

3.3. Self‐dual GUE to GOE via $2N\text{}\times \text{}2N$ real‐symmetric matrices

We consider now the case of large random matrices ($N\gg 1$). A generalization of the reasoning presented in previous subsection brought as to the unperturbed Hamiltonian $H_0$ with the block structure as given by the last equality in Eq. (13), but $A=A^{\text{T}}$ and $B=-B^{\text{T}}$ are now $N\times N$ random matrices. The elements of each block are independently generated according to a Gaussian distribution with zero mean and the variance $\text{Var}(A_{ij})=(1+{\delta }_{ij})/2N$ and $\text{Var}(B_{ij})=(1-{\delta }_{ij})/2N$, respectively. In turn, $H_0$ can be unitary mapped onto the matrix ${\tilde{H}}_0$ given by Eq. (8) with $C$ being an $N\times N$ member of GUE. The additive random‐matrix model $H(\lambda )$ is complemented with the perturbation $V$ being a $2N\times 2N$ member of GOE.

Ensembles of large pseudo‐random Hamiltonians $H(\lambda )$ were generated and diagonalized numerically, to check whether the standard scaling law ${\lambda }_{\text{fit}}\simeq {(2N)}^{1/2}\lambda$ [44] applies to spacings distribution of such matrices. Our presentation is limited to the matrix sizes $2N=200$, $400$, and $1000$; the statistical ensemble consists of the total amount of ${10}^6$, ${10}^5$, or ${10}^4$ matrices (respectively), same for each considered value of the parameter $\lambda$. To avoid the boundary effects, we limit our numerical study to about $30\%$ of the energy levels such that $\left|E\right|\le 0.5$. Selected examples are presented in Figure 2.

We find that nearest‐neighbor level spacings of large matrix $H(\lambda )$ follow the empirical distribution having the general form as given by Eq. (14).

$P_{\overline{\alpha },\overline{\kappa }}(\tilde{\lambda },S)=P( {\overline{\alpha }}_{N\gg 1}(\tilde{\lambda }),{\overline{\kappa }}_{N\gg 1}(\tilde{\lambda });S),$E17

with the empirical relations of Eqs. (15) and (16) [see blue solid lines in Figure 3] now replaced by

${\overline{\alpha }}_{N\gg 1}\left(\tilde{\lambda }\right)=1.114\times {\left[\sqrt[3]{1+{(0.60/\tilde{\lambda })}^3}\right]}^{0.98},$E18

Figure 2.

Level‐spacing distributions for ${10}^5$ randomly‐generated Hamiltonians $H(\lambda )$ with the size $2N=400$ (datapoints). The scaling parameters λ is varied between the panels. The least‐squares fitted functions $P(\alpha ,\kappa ;S)$ defined by Eq. (14) are also shown (solid lines).

and

${\overline{\kappa }}_{N\gg 1}(\tilde{\lambda })=\{\begin{array}{ll}\sqrt{{\left[\frac{1+{\tilde{\lambda }}^{-2}}{1+{(0.27)}^{-2}}\right]}^{0.29}-1}\hfill & \text{if} \tilde{\lambda }<0.27,\hfill \\ 0\hfill & \text{if} \tilde{\lambda }\ge \mathrm{0.27.}\hfill \end{array}$E19

The above formulae are marked in Figure 3 with red dashed lines. We also find that the scaling law $\tilde{\lambda }={\lambda }_{\text{fit}}\simeq {(2N)}^{1/2}\lambda$ [with $\lambda$ being the original parameter of $H(\lambda )$] is satisfied for the matrices considered with a surprising accuracy (see Figure 4).

Figure 3.

Least‐squares fitted parameters of $P(\alpha ,\kappa ;S)$ Eq. (14) for different values of $2N$ as functions of the scaled model parameter ${(2N)}^{1/2}\lambda$ (datapoints). The empirical relations ${\overline{\alpha }}_{4\times 4}(\lambda )$ Eq. (15) and ${\overline{\kappa }}_{4\times 4}(\lambda )$ Eq. (16) valid for $2N=4$ are shown with blue solid lines; the relations ${\overline{\alpha }}_{N\gg 1}(\lambda )$ Eq. (18) and ${\overline{\kappa }}_{N\gg 1}(\lambda )$ Eq. (19) for large matrices are shown with red dashed lines.

Figure 4.

Scaling law for the best fitted parameters $\tilde{\lambda }={\lambda }_{\text{fit}}$ in the distribution $P_{\overline{\alpha },\overline{\kappa }}(\tilde{\lambda },S)$ Eq. (17) approximating $P(S)$ obtained numerically for random Hamiltonians $H(\lambda )$ with $2N=200$, $400$, and $1000$. [See the main text for details.] Blue solid line marks ${\lambda }_{\text{fit}}=\lambda$.

4. Consequences for graphene nanoflakes

4.1. Level‐spacing distributions revisited

In this section, the empirical distribution $P_{\overline{\alpha },\overline{\kappa }}(\lambda ,S)$ (17) with least‐square fitted $\lambda ={\lambda }_{\text{fit}}$ is utilized to rationalize level‐spacing distributions for triangular graphene nanoflakes with zigzag edges.

At zero magnetic field, the tight‐binding Hamiltonian for weakly‐disordered graphene can be written as

${\mathscr{H}}_{\text{TBA}}={\displaystyle \sum }_{ij}\left[t_{⟨ij⟩}\left|i⟩⟨j\right|+\text{h}.\text{c}.\right]+ {\displaystyle \sum }_i\left[M_V(r_i)+U_{\text{imp}}(r_i)\right]\left|i⟩⟨j\right|,$E20

where $t_{ij}=-t$ if the orbitals $|i$ and $|j$ are nearest neighbors on the honeycomb lattice (with $t=\frac{2}{3}\sqrt{3}ħv_{\text{F}}/a\approx 3$ eV, and $a=0.246$ nm being the lattice spacing), otherwise $t_{ij}=0$. (The symbol

denotes that each pair $⟨ij⟩$ is counted only once.) The terms $M_{\text{V}}(r_i)$ and $U_{\text{imp}}(r_i)$ represent the potentials abruptly and slowly varying on the scale of atomic separation (respectively). Here, we put $M_V(r_i)=0.7t$ if $r_i$ is the outermost atom position at zigzag edge, otherwise $M_{\text{V}}(r_i)=0$. The random contribution $U_{\text{imp}}(r_i)$ is generated in as follows: first, we randomly choose $N_{\text{imp}}$ lattice sites $R_n$ ($n=1,\dots ,N_{\text{imp}}$) out of $N_{\text{tot}}$. Next, the amplitudes $U_n\in (-\delta ,\delta )$ are randomly generated. Finally, the potential is smoothed over a distance $\xi =\sqrt{3}a$ by convolution with a Gaussian, namely

$U_{\text{imp}}\left(r\right)={\displaystyle \sum }_{n=1}^{N_{\text{imp}}}{U}_n\exp \left(-\frac{|r-R_n{|}^2}{2{\xi }^2}\right).$E21

A model of substrate‐induced disorder, constituted by Eqs. (20) and (21), was widely used to reproduce numerically several transport properties of disordered graphene samples [45‒48]. Here, we revisit the spectra of closed graphene flakes considered in Ref. [25], within a simplified empirical model $P_{\overline{\alpha },\overline{\kappa }}({\lambda }_{\text{fit}},S)$ (17), in order discuss the consequences for prospective experimental observation of the zero‐field time‐reversal symmetry breaking in such systems.

A compact measure of the disorder strength is given by the dimensionless correlator

$K_0=\frac{\mathcal{A}}{{\left(ħv_F\right)}^2}\frac{1}{N_{tot}^2}{\displaystyle \sum }_{i=1}^{N_{\text{tot}}}{\displaystyle \sum }_{j=1}^{N_{\text{tot}}}⟨U_{\text{imp}}(r_i)U_{\text{imp}}(r_j)⟩,$E22

where the system area

, and the averaging takes place over possible realizations of the disorder in Eq. (21). For $\xi \gg a$, Eq. (22) leads to

$K_0=\frac{64{\pi }^2\sqrt{3}}{27} \frac{N_{imp}}{N_{tot}}{\left(\frac{\delta }{t}\right)}^2{\left(\frac{\xi }{a}\right)}^4.$E23

For $\xi =\sqrt{3}a$, used for numerical demonstration in the remaining of this article, Eq. (23) still provides a good approximation of the actual value of $K_0$ and can be rewritten as

$K_0\approx 364.7\times \frac{N_{imp}}{N_{tot}}{\left(\frac{\delta }{t}\right)}^2.$E24

Figure 5.

Left: Level‐spacing distributions $P(S)$ for triangular graphene nanoflakes with zigzag edges. The flake area is , the disorder strength is $K_0\approx 0.125$, the number of edge vacancies $N_{\text{vac}}$ is varied betweenthe panels. Numerical results (replotted with permission from Ref. [25]) are shown with black solid lines. The other lines correspond to empirical distributions $P_{\overline{\alpha },\overline{\kappa }}(\lambda ,S)$ Eq. (17) with $\lambda ={\lambda }_{\text{fit}}$ (red solid line), $\lambda =0$ (blue dashed line) or $\lambda =\infty$ (blue dotted line). Right: Least‐squares fitted parameters for different numbers of edge vacancies $1\le N_{\text{vac}}\le 30$ and the flake areas (open symbols) and (closed symbols), corresponding to the total number of terminal atoms $N_{\text{edge}}=270$ and $540$ (respectively). Solid line depicts the approximating power‐law relation given by Eq. (25).

The numerical results are presented in Figure 5, where we have fixed the remaining disorder parameters at $\delta /t=0.1$ and $N_{\text{imp}}/N_{\text{tot}}=0.034$ leading to $K_0=0.125$.³

The disorder parameters are actually same as in Figures 8 and 9 of Ref. [25], where we have mistakenly omitted the factor $\pi$ in the numerical evaluation of $K_0$.

Level‐spacing distributions $P(S)$ obtained numerically for triangular nanoflakes with zigzag edges [see left panels in Figure 5, black solid lines] are replotted with permission from Ref. [25], where we used approximately $1500$ energy levels with energies $0.1\le \left|E\right|/t\le 0.5$ out of the total number of $N_{\text{tot}}(N_{\text{vac}})=32758-N_{\text{vac}}$ (corresponding the flake area

, with $N_{\text{vac}}$ being the number of vacancies, randomly distributed along the system boundary. Typically, best‐fitted parameters $\lambda ={\lambda }_{\text{fit}}$ of the simplified distribution $P_{\overline{\alpha },erline\kappa }(\lambda ,S)$ (17) coincide with given in Ref. [25] up to a second decimal place. New values of ${\lambda }_{\text{fit}}$ for $1\le N_{\text{vac}}\le 30$ and two flake sizes $N_{\text{tot}}(0)=8278$ and $N_{\text{tot}}(0)=32,758$ are displayed in the right panel of Figure 5. The dependence of ${\lambda }_{\text{fit}}$ on $N_{\text{vac}}$ and $N_{\text{tot}}$ can be rationalize within a power‐law

${\lambda }_{\text{fit}}\approx 0.103\times {\left(N_{\text{vac}}\sqrt{N_{\text{edge}}}\right)}^{0.34},$E25

where the total number of terminal sites

$N_{\text{edge}}=3\sqrt{N_{\text{tot}}+N_{\text{vac}}+3}-3.$E26

4.2. Phase diagram for triangular flakes with zigzag edges

Eq. (25) is now employed to estimate the maximal system size $N_{\text{tot}}$, and the maximal number of edge vacancies $N_{\text{vac}}$, for which signatures of TRS breaking still can be identified in the spectrum. This is possible as long as ${\lambda }_{\text{fit}}<{\lambda }_{\star }=0.27$ (see Eq. (19)), as for any ${\lambda }_{\text{fit}}\ge {\lambda }_{\star }$ we have $\overline{\kappa }({\lambda }_{\text{fit}})=0$ and level‐spacing distribution simply evolves from that characterizing GOE matrix with approximate twofold degeneracy of each level toward GOE without such a degeneracy. For instance, we obtain

$N_{\text{tot}}\lesssim 9500 \text{for} N_{\text{vac}}=1,$E27

$N_{\text{tot}}\lesssim 630 \text{for} N_{\text{vac}}=2,$E28

$N_{\text{tot}}\lesssim 130 \text{for} N_{\text{vac}}=3.$E29

On the other hand, system size and the number of energy levels taken into account must be large enough to distinguish between spectral fluctuations of GUE and spectral fluctuations of other ensembles.

Density of states (per one direction of spin) for bulk graphene reads

${\rho }_{\text{bulk}}(E)=\frac{\mathcal{A}}{\pi {(ħv_F)}^2}\left|E\right|=\frac{N_{tot}}{\sqrt{3}\pi t^2}\left|E\right|.$E30

The number of energy levels $N_{\text{lev}}$ in the interval $(0,E_{\text{max}})$ can be approximated by

$N_{\text{lev}}\approx \underset{0}{\overset{E_{\text{max}}}{{\displaystyle \int }}}{\rho }_{\text{bulk}}\left(E\right)dE\approx 0.0919\times N_{\text{tot}}{\left(\frac{E_{max}}{t}\right)}^2.$E31

Physically, occupying $N_{\text{lev}}$ electronic levels above the Dirac point one produces the electric charge $Q=-2e{N}_{\text{lev}}$, resulting in a typical experimental limit of $E_{\text{max}}=0.2-0.3$ eV for graphene nanostructures on SiO₂‐based substrates [49].

Level‐spacing distributions $P(S)$ are normalized such that

. In turn, the variance

raises as the lowest moment allowing one to distinguish between different distributions. In particular, we have

$\text{Var}{\left\{S\right\}}_{\text{GOE}}=\frac{4}{\pi }-1\approx 0.273, \text{Var}{\left\{S\right\}}_{\text{GUE}}=\frac{3\pi }{8}-1\approx 0.178,$E32

$\text{and} \text{Var}{\left\{S\right\}}_{\lambda \to 0}=\frac{3\pi }{4}-1\approx 1.356,$E33

where Eq. (33) refers to the empirical distribution $P_{\overline{\alpha },\overline{\kappa }}(\lambda ,S)$ given by Eq. (17) with $\lambda \to 0$. Similar calculation for arbitrary $\lambda$ is straightforward, but the resulting formula is too lengthy to be presented.⁴

We use the property of $m$-th cumulant of the distribution $P(S)=\frac{1}{2}[a{P}_1(aS)+b{P}_2(bS)]$, which is equal to ${⟨S^m⟩}_P=1/2(a^{(-m)}{⟨S^m⟩}_{(P_1)}+b^{(-m)}{⟨S^m⟩}_{(P_2)}).$. For $P_1=P_{\text{GOE}}$ and $P_2=P_{\text{GOE}-\text{GUE}}$, see Eqs. (6) and (3), necessary integrals for $m=2, 3, 4$ can be calculated analytically.

When $\text{Var}\left\{S\right\}$ is calculates for a large but finite collection of spacings $N_{\text{spc}}=N_{\text{lev}}\text{}-\text{}1$, it becomes a random variable itself, with a variance which can be approximated by

$\text{Var}\left\{\text{Var}\left\{S\right\}\right\}\approx \frac{1}{N_{\text{lev}}}\left({\mu }_4-{\sigma }^2\right)=\frac{1}{N_{\text{lev}}}\left(⟨S^4⟩-4⟨S^3⟩+8⟨S^2⟩-{⟨S^2⟩}^2-4\right),$E34

where ${\mu }_4={⟨\left(⟨S⟩-S\right)⟩}^4$ denotes the fourth central moment and we have used the normalization $⟨S⟩=1$. In turn, for spacings following the distribution $P_{\overline{\alpha },\overline{\kappa }}(\lambda ,S)$ (17) one can find they do not follow GOE if

$N_{\text{lev}}\gtrsim {\left[\frac{3\sqrt{{\left({\mu }_4-{\sigma }^2\right)}_{\lambda }}}{\text{Var}{\left\{S\right\}}_{\lambda }-\text{Var}{\left\{S\right\}}_{\text{GOE}}}\right]}^2,$E35

where the factor $3$ in the nominator corresponds to the $3\sigma$ level of significance. Substituting Eq. (31) one can rewrite the above as

$E_{\text{max}}\sqrt{N_{\text{tot}}} \gtrsim 9.9t\times \frac{\sqrt{{\left({\mu }_4-{\sigma }^2\right)}_{\lambda }}}{\text{Var}{\left\{S\right\}}_{\lambda }-0.273}.$E36

For $\lambda \to 0$, we have ${\mu }_4-{\sigma }^2\to \frac{21}{16}{\pi }^2-2\pi -4\approx 2.671$, leading to

$N_{\text{tot}} \gtrsim 223\times {\left(\frac{t}{E_{\text{max}}}\right)}^2\text{for}{N}_{\text{vac}}=0.$E37

Figure 6.

Phase diagram for triangular graphene nanoflakes with zigzag edges. Grey solid line in left panel (replotted as a dashed line in right panel) corresponding to Eq. (37) for $N_{\text{vac}}=0$ splits the region where number of available energy levels is insufficient to determine the class of spectral fluctuation (below the line) and the region where one should be able to identify the unitary class with approximate twofold degeneracy (2 × GUE). Blue solid line in right panel is same as solid line in left panel, but for $N_{\text{vac}}=1$, calculated numerically from Eq. (13) for $\lambda ={\lambda }_{\text{fit}}(N_{\text{tot}})$ (see Eq. (25)). Vertical red line in right panel marks the limit given by Eq. (27), above which the orthogonal class with gradual degeneracy splitting (2 × GOE $\to$ GOE) appears.

Limiting values of $N_{\text{tot}}$ and $E_{\text{max}}$, following from Eqs. (27), (36), and (37) are in depicted Figure 6, presenting the central results of this work. In the absence of edge vacancies ($N_{\text{vac}}=0$), the attainable Fermi energy $E_{\text{max}}=0.25$ eV should make it possible to detect TRS breaking in nanoflakes containing $N_{\text{tot}}\gtrsim 3{10}^4$ carbon atoms, corresponding to the physical diameter of

. For $N_{\text{vac}}=1$, the limit of $N_{\text{tot}}\lesssim 9500$ (see Eq. (27)) implies $E_{\text{max}}\gtrsim 0.8$ eV is required, slightly exceeding current experimental limits for graphene nanostructures.

5. Concluding remarks

We have revisited level‐spacing statistics of triangular graphene nanoflakes with zigzag edges, subjected to weak substrate‐induced disorder. Our previous study of the system is complemented by comparing the spectral fluctuations with these of large random matrices belonging to a mixed ensemble interpolating between GUE with self‐dual symmetry and generic GOE. The results show that for a fixed value of maximal Fermi energy $E_{\text{max}}$ (in typical experiment, the Fermi energy is tuned in the range $-E_{\text{max}}<E<E_{\text{max}}$ by top gate electrode), the system size required to detect signatures of the time‐reversal symmetry breaking at zero magnetic field is bounded from the bottom by the condition for minimal number of quantum‐dot energy levels allowing one to distinguish between different classes of spectral fluctuations. A finite number of vacancies at the system boundary may lead to intervalley scattering restoring TRS, resulting in additional, upper limit for the system size.

In conclusion, we expect that triangular graphene flakes with perfect zigzag edges may show signatures of TRS breaking starting from physical sizes exceeding $15$ nm. For a finite number of atomic‐scale defects (starting from a single edge vacancy), one should search for signatures of the unitary symmetry class in artificial graphene‐like systems rather then in real graphene nanoflakes.

Acknowledgments

The author thanks to Huang Liang for the correspondence. The work was supported by the National Science Centre of Poland (NCN) via Grant no. 2014/14/E/ST3/00256. Computations were partly performed using the PL‐Grid infrastructure.