**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

Nonstandard Transition GUE‐GOE for Random Matrices and Spectral Statistics of Graphene Nanoflakes 595 http://dx.doi.org/10.5772/64240

$$H(\lambda) = \frac{H\_0 + \lambda V}{\sqrt{1 + \lambda^2}} \tag{2}$$

where *H*<sup>0</sup> =(*H*0)† and *<sup>V</sup>* <sup>=</sup>*<sup>V</sup>* † are members of different Gaussian ensembles,2 and the parameter *λ* ∈ 0, *∞* .

For instance, if elements of *H*<sup>0</sup> are real numbers chosen to follow a Gaussian distribution with zero mean and the variance (*H*0) *ij* <sup>2</sup> =(1 + *δij* ) / *N* , where *δ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 (Re*Vij* )<sup>2</sup> =(1 + *δij* ) / 2*N* , (Im*Vij* )<sup>2</sup> =(1−*δij* )/ 2*N* (respectively), the Hamiltonian *H* (*λ*) (2) refers to transition GOE‒GUE. For *N* =2, statistical distribution of the spacing between energy levels *S* = | *E*<sup>1</sup> − *E*<sup>2</sup> | can be found analytically [40], and reads

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

where erf(*x*) is the error function, i.e., erf(*x*)=(2 / *<sup>π</sup>*)*∫* 0 *x* exp (−*t* <sup>2</sup> )*dt*, and

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

The above follows from the normalization condition

[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

**Figure 1.** Transitions between symmetry classes and random matrix ensembles relevant for *closed* nanosystems in gra‐ phene characterized by the disorder strength, the intervalley scattering rate, and (optionally) placed in the weak mag‐

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,

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

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**

Hamiltonian of the form

signatures of the symplectic symmetry class were reported [36].

944 Recent Advances in Graphene Research

netic field *B*. (Reprinted with permission from Ref. [25].)

$$\mathcal{S}\left\langle\mathcal{S}\right\rangle = \int\_0^\infty \mathcal{S} \, P\_{\text{GOE} - \text{GUE}}\left\langle\mathcal{S}; \mathcal{S}\right\rangle dS = 1 \quad \text{for} \quad 0 < \lambda < \infty. \tag{5}$$

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}\mathcal{S}\exp\left(-\frac{\pi S^2}{4}\right) \equiv P\_{\text{GOE}}\left\{\mathcal{S}\right\},\tag{6}$$

$$P\_{\text{GUE}\\_\text{GUE}}\left(\lambda \to \infty; \mathcal{S}\right) = \frac{32}{\pi^2} \mathcal{S}^2 \exp\left(-\frac{4S^2}{\pi}\right) \equiv P\_{\text{GUE}}\left(\mathcal{S}\right),\tag{7}$$

coinciding with well‐known Wigner surmises for GOE and GUE, respectively [32]. For *N* ≫1, it was also shown that actual spacing distributions obtained numerically can be

<sup>2</sup> 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* <sup>2</sup> =*δi*j.

approximated (with an astonishing accuracy) by *P*GOE−GUE(*λ*fit;*S*), where the empirical param‐ eter *λ*fit∝*λ 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 = \begin{pmatrix} C & 0 \\ 0 & C^T \end{pmatrix}. \tag{8}
$$

where random matrix *C* is an *N* × *N* member of one of Gaussian ensembles, *C* T denotes the transposition of *C*. The matrix *V* in Eq. (2) is a generic 2*N* ×2*N* member of the other ensemble (hereinafter, we redefine the *H* (*a*) size as 2*N* ). In turn, for *λ* =0, each eigenvalue is doubly degenerate. For *λ* ≠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<*λ* <*∞* 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‐dimen‐ sional 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* (*λ*), 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×4 real‐symmetric matrices**

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

For *<sup>N</sup>* =2, the matrix *H*˜ <sup>0</sup> can be written as

$$
\hat{H}\_0^{4 \times 4} = \begin{pmatrix} a & c+id & 0 & 0 \\ c-id & b & 0 & 0 \\ 0 & 0 & a & c-id \\ 0 & 0 & c+id & b \end{pmatrix} \tag{9}
$$

where *a* and *b* are real random numbers following Gaussian distribution with zero mean and the variance *a* <sup>2</sup> = *b* <sup>2</sup> =1 / 2, whereas *c* and *d* are real random numbers following Gaussian distribution with zero mean and the variance *c* <sup>2</sup> = *d* <sup>2</sup> =1 / 4. Exchanging the second row with the third row, as well as the second column with the third column, we find the matrix *H*˜ 0 4×4 is equivalent, up to an orthogonal transformation, to

Nonstandard Transition GUE‐GOE for Random Matrices and Spectral Statistics of Graphene Nanoflakes 797 http://dx.doi.org/10.5772/64240

$$
\widetilde{H}\_0^{4 \times 4} \overset{0}{\leftrightarrow} \begin{pmatrix} a & 0 & c+id & 0 \\ 0 & a & 0 & c-id \\ c-id & 0 & b & 0 \\ 0 & c+id & 0 & b \end{pmatrix} . \tag{10}
$$

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

$$U\begin{pmatrix} a & 0 & c+id & 0\\ 0 & a & 0 & c-id\\ c-id & 0 & b & 0\\ 0 & c+id & 0 & b \end{pmatrix} U^\dagger = \begin{pmatrix} a & 0 & c & d\\ 0 & a & -d & c\\ c & -d & b & 0\\ d & c & 0 & b \end{pmatrix} \tag{11}$$

where

(8)

(9)

0 4×4 is

approximated (with an astonishing accuracy) by *P*GOE−GUE(*λ*fit;*S*), where the empirical param‐ eter *λ*fit∝*λ N* [39]. Similar scaling laws apply generically to all transitions between basic

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

where random matrix *C* is an *N* × *N* member of one of Gaussian ensembles, *C* T denotes the transposition of *C*. The matrix *V* in Eq. (2) is a generic 2*N* ×2*N* member of the other ensemble (hereinafter, we redefine the *H* (*a*) size as 2*N* ). In turn, for *λ* =0, each eigenvalue is doubly degenerate. For *λ* ≠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<*λ* <*∞* 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‐dimen‐

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

We focus now on the situation, when the matrix *C* in Eq. (8) is chosen to be an *N* × *N* member

where *a* and *b* are real random numbers following Gaussian distribution with zero mean and the variance *a* <sup>2</sup> = *b* <sup>2</sup> =1 / 2, whereas *c* and *d* are real random numbers following Gaussian distribution with zero mean and the variance *c* <sup>2</sup> = *d* <sup>2</sup> =1 / 4. Exchanging the second row with the third row, as well as the second column with the third column, we find the matrix *H*˜

symmetry classes.

966 Recent Advances in Graphene Research

transformation) to the matrix having a block structure

sional integrals to be evaluated numerically.

**3.2. Self‐dual GUE to GOE via 4×4 real‐symmetric matrices**

of GUE, whereas *V* in Eq. (2) is a 2*N* ×2*N* member of GOE.

equivalent, up to an orthogonal transformation, to

<sup>0</sup> can be written as

obtained numerically.

For *<sup>N</sup>* =2, the matrix *H*˜

$$U = \frac{1}{\sqrt{2}} \mathbb{I}\_{2 \times 2} \otimes \begin{pmatrix} 1 & 1 \\ i & -i \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 & 0 & 0 \\ i & -i & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & i & -i \end{pmatrix} . \tag{12}$$

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

$$H\_0^{4 \times 4} = \begin{pmatrix} a & c & 0 & d \\ c & b & -d & 0 \\ 0 & -d & a & c \\ d & 0 & c & b \end{pmatrix} = \begin{pmatrix} A & B \\ -B & A \end{pmatrix},\tag{13}$$

where the blocks *A* and *B* are real‐symmetric (*A*<sup>T</sup> = *A*) and skew‐symmetric (*B* <sup>T</sup> = − *B*) random matrices.

Spectral statistics of the Hamiltonian *H* (*λ*)=(*H*<sup>0</sup> 4×4 + *λV* 4×4 ) / 1 + *λ* <sup>2</sup> , with *V* 4×4 being a 4×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(\boldsymbol{a},\kappa;\mathcal{S}) = \frac{aP\_{\text{GOL}}\left(a\mathcal{S}\right) + \beta\left(a\right)P\_{\text{GOL}}\\_{\text{GOL}}\left(\kappa;\beta\left(a\right)\mathcal{S}\right)}{2},\tag{14}$$

with *β*(*α*)=*α* / (2*α* −1), *P*GOE(*S*) given by Eq. (6), *P*GOE−GUE(*κ*;*S*) given by Eq. (3), and the param‐ eters *α* and *κ* which can be approximated by empirical functions

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

and

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

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 2***N* **×2***N* **real‐symmetric matrices**

We consider now the case of large random matrices (*N* ≫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*T and *B* = − *B* T are now *N* × *N* random matrices. The elements of each block are independently generated according to a Gaussian distribution with zero mean and the variance Var(*Aij* )=(1 + *δij* ) / 2*N* and Var(*Bij* )=(1−*δij* )/ 2*N* , respectively. In turn, *H*0 can be unitary mapped onto the matrix *H*˜ 0 given by Eq. (8) with *C* being an *N* × *N* member of GUE. The additive random‐matrix model *H* (*λ*) is complemented with the perturbation *V* being a 2*N* ×2*N* member of GOE.

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

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

$$P\_{\overline{a,\mathsf{K}}}(\vec{\lambda},\mathcal{S}) = P(\mathsf{Z}\_{\mathbb{N}\gg 1}(\vec{\lambda}), \overline{\mathsf{x}}\_{\mathbb{N}\gg 1}(\vec{\lambda}); \mathcal{S}),\tag{17}$$

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

$$\overline{a}\_{N\gg1}\{\tilde{\lambda}\} = 1.114 \times \left[ \sqrt[3]{1 + (0.60/\tilde{\lambda})^3} \right]^{0.98},\tag{18}$$

Nonstandard Transition GUE‐GOE for Random Matrices and Spectral Statistics of Graphene Nanoflakes 999 http://dx.doi.org/10.5772/64240

**Figure 2.** Level‐spacing distributions for 10<sup>5</sup> randomly‐generated Hamiltonians *H* (*λ*) with the size 2*N* =400 (data‐ points). The scaling parameters λ is varied between the panels. The least‐squares fitted functions *P*(*α*, *κ*;*S*) defined by Eq. (14) are also shown (solid lines).

and

(16)

Eqs. (15) and (16) represent simplified versions of the corresponding formulas given in Ref.

We consider now the case of large random matrices (*N* ≫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*T and *B* = − *B* T are now *N* × *N* random matrices. The elements of each block are independently generated according to a Gaussian

being an *N* × *N* member of GUE. The additive random‐matrix model *H* (*λ*) is complemented

Ensembles of large pseudo‐random Hamiltonians *H* (*λ*) were generated and diagonalized

distribution of such matrices. Our presentation is limited to the matrix sizes 2*N* =200, 400, and

tively), same for each considered value of the parameter *λ*. To avoid the boundary effects, we limit our numerical study to about 30*%* of the energy levels such that |*E* | ≤0.5. Selected

We find that nearest‐neighbor level spacings of large matrix *H* (*λ*) follow the empirical

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

)=(1 + *δij*

) / 2*N* and Var(*Bij*

1/2

, 10<sup>5</sup>

, or 10<sup>4</sup>

)=(1−*δij*

0 given by Eq. (8) with *C*

*λ* [44] applies to spacings

matrices (respec‐

)/ 2*N* ,

(17)

(18)

[43]. A comparison with the numerical will be given later in this section.

**3.3. Self‐dual GUE to GOE via 2***N* **×2***N* **real‐symmetric matrices**

respectively. In turn, *H*0 can be unitary mapped onto the matrix *H*˜

numerically, to check whether the standard scaling law *λ*fit≃(2*N* )

1000; the statistical ensemble consists of the total amount of 10<sup>6</sup>

distribution having the general form as given by Eq. (14).

examples are presented in **Figure 2**.

by

distribution with zero mean and the variance Var(*Aij*

with the perturbation *V* being a 2*N* ×2*N* member of GOE.

and

988 Recent Advances in Graphene Research

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

The above formulae are marked in **Figure 3** with red dashed lines. We also find that the scaling law *λ*˜ <sup>=</sup>*λ*fit≃(2*<sup>N</sup>* ) 1/2 *λ* [with *λ* being the original parameter of *H* (*λ*)] is satisfied for the matrices considered with a surprising accuracy (see **Figure 4**).

**Figure 3.** Least‐squares fitted parameters of *P*(*α*, *κ*;*S*) Eq. (14) for different values of 2*N* as functions of the scaled model parameter (2*N* ) 1/2 *<sup>λ</sup>* (datapoints). The empirical relations *α*¯ <sup>4</sup>×4(*λ*) Eq. (15) and *κ*¯ <sup>4</sup>×4(*λ*) Eq. (16) valid for <sup>2</sup>*<sup>N</sup>* <sup>=</sup><sup>4</sup> are shown with blue solid lines; the relations *α*¯ *<sup>N</sup>* <sup>≫</sup>1(*λ*) Eq. (18) and *κ*¯ *<sup>N</sup>* <sup>≫</sup>1(*λ*) Eq. (19) for large matrices are shown with red dashed lines.

**Figure 4.** Scaling law for the best fitted parameters *λ*˜ <sup>=</sup>*λ*fit in the distribution *Pα*¯,*κ*¯(*λ*˜, *<sup>S</sup>*) Eq. (17) approximating *P*(*S*) obtained numerically for random Hamiltonians *H* (*λ*) with 2*N* =200, 400, and 1000. [See the main text for details.] Blue solid line marks *λ*fit=*λ*.

#### **4. Consequences for graphene nanoflakes**

#### **4.1. Level‐spacing distributions revisited**

The above formulae are marked in **Figure 3** with red dashed lines. We also find that the scaling

**Figure 3.** Least‐squares fitted parameters of *P*(*α*, *κ*;*S*) Eq. (14) for different values of 2*N* as functions of the scaled

<sup>2</sup>*<sup>N</sup>* <sup>=</sup><sup>4</sup> are shown with blue solid lines; the relations *α*¯ *<sup>N</sup>* <sup>≫</sup>1(*λ*) Eq. (18) and *κ*¯ *<sup>N</sup>* <sup>≫</sup>1(*λ*) Eq. (19) for large matrices are

**Figure 4.** Scaling law for the best fitted parameters *λ*˜ <sup>=</sup>*λ*fit in the distribution *Pα*¯,*κ*¯(*λ*˜, *<sup>S</sup>*) Eq. (17) approximating *P*(*S*) obtained numerically for random Hamiltonians *H* (*λ*) with 2*N* =200, 400, and 1000. [See the main text for

*<sup>λ</sup>* (datapoints). The empirical relations *α*¯ <sup>4</sup>×4(*λ*) Eq. (15) and *κ*¯ <sup>4</sup>×4(*λ*) Eq. (16) valid for

*λ* [with *λ* being the original parameter of *H* (*λ*)] is satisfied for the matrices

law *λ*˜ <sup>=</sup>*λ*fit≃(2*<sup>N</sup>* )

10010 Recent Advances in Graphene Research

model parameter (2*N* )

shown with red dashed lines.

details.] Blue solid line marks *λ*fit=*λ*.

1/2

1/2

considered with a surprising accuracy (see **Figure 4**).

In this section, the empirical distribution *Pα*¯,*κ*¯(*λ*, *S*) (17) with least‐square fitted *λ* =*λ*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

$$\mathcal{H}\_{\text{TBA}} = \Sigma\_{\langle ij \rangle} \left[ t\_{ij} \langle i \rangle \langle j \rangle + \text{h.c.} \right] + \Sigma\_{\langle} \left[ M\_V(\mathbf{r}\_i) + U\_{\text{imp}}(\mathbf{r}\_l) \right] \left| i \rangle \langle i \rangle. \tag{20}$$

where *tij* = −*t* if the orbitals |*i* and | *j* are nearest neighbors on the honeycomb lattice (with *<sup>t</sup>* <sup>=</sup> <sup>2</sup> <sup>3</sup> 3*ħv*<sup>F</sup> / *a* ≈3 eV, and *a* =0.246 nm being the lattice spacing), otherwise *tij* =0. (The symbol denotes that each pair *ij* is counted only once.) The terms *M*V(*ri* ) and *U*imp(*ri* ) represent the potentials abruptly and slowly varying on the scale of atomic separation (respectively). Here, we put *MV* (*ri* )=0.7 *t* if *ri* is the outermost atom position at zigzag edge, otherwise *M*V(*ri* )=0. The random contribution *U*imp(*ri* ) is generated in as follows: first, we randomly choose *N*imp lattice sites *Rn* (*n* =1, …, *N*imp) out of *N*tot. Next, the amplitudes *Un* ∈(−*δ*, *δ*) are randomly generated. Finally, the potential is smoothed over a distance *ξ* = 3 *a* by convolution with a Gaussian, namely

$$U\_{\rm imp} \left( r \right) = \sum\_{n=1}^{N\_{\rm imp}} U\_n \exp \left( -\frac{|r - \mathbb{R}\_n|^2}{2\xi^2} \right). \tag{21}$$

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 simpli‐ fied empirical model *Pα*¯,*κ*¯(*λ*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}}{(\hbar v\_F)^2} \frac{1}{N\_{tot}^2} \sum\_{i=1}^{N\_{tot}} \sum\_{j=1}^{N\_{tot}} \{U\_{\rm imp} \left(\mathbf{r}\_i\right) U\_{\rm imp} \left(\mathbf{r}\_f\right)\},\tag{22}$$

where the system area , and the averaging takes place over possible realizations of the disorder in Eq. (21). For *ξ*≫*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. \tag{23}$$

For *ξ* = 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 3.64.7 \times \frac{N\_{imp}}{N\_{tot}} \left(\frac{\delta}{t}\right)^2. \tag{24}$$

**Figure 5.** *Left:* Level‐spacing distributions *P*(*S*) for triangular graphene nanoflakes with zigzag edges. The flake area is , the disorder strength is *K*<sup>0</sup> ≈0.125, the number of edge vacancies *N*vac is varied betweenthe panels. Numerical results (replotted with permission from Ref. [25]) are shown with black solid lines. The other lines corre‐ spond to empirical distributions *Pα*¯,*κ*¯(*λ*, *S*) Eq. (17) with *λ* =*λ*fit (red solid line), *λ* =0 (blue dashed line) or *λ* =*∞* (blue dotted line). *Right:* Least‐squares fitted parameters for different numbers of edge vacancies 1≤ *N*vac≤30 and the flake areas (open symbols) and (closed symbols), corresponding to the total number of terminal atoms *N*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 disor‐ der parameters at *δ* / *t* =0.1 and *N*imp / *N*tot=0.034 leading to *K*<sup>0</sup> =0.125. 3 Level‐spacing distribu‐ tions *P*(*S*) obtained numerically for triangular nanoflakes with zigzag edges [see left panels

<sup>3</sup> The disorder parameters are actually same as in Figures 8 and 9 of Ref. [25], where we have mistakenly omitted the factor *π* in the numerical evaluation of *K*0.

in **Figure 5**, black solid lines] are replotted with permission from Ref. [25], where we used approximately 1500 energy levels with energies 0.1≤ |*E* | / *t* ≤0.5 out of the total number of *N*tot(*N*vac)=32758− *N*vac (corresponding the flake area , with *N*vac being the number of vacancies, randomly distributed along the system boundary. Typically, best‐fit‐ ted parameters *λ* =*λ*fit of the simplified distribution *Pα*¯,*erlineκ*(*λ*, *S*) (17) coincide with given in Ref. [25] up to a second decimal place. New values of *λ*fit for 1≤ *N*vac≤30 and two flake sizes *N*tot(0)=8278 and *N*tot(0)=32, 758 are displayed in the right panel of **Figure 5**. The depend‐ ence of *λ*fit on *N*vac and *N*tot can be rationalize within a power‐law

$$
\lambda\_{\rm flt} \approx 0.103 \times \left( N\_{\rm vac} \sqrt{N\_{\rm edge}} \right)^{0.34}, \tag{25}
$$

where the total number of terminal sites

For *ξ* = 3 *a*, used for numerical demonstration in the remaining of this article, Eq. (23) still

**Figure 5.** *Left:* Level‐spacing distributions *P*(*S*) for triangular graphene nanoflakes with zigzag edges. The flake area is

Numerical results (replotted with permission from Ref. [25]) are shown with black solid lines. The other lines corre‐ spond to empirical distributions *Pα*¯,*κ*¯(*λ*, *S*) Eq. (17) with *λ* =*λ*fit (red solid line), *λ* =0 (blue dashed line) or *λ* =*∞* (blue dotted line). *Right:* Least‐squares fitted parameters for different numbers of edge vacancies 1≤ *N*vac≤30 and the flake areas (open symbols) and (closed symbols), corresponding to the total number of terminal atoms *N*edge=270 and 540 (respectively). Solid line depicts the approximating power‐law relation given by

The numerical results are presented in **Figure 5**, where we have fixed the remaining disor‐

tions *P*(*S*) obtained numerically for triangular nanoflakes with zigzag edges [see left panels

<sup>3</sup> The disorder parameters are actually same as in Figures 8 and 9 of Ref. [25], where we have mistakenly omitted the

der parameters at *δ* / *t* =0.1 and *N*imp / *N*tot=0.034 leading to *K*<sup>0</sup> =0.125.

Eq. (25).

factor *π* in the numerical evaluation of *K*0.

, the disorder strength is *K*<sup>0</sup> ≈0.125, the number of edge vacancies *N*vac is varied betweenthe panels.

3

Level‐spacing distribu‐

(24)

provides a good approximation of the actual value of *K*0 and can be rewritten as

10212 Recent Advances in Graphene Research

$$N\_{\rm edge} = \Im \sqrt{N\_{\rm tot} + N\_{\rm vac} + 3} - 3. \tag{26}$$

#### **4.2. Phase diagram for triangular flakes with zigzag edges**

Eq. (25) is now employed to estimate the maximal system size *N*tot, and the maximal number of edge vacancies *N*vac, for which signatures of TRS breaking still can be identified in the spectrum. This is possible as long as *λ*fit<*λ*<sup>⋆</sup> =0.27 (see Eq. (19)), as for any *λ*fit≥*λ*⋆ we have *<sup>κ</sup>*¯(*λ*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 \quad \text{for} \quad N\_{\text{vac}} = 1,\tag{27}$$

$$N\_{\text{tot}} \lesssim 630 \quad \text{for} \quad N\_{\text{vac}} = 2,\tag{28}$$

$$N\_{\text{tot}} \lesssim 130 \quad \text{for} \quad N\_{\text{vac}} = 3. \tag{29}$$

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}}\left(E\right) = \frac{\mathcal{A}}{\pi (\hbar v\_{\text{F}})^2} \left| E\right| = \frac{N\_{\text{tot}}}{\sqrt{3} \,\pi t^2} \left| E\right|. \tag{30}$$

The number of energy levels *N*lev in the interval (0, *E*max) can be approximated by

$$N\_{\rm lev} \approx \int\_0^{E\_{\rm max}} \rho\_{\rm bulk} \text{ ( $E$ )} dE \approx 0.0919 \times N\_{\rm tot} \left(\frac{E\_{\rm max}}{\rm t}\right)^2. \tag{31}$$

Physically, occupying *N*lev electronic levels above the Dirac point one produces the electric charge *Q* = −2*eN*lev, resulting in a typical experimental limit of *E*max =0.2−0.3 eV for graphene nanostructures on SiO2‐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}\{\mathcal{S}\}\_{\text{GOE}} = \frac{4}{\pi} - 1 \approx 0.273, \quad \text{Var}\{\mathcal{S}\}\_{\text{GOE}} = \frac{3\pi}{8} - 1 \approx 0.178,\tag{32}$$

$$\text{Var}\{\mathcal{S}\}\_{\lambda \to 0} = \frac{3\pi}{4} - 1 \approx 1.356,\tag{33}$$

where Eq. (33) refers to the empirical distribution *Pα*¯,*κ*¯(*λ*, *S*) given by Eq. (17) with *λ* →0. Similar calculation for arbitrary *λ* is straightforward, but the resulting formula is too lengthy to be presented.<sup>4</sup> When Var{*S*} is calculates for a large but finite collection of spacings *N*spc= *N*lev −1, it becomes a random variable itself, with a variance which can be approximated by

$$\text{Var}\{\text{Var}\{\text{S}\}\} \approx \frac{1}{N\_{\text{lev}}} \left(\mu\_4 - \sigma^2\right) = \frac{1}{N\_{\text{lev}}} \left(\langle\mathcal{S}^4\rangle - 4\langle\mathcal{S}^3\rangle + 8\langle\mathcal{S}^2\rangle - \langle\mathcal{S}^2\rangle^2 - 4\right),\tag{34}$$

where *μ*<sup>4</sup> <sup>=</sup> ( *<sup>S</sup>* <sup>−</sup>*S*) <sup>4</sup> denotes the fourth central moment and we have used the normalization *S* =1. In turn, for spacings following the distribution *Pα*¯,*κ*¯(*λ*, *S*) (17) one can find they do not follow GOE if

$$N\_{\rm lev} \gtrsim \left[\frac{3\sqrt{(\mu\_4 - \sigma^2)\_\lambda}}{\text{Var}\,\{S\}\_\lambda - \text{Var}\,\{S\}\_{\rm OOE}}\right]^2,\tag{35}$$

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

$$E\_{\text{max}}\sqrt{N\_{\text{tot}}} \gtrsim 9.9 \text{ t} \times \frac{\sqrt{(\mu\_4 - \sigma^2)\_\lambda}}{\text{Var}\,\{\mathcal{S}\}\_\lambda - 0.273}}. \tag{36}$$

<sup>4</sup> We use the property of *m*-th cumulant of the distribution *P*(*S*)= <sup>1</sup> <sup>2</sup> *aP*1(*aS*) + *bP*2(*bS*) , which is equal to *<sup>S</sup> <sup>m</sup> <sup>P</sup>* =1 / 2(*<sup>a</sup>* (−*m*) *<sup>S</sup> <sup>m</sup>* (*P*1) <sup>+</sup> *<sup>b</sup>* (−*m*) *<sup>S</sup> <sup>m</sup>* (*P*2)).. For *P*<sup>1</sup> <sup>=</sup>*P*GOE and *P*<sup>2</sup> <sup>=</sup>*P*GOE−GUE, see Eqs. (6) and (3), necessary integrals for *m*=2, 3, 4 can be calculated analytically.

Nonstandard Transition GUE‐GOE for Random Matrices and Spectral Statistics of Graphene Nanoflakes 105 15 http://dx.doi.org/10.5772/64240

For *<sup>λ</sup>* <sup>→</sup>0, we have *μ*<sup>4</sup> <sup>−</sup>*<sup>σ</sup>* <sup>2</sup><sup>→</sup> <sup>21</sup> <sup>16</sup> *<sup>π</sup>* <sup>2</sup> <sup>−</sup>2*<sup>π</sup>* <sup>−</sup>4≈2.671, leading to

(31)

(32)

(33)

(34)

(35)

(36)

<sup>2</sup> *aP*1(*aS*) + *bP*2(*bS*) , which is equal to

Physically, occupying *N*lev electronic levels above the Dirac point one produces the electric charge *Q* = −2*eN*lev, resulting in a typical experimental limit of *E*max =0.2−0.3 eV for graphene

Level‐spacing distributions *P*(*S*) are normalized such that . In turn, the variance

where Eq. (33) refers to the empirical distribution *Pα*¯,*κ*¯(*λ*, *S*) given by Eq. (17) with *λ* →0. Similar calculation for arbitrary *λ* is straightforward, but the resulting formula is too lengthy to be

where *μ*<sup>4</sup> <sup>=</sup> ( *<sup>S</sup>* <sup>−</sup>*S*) <sup>4</sup> denotes the fourth central moment and we have used the normalization *S* =1. In turn, for spacings following the distribution *Pα*¯,*κ*¯(*λ*, *S*) (17) one can find they do not

where the factor 3 in the nominator corresponds to the 3*σ* level of significance. Substituting

*<sup>S</sup> <sup>m</sup> <sup>P</sup>* =1 / 2(*<sup>a</sup>* (−*m*) *<sup>S</sup> <sup>m</sup>* (*P*1) <sup>+</sup> *<sup>b</sup>* (−*m*) *<sup>S</sup> <sup>m</sup>* (*<sup>P</sup>*2)).. For *P*<sup>1</sup> <sup>=</sup>*P*GOE and *P*<sup>2</sup> <sup>=</sup>*P*GOE−GUE, see Eqs. (6) and (3), necessary

it becomes a random variable itself, with a variance which can be approximated by

When Var{*S*} is calculates for a large but finite collection of spacings *N*spc= *N*lev −1,

raises as the lowest moment allowing one to distinguish between

nanostructures on SiO2‐based substrates [49].

10414 Recent Advances in Graphene Research

different distributions. In particular, we have

presented.<sup>4</sup>

follow GOE if

4

Eq. (31) one can rewrite the above as

We use the property of *m*-th cumulant of the distribution *P*(*S*)= <sup>1</sup>

integrals for *m*=2, 3, 4 can be calculated analytically.

$$N\_{\rm tot} \gtrsim 223 \times \left(\frac{t}{\varepsilon\_{\rm max}}\right)^2 \qquad \text{for} \quad N\_{\rm vac} = 0. \tag{37}$$

**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*vac=0 splits the region where number of available en‐ ergy 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*vac=1, calculated numerically from Eq. (13) for *λ* =*λ*fit(*N*tot) (see Eq. (25)). Vertical red line in right panel marks the limit given by Eq. (27), above which the orthogonal class with grad‐ ual degeneracy splitting (2 × GOE → GOE) appears.

Limiting values of *N*tot and *E*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*vac=0), the attainable Fermi energy *E*max =0.25 eV should make it possible to detect TRS breaking in nanoflakes containing *N*tot≳<sup>3</sup> <sup>10</sup><sup>4</sup> carbon atoms, corresponding to the physical diameter of . For *N*vac=1, the limit of *N*tot≲9500 (see Eq. (27)) implies *E*max≳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 comple‐ mented 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*max (in typical experiment, the Fermi energy is tuned in the range −*E*max <*E* <*E*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.
