Splitting of critical energies in the $n$=0 Landau level of graphene

The lifting of the degeneracy of the states from the graphene $n$=0 Landau level (LL) is investigated through a non-interacting tight-binding model with random hoppings. A disorder-driven splitting of two bands and of two critical energies is observed by means of density of states and participation ratio calculations. The analysis of the probability densities of the states within the $n$=0 LL provides some insights into the interplay of lattice and disorder effects on the splitting process. An uneven spatial distribution of the wave function amplitudes between the two graphene sublattices is found for the states in between the two split peaks. It is shown that as the splitting is increased (linear increasing with disorder and square root increasing with magnetic field), the two split levels also get increasingly broadened, in such a way that the proportion of the overlapped states keeps approximately constant for a wide range of disorder or magnetic field variation.


Abstract.
The lifting of the degeneracy of the states from the graphene n=0 Landau level (LL) is investigated through a non-interacting tight-binding model with random hoppings. A disorder-driven splitting of two bands and of two critical energies is observed by means of density of states and participation ratio calculations. The analysis of the probability densities of the states within the n=0 LL provides some insights into the interplay of lattice and disorder effects on the splitting process. An uneven spatial distribution of the wave function amplitudes between the two graphene sublattices is found for the states in between the two split peaks. It is shown that as the splitting is increased (linear increasing with disorder and square root increasing with magnetic field), the two split levels also get increasingly broadened, in such a way that the proportion of the overlapped states keeps approximately constant for a wide range of disorder or magnetic field variation.

Introduction
Since the first observation of the anomalous integer quantum Hall effect in graphene [1,2], there has been an enormous interest in the quantum Hall regime of graphene. The sequence of quantized Hall plateaus in graphene is shifted by half-integer if compared to the well-known quantum Hall effect in usual two-dimensional electron gases. This is understood to be due to the unique features of the Dirac-like band structure of graphene around Fermi energy [3], which give origin to the following Landau level (LL) quantization: E n = ±ν F 2ehB|n|, where n is the Landau level index and ν F is the Fermi velocity. In this way, the energy dependence for the LLs with magnetic field B is not linear anymore, but goes with the square root of B (positive energies for electrons and negative energies for holes), and there is also a LL at zero-energy (the n=0 LL), which is shared by electrons and holes. In addition to the spin degree of freedom, graphene LLs exhibit a valley (sublattice) degree of freedom, giving a fourfold degeneracy to each LL.
Of special interest in the recent literature is the question of the lifting of the n = 0 LL degeneracies, observed experimentally [4,5,6,7,8] and discussed in many theoretical works [9][10][11][12][13][14][15][16][17][18][19][20][21][22]. Amongst these works, there are various different proposed explanations for the origins of the observed splittings, however a consensus is still missing. Two of the experimental works [4,5] show a complete lifting of the n=0 fourfold degeneracy occurring for high magnetic fields up to 45T and low temperatures, and there are evidences that the extra quantum Hall plateaus observed at filing factors ν = 0 and ν = ±1 should be due to the lifting of the spin and sublattice degeneracy, respectively [5]. The other three experiments [6,7,8] show the lifting of only one of the degeneracies. The authors from ref. [8] suggest that the opening of the gap in the n=0 observed by them might be spin related (Zeeman splitting), because of a linear variation of the gap with B, however, in ref. [7] the authors believe that the gap they measure should be driven by a broken sublattice symmetry. The physical mechanisms leading to the sublattice symmetry breaking are still not completely understood.
In this work, a splitting of two critical energies in the n=0 LL is observed within a nearest-neighbor tight-binding model for the graphene lattice, considering random hoppings as the only source of disorder. The splitting is inferred by means of participation ratio calculations, which show clearly the presence of two split peaks indicating two critical energies inside the n=0 Landau band, with localized states in between the peaks, around E=0. The splitting is found to have a linear dependence with the disorder strength and a square root dependence with magnetic field, in agreement with the numerical calculations from ref. [17] and with the magnetic field dependence observed for the gap in ref. [5]. As the model used in this work considers a single particle picture (no electron-electron interaction is taken into account) and there are no spin-dependent terms in the Hamiltonian, the splitting observed is related only to the effects of the random hopping disorder, which promotes intervalley mixing. To help in the understanding of these effects, I look here to the wave functions (probability densities) of the states within the n=0 LL, inspecting in particular the distribution of amplitudes in each of the sublattices, which gives interesting information about the process of disorder-driven sublattice mixing.
A unit cell of the graphene hexagonal lattice contains two carbon atoms, defining two sublattices, A and B. For the perfect graphene lattice (in absence of any kind of disorder), the wave functions of all the states within the n=0 LL have non-zero amplitudes in only one of the sublattices. "How does disorder affect this picture?" was the question addressed in Ref. [23], where it is shown how the increasing inclusion of a diagonal disorder (both white-noise or a Gaussian-correlated disorder model) increasingly promotes the spread of the wave functions to both sublattices in the n=0 LL, a effect that revealed to be directly connected to an anomalous localization: an enhancement of the participation ratio of the states with increasing disorder [23,24]. However, for diagonal disorder models (on-site energy fluctuations) the splitting of critical energies does not occur. Here the same question is addressed for a off-diagonal disorder model, for which a completely different picture emerges. The probability densities of the states from the n=0 LL are not anymore concentrated predominantly in one of the the sublattices, even when small disorders are considered. However, the interesting find is that for the states closest to E=0 (which get localized for the present disorder model, contributing to the critical energies splitting process), the spatial distribution of the wave functions in one sublattice is completely different and independent of the distribution in the other sublattice. These are states corresponding to the region where the two split bands have an important overlapping.

Model: Tight-Binding Hamiltonian with Random Hopping
The tight-binding Hamiltonian model for the honeycomb lattice used in this work is defined by: where c i is the fermionic operator on site i and the sum runs over nearest-neighbor bonds. The hopping parameters t ij randomly fluctuate (white-noise) from bond to bond around the average value t: Thus W represents the width of the uniform off-diagonal disorder distribution. t≈2.7eV for graphene, however note that the results shown here are parameterized by t. The graphene lattices considered have M ×N atoms (M zig-zag chains, each containing N atoms), which are repeated by means of periodic boundary conditions. In this way, the dimensions of the lattices that are periodically repeated are given by: where a=2.46Å is the lattice constant for graphene. The perpendicular magnetic field B is included by means of a Peierls' substitution (a complex phase factor in the hopping parameter): φ ij = 2π(e/h) i j A · dl . In the Landau gauge, φ ij = 0 along the x direction and φ ij = ±π(x/a)Φ/Φ 0 along the ∓y direction, with Φ/Φ 0 = Ba 2 √ 3e/(2h). The focus is on the low-flux limit (Φ/Φ 0 < 0.05), where the graphene LLs are well defined in this model.

Density of States and the n=0 LL Splitting
In Figure 1  defining the disordered unit cell that is periodically repeated (periodic boundary conditions). The spectrum is symmetric about zero energy, thus the n < 0 Landau levels are not shown here. The DOS is an energy histogram accumulated after the Hamiltonian diagonalization for hundreds of disorder realizations -in the case of Fig.1, 360 realizations were taken. Here we see that the n=0 LL has about the same broadening width of the higher levels, which is different from the result obtained in Ref. [17] for a random magnetic flux disorder. However, the dependence of the n=0 LL broadening with the correlation length of the random hopping was already elucidated in Ref. [20], which shows a DOS with broadenings in accordance to the results shown here for spatially uncorrelated random bonds. Fig.1(b) shows a zoom of the DOS for the region of the n=0 LL, allowing the observation in more details of the shape of this central band, which has many differences when compared to the n=0 Landau band shape obtained for a diagonal disorder model [23,24] and even when compared to the higher levels. One can see in Fig1(b) that the broadening of this central LL does not resemble the usual Gaussian or semielliptic LL line shape, but instead, has a lowering in the density of states around the center (apart from the presence of a small peak observed in the DOS for the states closer to E = 0). It is important to note that these characteristics for the n=0 LL are kept for many other values of disorder and magnetic fluxes investigated in this work.
The localization properties of the states within the n=0 LL are inspected, as shown in Fig.1(c), throughout the calculation of the participation ratio (PR) [25], defined by: where ψ i is the amplitude of the normalized wave function on site i, and N = M × N is the total number of lattice sites. The PR gives therefore the proportion of the lattice sites over which the wave function is spread. In this way, the PR for a localized state vanishes in the thermodynamic limit, while peaks in the PR indicate the presence of extended states (critical energies). Fig.1(c) allows the observation of two peaks in the PR, clearly indicating the splitting of two critical energies. Note that this splitting occurs exclusively for the n=0 LL in the considered model: the PR calculated for the higher levels (not shown here) shows always only one critical energy (one PR peak) per Landau band. It is also important to point out that the splitting of the two critical energies is well defined even considering that there is no clear opening of an energy gap in the DOS (see section 6 for further discussion on this). As observed in the comparison between Fig.1(b) and Fig.1(c), there are states contributing to the DOS in the energy range between the two PR peaks, however, the closer these states are to E=0, the more localized they are, i.e., the smaller is their localization length.
Throughout this work, various PR calculations are shown, corresponding to different system sizes, different values of magnetic flux and different disorder strengths. In order to examine the question of possible finite lattice-size effects on the participation ratio calculations and get a better understanding of the PR peaks, in Fig.2 the PR is shown for three different lattice sizes. The dimension L x =19.8nm is kept constant while L y varies from L y =3.3nm, to L y =4.8nm and to L y =11.4nm (more details are specified in the figure caption). Disorder and magnetic fluxes values are the same for the three cases: W/t=0.1 and Φ/Φ 0 =2/94≈0.021. The magnetic length for this flux is l B = h/(eB)=6.3Å. Fig.2 shows how the peaks in the PR are better defined as the system size is increased. As expected, the PR values corresponding to the localized states (at the LL tails and LL center) decrease with increasing system size [25]. The periodical boundary conditions impose restrictions over the values of magnetic flux in order to guarantee a commensurability between the wave function periodicity (the phases in the hopping) and the periodicity of the lattice considered. For a system with M ×N sites, and considering the Landau gauge, the possible sequence of flux values is: Φ/Φ 0 = 2/M, 4/M, 6/M, .... Consequently, to consider small values of magnetic flux it is necessary to increase M (i.e., increase the L x dimension), however, computational limitations restrict the maximum matrix sizes to be diagonalized. To analyze small fluxes it is then usually necessary to consider rectangular lattices, which are bigger in L x than in L y . However, we have to keep in mind that L y should still be many times greater than l B in order to avoid finite size effects and to the PR peak to be well defined, as Fig.2 suggests. It is worth to note that the increase of N in Fig.2 determines also the increase of the number of states per Landau level. The degeneracy is given by d = M ×N ×Φ/Φ 0 and, due to the flux restrictions discussed above, for a flux Φ/Φ 0 = 2/M , we have d = 2N states per LL.

Probability Densities in each of the Sublattices
In order to try to get more information about the states within the n=0 LL and, consequently, about the origin of the critical energies splitting observed, the wave functions of these states are now examined. Special attention is paid to how the wave  and |Ψ B | 2 ), for three selected states indicated by the arrows. While for the states 2 and 3 the probability density has the very same spatial distribution in each of the two sublattices, the state 1 is peculiar because it has completely different and independent probability densities in each sublattice (all the states close to E=0, from the region inside the circle, show this behavior). function amplitudes are distributed in each of the two sublattices, once in absence of disorder all the eigenstates within the n=0 LL have non-zero amplitudes only on one of the sublattices [9].
In Fig. 3 it is shown the participation ratio calculated for a graphene lattice with M ×N = 60×108 atoms (L x ×L y = 12.6nm×13.2nm), considering a magnetic flux Φ/Φ 0 = 2/60 ≈ 0.033 and a random hopping amplitude W/t=0.1 averaged over 500 disorder realizations. It is also shown, for three selected states indicated by the arrows, the wave function probability densities, |Ψ| 2 , plotted over all the 60×108 lattice sites, and then, bellow, the probability densities over each sublattice are plotted separately. After calculating the eigenvectors from the exact Hamiltonian diagonalization, we have the probability density |Ψ (i,j) | 2 for all the lattice sites (i, j). Then, the decomposition shown consists simply of plotting separately the probability densities for only sites (i, j) which belong to one of the sublattices: A (|Ψ A | 2 ) or B (|Ψ B | 2 ). Each one of the three selected states corresponds to a different characteristic region from the n=0 LL: state 1 is a localized state very close to E=0, state 2 is a state from the region of critical (extended) states at the PR peak, and state 3 is a localized state from the band tail.
The most interesting observation refers to the state 1, for which the spatial distribution of the probability density on one sublattice is completely different from that on the other one. Observing |Ψ A | 2 and |Ψ B | 2 separately, after the sublattice decomposition, one can see that for the state 1 the distribution of maxima and minima amplitudes on one sublattice occurs in different positions and is completely independent from the distribution on the other. All the states in the region delimited by the dashed circle, states close to E=0, have this absolutely uneven spatial distribution of amplitudes between the sublattices. Observe that this is the region where the two Landau bands (split after the symmetry breaking) have a more significant overlap. For the states in the region right after the limits of the circle, some position correlation between the two sublattices starts to appear in the probability density (not shown here). But when reaching the proximities of the PR peak, and after that until the band tail, we observe that all the states show a symmetric distribution between the sublattices, as can be observed from the decompositions of states 2 and 3: for these states, the amplitudes on one sublattice looks pretty much like the amplitudes on the other.
It is worth to note that although having distributions that are asymmetric between the sublattices, it is not the case that the wave function is more concentrated in one sublattice than in the other. The behavior observed here is different from the observed when the disorder is introduced through on-site energy fluctuations (diagonal disorder models), for which the amplitudes over one sublattice are more significant than over the other [23,24]. In fact, both sublattices have equally significant amplitudes (note that the vertical scales of the graphics are kept the same in the decomposition): the sum of all amplitudes over the sublattice A was calculated and is pratically equal (with differences smaller than 5%) to the sum over the sublattice B.
Observing the wavefunction amplitudes in Fig. 3 one sees that the localization length of state 1 is larger than the one of state 3. Nevertheless, this difference reflects exactly the difference in the PR values of these states: for this system size, the PR for state 1 is not as small as for state 3 (remember however that the lattice size analysis shown in Fig. 2 indicated that both state 1 and state 3 are truly localized states). The reason for the larger localization length for the state 1 compared to state 3 is most probably due to the mixing of states in this region between the PR peaks, where the split LLs are overlapped (see section 6).
Another point to be noted is that, for the disorder model considered (only offdiagonal, with random nearest-neighbor hoppings), there is a perfect symmetry around E=0: the probability density of the state at an specific energy E is exactly equal to that of the state at −E. One can also observe this through the perfect symmetry of the participation ratio curve around E = 0.

Splitting as a function of Disorder and Magnetic Field
The evolution of the critical energies splitting in the n=0 LL with increasing disorder strength and magnetic field is investigated here, throughout the results shown in Figures  4 and 5, respectively. Fig.4(a) shows the PR of the states from the n=0 LL, for three different disorder strengths: W/t=0.2, W/t=0.4 and W/t=0.6. One can see that the energy splitting ∆E/t between the two critical energies (between the two PR peaks) is increased as the disorder increases the Landau level broadening. In Fig.4  the data for the wide interval of disorder considered: the red line passes through origin and has angular coefficient 0.1306 ±0.0006, which is valid for the specific magnetic flux considered here ( Φ/Φ 0 =2/66≈0.030): The results shown in Fig.4 are calculated for lattices having size 13.9nm×9.2nm (M ×N =66×76) and for a magnetic flux Φ/Φ 0 =2/66≈0.030, averaging the results for each W/t over 400 disorder realizations.
In Fig.5 the dependence of the splitting with magnetic field is examined. Fig.5(a) shows four examples of PR calculations within the n=0 LL, for fixed disorder (W/t = 0.1) and different magnetic fluxes. An increasing energy splitting of the two PR peaks is observed with increasing magnetic flux. In Fig.5(b) this splitting is plotted as a function os magnetic flux for several magnetic fluxes, and we observe that the data is well adjusted by a square root fit (the specific coefficient value is valid for the specific disorder W/t=0.1 considered): The dimensions of the lattices considered for each flux, as well as the magnetic lengths and number of disorder realizations undertaken are listed in table 1.
Regarding the square root dependence of the splitting with magnetic field, observe that the critical energies from the n=0 LL follow then the same dependence with B presented by the higher Landau levels. Table 1. Lattice dimensions, magnetic length l B and number of disorder realizations, for each of the magnetic fluxes considered to define the data points in Fig.5. It is important to note that these results are in agreement with reference [17], where a linear dependence of the energy splitting with the disorder, and a square root dependence with the magnetic field were observed for a random magnetic-field disorder model, and where the splitting was determined by means of the two-terminal conductance peaks.

Interplay between Splitting and LL Broadening
It is observed in Fig.1 that the splitting is much better defined trough the localization properties of the states (definition of two clearly split PR peaks) than trough the density of states. Nevertheless, the DOS shape observed in Fig.1(b) suggests a superposition of two split bands. Indeed, using a multi-peak fitting procedure, shown in Fig.6 (a) and (b), it is found that the DOS obtained are reasonably well fitted by two overlapping Gaussian curves (apart from the small peak at E = 0, whose origin is not understood). In Fig.6(c) it is possible to compare overall broadenings of the two n=0 LL DOS, which are calculated for two different values of disorder, W/t=0.05 and W/t=0.6, keeping constant all the other parameters: lattice size is M ×N =66×76 and the magnetic flux considered is Φ/Φ 0 =2/66=0.0303. Fig.6 (d) shows that two Gaussian peaks fitted to the DOS have a good coincidence with the calculated PR peaks positions.
The splitting is observed in the previous section, through the separations of the two PR peaks, to have a square root dependence with magnetic field, and to increases linearly with the disorder. However, what is important to discuss now is that, even when the splitting increases, the overlapping between the two degeneracy broken levels shows to be pretty much about the same, for a wide range of parameters observed, as the comparison between Fig.6(a) and Fig.6(b) indicates. This is because the same increasing disorder or magnetic field that causes the increasing of the critical energies splitting, also causes the increasing broadening of each Landau band. In other words, we would need that the rate on which the splitting increases with disorder and with magnetic field were much higher than the rate on which the LLs get broadened by them, to be able to observe an increasingly well defined separation of the DOS in two independent (not overlapped) bands. The interplay between both of these effects (splitting and broadening of the bands) produces in the end an overall DOS for the n=0 LL having always the same appearance (the difference being only the broadening, but the main features are kept the same). This indicates a direct connection between the effects of LL broadening and LL splitting. To the experiments, one consequence of this behavior is that it might be equally difficult to observe a real energy gap throughout DOS measurements, independent of the amount of disorder or of the intensity of the magnetic field. The random hoppings considered here are spatially uncorrelated. A further desirable extension for this work is to consider a spatial correlation for the hoppings, a situation that closer emulates the ripples in real graphene sheets. In Ref. [20] it is observed that the introduction of an increasingly higher correlation length for the random hoppings in graphene results in an increasingly thinner n=0 LL compared to the broadening of the higher LLs. On the other hand, the splitting of two critical energies within the n=0 LL is observed in Ref. [17] for a random flux model, for which the n=0 LL is much thinner than the higher levels (even being thinner, the DOS shows a n=0 LL shape very similar to the observed here). The results of these two works [17,20], together with the results shown here, suggest that the splittings in the central LL should survive when a finite correlation length for the hoppings is taken into account in the model, however with the energy splittings being scaled down with the width of the Landau band, in a similar fashion to the observed in Fig. 6.

Conclusions
The breaking of the degeneracy of the graphene n=0 LL is investigated in this work through a simple numerical model, considering a non-interacting tight-binding model for a hexagonal lattice with random nearest-neighbors hoppings. Inferring the localization properties of the states by means of participation ratio calculations, two clearly split critical energies are observed. The origin for this splitting has to be intrinsic to the disorder model considered, which introduces inter-valley mixing to the system.
The energy splitting is observed to have a linear dependence with the disorder strength, and a square root dependence with the magnetic field, confirming the results obtained in ref. [17]. The analysis of the probability densities of the wave functions within the n=0 LL shows that there is a region, for the states closer to E=0, where there is an important asymmetry in the distribution of the wave function amplitudes between the two sublattices. In this region, although there are amplitudes over both sublattices, there is no matching of the spatial positions of the amplitudes over each sublattice. It is like as each sublattice has its own probability density, completely independent from the other. This may lead to a diminish in the state percolation over the lattice, increasing even further the localized character of these states and also influencing on anomalous localization properties. It is also observed that this region of states with special wave function distribution between the sublattices coincides with the region where the two degeneracy broken Landau levels have a more important overlapping. Another interesting observation is that this overlapping keeps approximately constant for different disorder strengths, due to the interplay between LL broadenings and LL splitting with increasing disorder or magnetic flux. These results may help in elucidating the physics involved in the splitting of the n=0 LL due to valley/sublattice symmetry breaking.