Perturbative/nonperturbative aspects of Bloch electrons in a honeycomb lattice

We revisit the spectral problem for Bloch electrons in a two-dimensional bipartite honeycomb lattice under a uniform magnetic field. It is well-known that such a honeycomb structure is realized in graphene. We present a systematic framework to compute the perturbative magnetic flux expansions near two distinct band edges. We then analyze the nonperturbative bandwidth of the spectrum. It turns out that there is a novel similarity between the spectrum near the Dirac point in the honeycomb lattice and the spectrum in the supersymmetric sine-Gordon quantum mechanics. We finally confirm a nontrivial vacuum-instanton-bion threesome relationship. Our analysis heavily relies on numerical experiments.


Introduction
In his seminal paper [1], Hofstadter considered a simple two-dimensional lattice model for Bloch electrons in the presence of a uniform magnetic field. In spite of the simple physical setup, its electron spectrum turned out to be remarkably rich. Now, the Hofstadter model is recognized to be significant both in theoretical physics and in experimental physics (and also in mathematical physics). Recently, a new aspect of this model was revealed in the context of topological string theory on toric Calabi-Yau threefolds [2], based on [3] that is far from condensed matter physics originally. In [4], the correspondence was generalized to the triangular lattice system [5] and another Calabi-Yau geometry.
Motivated by these developments, we here revisit a more interesting model, Bloch electrons on a hexagonal honeycomb lattice. As is well-known, graphene has such a honeycomb structure, and thus the honeycomb lattice is one of the most exciting 2d electron systems. A characteristic property of the honeycomb lattice is that it has a zero energy gap. Around the zero-gap energy, the dispersion behaves linearly, and it forms the Dirac cone, where electrons are effectively described by massless Dirac fermions. Such singular points (referred to as Dirac points below) make graphene a quite rich matter. The spectrum of Bloch electrons in the honeycomb lattice under a magnetic field was studied by Rammal in [6].
In this paper, we focus on perturbative and nonperturbative aspects of the honeycomb lattice. A relation to quantum geometry in a toric Calabi-Yau threefold will be studied in detail in [7]. We here analyze the spectrum in the weak magnetic flux regime. It is well-known that in this regime, the spectrum is characterized by Landau levels. In the case of the honeycomb lattice, there are two energy regions in which we can use perturbation theory. One is around the Dirac point, E = 0. This is the bottom (or top) of the upper (or lower) band corresponding to the zero-gap in the zero magnetic field limit. We can also consider the perturbation around the other band edge at the top. In our convention, this corresponds to E = 3 (or E = −3 equivalently). By turning on the magnetic field, the spectrum is approximated by the perturbative expansion of the (rescaled) magnetic flux φ := 2πΦ/Φ 0 , where Φ 0 = hc/e is the magnetic flux unit and Φ is the magnetic flux per unit cell. We find that the perturbative expansion around the top E = 3 is given by On the other hand, the perturbative expansion near the Dirac point is more involved. We find the following expansion: In these equations, the non-negative integer n labels the Landau level of the spectrum. Note that the singular behavior E ∼ ±3 1/4 √ nφ at the Dirac point was first found by McClure long time ago [8]. We apply an idea by Bender and Wu in [9] to the Bloch electron system, and compute these perturbative expansions systematically. See also [10] on the similar approach. Such an efficient way is helpful to study nonperturbative physics in the sprit of "resurgence theory".
One notices that the lowest Landau level n = 0 at the Dirac point is very special. Its perturbative expansion seems vanishing, E pert Dirac (0, φ) = 0, (1.3) to all orders in perturbation theory. In fact, one can check that E = 0 is always one of band edges for any rational Φ/Φ 0 . The band edge E = 0 does not receive any "quantum corrections" in φ. It is protected due to a certain symmetry like supersymmetry. 1 We are interested in nonperturbative corrections to the perturbative expansions above. A good quantity to explore them is the bandwidth of the spectrum since the bandwidth is purely nonperturbative in φ, related to quantum tunneling effects. Based on the numerical analysis, we observe that the bandwidth for φ = 2π/Q with integer Q actually scales as where the constant A = 10.149 · · · is exactly given by an integral form (3.18). Interestingly, we find that the spectrum near the Dirac point is very similar to that in the supersymmetric sine-Gordon quantum mechanics [12,13]. 2 This suggests that "Cheshire Cat Resurgence" in [13] is useful in exploring more details on the nonperturbative structure in the weak flux regime. Finally, we check the following threesome relation among the fluctuations around the vacuum, one-instanton and bion (i.e., instanton-anti-instanton) saddles: where n is a quantum number of the spectrum, and these fluctuations are normalized appropriately. This is naturally expected by the so-called "perturbative/nonperturbative relation," originally found in quantum mechanics [14][15][16][17]. We would like to emphasize that the relation (1.5) holds universally in a wide class of quantum mechanical systems including 2d Bloch electron systems, whose resurgent property has not yet been understood well. 3 The organization of this paper is as follows. In the next section, we start by reviewing Bloch electrons in a bipartite honeycomb lattice. We explain how to compute the spectrum if the magnetic field is turned on. Section 3 is a main part of this paper. We consider the weak magnetic flux limit, and present a systematic method to compute the weak magnetic flux expansions near E = 3 and E = 0. We also discuss nonperturbative corrections in the bandwidth, based on numerical experiments. We check a perturbative-instanton-bion triangle relation. In section 4, we give comments on nonperturbative corrections in the supersymmetric sine-Gordon quantum mechanics. In section 5, we give final remarks. In appendix A, we explain our numerical technique used in this paper.

Bloch electrons in a honeycomb lattice
We will start with a review of an electron system in a two-dimensional honeycomb lattice. We fix our convention by following [6].

No magnetic flux
In this subsection, we consider the case of no magnetic field. The honeycomb lattice we are now treating is a bipartite system with two sublattices. We refer to them as A and B. See the left of figure 1. The primitive translation vectors are where a is a lattice spacing constant. For later convenience, we also introduce three vectors in the real space, as shown in figure 1 (left), by In the reciprocal space, the primitive vectors are given by 3) The first Brillouin zone is a hexagon. Two corners of the hexagon are usually denoted by K and K , whose wave vectors are given by See the right of figure 1. Let us proceed to the dispersion relation. In the tight-binding approximation, the wave function with wave vector k is written as where R A = n 1 a 1 + n 2 a 2 + d 1 and R B = n 1 a 1 + n 2 a 2 ((n 1 , n 2 ) ∈ Z 2 ) label positions of sublattices A and B respectively. If considering only the nearest neighboring hopping, one obtains the following eigenvalue equations: where D(k) := e −ik·d 1 + e −ik·d 2 + e −ik·d 3 . Therefore the eigenvalues are given by E(k) = ±|D(k)|. Plugging (2.2), we get the dispersion relation (2.7) The range of the allowed energy is −3 ≤ E(k) ≤ 3. Obviously, it is symmetric under E → −E. For a given wave vector k, we have two eigenvalues E(k) = ±|D(k)|, but at the points K and K , these coincide, since E(K) = E(K ) = 0. Therefore these points are   zero-gap points. The dispersion is shown in figure 2 (left). Since the dispersion near the zero-gap point K behaves as the energy surface forms two cones around this point, as in figure 2 (right). These are called the Dirac cones, and we refer to the points K and K as the Dirac points, near which electrons are effectively described by relativistic massless fermions. Finally, the density of states is given by where m = (2|E|) −1 (1 + |E|) 3 (3 − |E|) and K(m) is the complete elliptic integral of first kind. We plot ρ(E) in figure 3. The density of states turns out to diverge at E = ±1 (a.k.a. Van Hove singularities). Also it behaves as ρ(E) = |E|/( √ 3π) + O(|E| 3 ) near E = 0.

Turning on magnetic flux
If the magnetic field is turned on, the spectrum of the electron becomes richer and more involved. In this case, one can effectively replace the wave vector by (2.10) Then the functions D(k) and D(k) * are replaced by difference operators. In the Landau gauge A = (0, Bx, 0), these difference operators take the form of where the phase factors e ±iφ/12 come from the Baker-Campbell-Hausdorff formula, and Since the area of the hexagonal unit cell is S unit cell = 3 √ 3a 2 /2, Φ measures the magnetic flux per unit cell. The eigenvalue equation (2.6) is now replaced by the following twodimensional difference equations: (2.13) In the y-direction, we can take the plane wave solution by Ψ X (x, y) = e ikyy ψ X (x), and the problem reduces to the following one-dimensional problem: (2.14) By eliminating the unknown function ψ A (x), one gets the difference equation for ψ B (x): . Then the above eigenvalue equation is finally reduced to so-called Harper's equation: where κ = − √ 3ak y /2. If φ/(2π) = Φ/Φ 0 is rational, the eigenvalue problem of Harper's equation (2.16) is reduced to the diagonalization of a finite dimensional matrix A. Let us assume φ = 2πa/b with coprime integers a and b. As shown in [6], the matrix elements of A are explicitly given by where θ 1 = 2κ and θ 2 is Bloch's angle due to the periodicity of Harper's equation. We can easily compute the eigenvalues of A for given angles θ 1 and θ 2 . These form finite b subbands of λ. We show the spectra of λ and E = ± √ λ + 3 for 0 ≤ φ ≤ 2π in figure 4. 4

Perturbative/nonperturbative aspects
In this section, we consider the weak magnetic flux limit φ ∼ 0 in the honeycomb system. As is well-known, in the weak flux limit, the spectrum is characterized by Landau levels.
We first discuss the perturbative expansion in φ, and then proceed to nonperturbative effects in the spectrum.

Perturbative expansions
To compute the perturbative expansion of the spectrum, we first rewrite the eigenvalue equation (2.15) as a more symmetric form. By shifting the argument The corresponding difference operator to the right hand side is After introducing a new variable q := φx/(3a), one obtains Recall that the eigenvalues of this operator gives λ = E 2 − 3 rather than E. Let us first consider the zero flux limit. The classical curve for (3.4) is λ = 4 cos q cos p 2 + 2 cos 2q. (3.5) Therefore λ takes the maximal value λ = 6 at p = q = 0 and the minimal value λ = −3 at p = 0, q = 2π/3. Here λ = 6 corresponds to the top of the band E = 3, while λ = −3 to the zero-gap energy E = 0. Around these points, H can be expanded as a perturbation of the harmonic oscillator. We further rescale the variables by q → gq and p → gp, where g := √ φ. We obtain (3.6) Expanding it around g = 0, we find This is the perturbative expansion around the top of the band. The coefficient at order g 2 is just the Hamiltonian of the harmonic oscillator with = m = 1 and ω = 2 √ 3. Therefore up to this order, the spectrum should be given by where the quantum number n represents the Landau level. Since the Hamiltonian (3.7) can be regarded as the perturbation of the harmonic oscillator, one can compute the higher order corrections by the textbook method. To compute it more efficiently, the idea of Bender and Wu in [9] is extremely useful. 5 According to [9], the normalizable perturbative wave function of the operator (3.7) takes the form of where P ( ) n (q) is a polynomial of q. For instance, the ground state wave function is given by with the perturbative spectrum: Assuming the Bender-Wu structure (3.9) to all orders in perturbation theory, one can determine P ( ) n (q) and the perturbative corrections to λ recursively. This enables us to put the computation on a computer. See [19,20] for example. Using this method, we find the perturbative expansion of λ: (3.12) Plugging this expansion into E = ± √ λ + 3, one obtains the expansion (1.1) of the electron energy.
To compute the expansion around λ = −3, we first shift the position by q → q − 2π/(3g), and then expand the Hamiltonian around g = 0. The result is In this case, we consider the perturbation of the harmonic oscillator with = 1, m = 2 and ω = √ 3. As in the same computation above, the perturbative expansion around the Dirac point results in (3.14) Note that in this expansion, the correction at order φ takes the form of n rather than n + 1/2. As a consequence, the electron energy behaves as E Dirac ∼ 3 1/4 √ nφ near the cone [8]. In figure 5, we show the Landau level splitting in the λ-spectrum. Our perturbative result shows good agreement with the weak magnetic spectrum. The perturbative expansion (3.14) is very similar to that in the supersymmetric sine-Gordon quantum mechanics [12,13] (see equation (4.5)). In particular, the lowest Landau level eigenvalue does not seem to receive any perturbative corrections: We have checked this property up to φ 60 . As we will see in the next subsections, the nonperturbative structure is also similar. Such a structure implies that the honeycomb lattice effectively has hidden supersymmetry near the Dirac point. This implication is consistent with an earlier work [11,21].

Nonperturbative bandwidth
In the previous subsection, we computed the perturbative expansions around two regimes λ = 6 and λ = −3. This is not the end of the story. The spectrum of Bloch electrons in the lattice system under the magnetic field has the subband structure rather than the discrete spectrum. The naive question is thus how wide these subbands are. This problem is not simple to answer because in the case of Harper's difference equation, the spectrum has a self-similar pattern. The subband structure becomes very complicated if the rational flux Φ/Φ 0 is complicated. In this paper, we restrict our analysis to the particular case of the flux with the form φ = 2π/Q, where Q is a large integer. We consider the weak flux limit Q → ∞, and observe a simple pattern of the width of subbands. As is expected, the bandwidth turns out to have an exponentially small scale e −αQ in the large Q limit due to quantum mechanical tunneling effects. Therefore the analysis of the bandwidth is a first step to understand   Figure 6. These figures show the bandwidths near λ = 6 (left) and near λ = −3 (right) against Q. In the right figure, there are two subbands whose bandwidths are almost same for each Landau level n ≥ 1. We denote these two by indices ±.
the nonperturbative effects. By using the numerical analysis explained in appendix A, we conjecture several analytic forms of the bandwidths.
As an example, we show in table 1 the numerical values of the edges of subbands of λ for Q = 30. In this case, there are 30 subbands, but we show only the five subbands near the top (λ = 6) and near the bottom (λ = −3), respectively, since we are interested in the bandwidth near these points. One can observe that these subbands are actually very narrow.
Let us first consider the bandwidth near the top. In figure 6 (left), we plot the bandwidths for the first five Landau levels in the range 10 ≤ Q ≤ 50. Clearly, these bandwidths are exponentially small. We can compute the exponential factor as follows. As is wellknown, the quantum mechanical tunneling effect is explained by the semiclassical WKB period integral for a classically forbidden orbit [22]. In our case, the classical curve of the system is given by (3.5). Then the period integral to explain the tunneling effect turns out to be where p (±) λ (q) are two branches in (3.5). By the careful choice of the numerical constant C to match the numerical result, we conclude that the correct exponential factor is where A is given by dq arccos 2 cos q − cos q = 10.149416064 · · · . (3.18) We also give a conjecture of the prefactor. By using the numerical fitting of the bandwidths for various Landau levels, we find the leading contribution where P inst top (n, 2π/Q) = 1 + O(Q −1 ) is the perturbative fluctuation around a nonperturbative saddle-point in the path integral perspective. A first few coefficients of P inst top (n, 2π/Q) will be conjectured in the next subsections. We expect that it is explained by the oneinstanton saddle. It is straightforward to rewrite this result into the width of the energy bands.
The analysis near the bottom edge is much more complicated. We first observe that for each Landau level n ≥ 1 there is a pair of subbands whose bandwidths are almost same, as shown in figure 6 (right). The gap of these two subbands is extremely narrow, and it is almost regarded as a zero-gap. See table 1 for instance. As will be seen in the next section, such a zero-gap structure also appears in the supersymmetric sine-Gordon system. We distinguish these two subbands by subscript ±. Then we again observe from numerics that each bandwidth of λ scales as ∆λ band Dirac,± (n, 2π/Q) ≈ where A is the same number given previously. The lowest Landau level n = 0 is special. As shown in figure 6 (right), the exponential decay for n = 0 is much faster than that for n ≥ 1. The careful numerical analysis reveals ∆λ band Since the exponential factor is the square of that for n ≥ 1, the fluctuation P bion Dirac should be explained by the two-instanton sector. By analogy with quantum mechanics, this contribution comes from the instanton-anti-instanton (i.e. bion) saddle. Let us translate this result into the energy spectrum. We use the fact that λ = −3 is always the bottom edge of the subband at the lowest Landau level. Hence the top edge of this subband is λ top edge Dirac (0, 2π/Q) = −3 + ∆λ band Dirac (0, 2π/Q). (3.22) The bandwidth of the lowest Landau energy is finally given by ∆E band Dirac (0, 2π/Q) = ∆λ band We conclude that the bandwidth of the energy for the lowest Landau level has the same order as those for n ≥ 1.

Testing a PNP threesome relation
The analysis in the previous subsection heavily relies on the numerical experiments. This is a kind of guess works. We would like to extract universal properties from these data.
Remarkably, recent resurgent analysis shows up new relations between the perturbative sector and nonperturbative sectors. Such a resurgent consideration also must tell us about relations among nonperturbative sectors. For the resurgent analysis in high energy physics, see [23][24][25] and references therein. Inspired by this development, we here check a simple relation among the perturbative, the one-instanton and the instanton-anti-instanton sectors, which is expected by the resurgent analysis in usual quantum mechanics. We find that the honeycomb lattice has the following beautiful relation where N is a normalization factor, which should be chosen as N Dirac = √ 3φ for the Dirac point λ Dirac = −3 and as N top = −2 √ 3φ for the top λ top = 6. For the reader who wonders why this relation is expected, see subsection 4.2. We have neither a rigorous proof nor strong evidence for this relation in the honeycomb system so far. Nevertheless we believe that (3.24) works both for the Dirac point and for the top. The reason is that this kind of relations seems to hold very widely not only for well-studied quantum mechanical systems but also for Hofstadter-type 2d electron systems [18].
We here give a few nontrivial checks, based on the numerical analysis. We focus on the Dirac point. As was seen before, the lowest Landau level, the bandwidth starts from the bion contribution. We observe that its fluctuation is given by Also the bandwidth for n ≥ 1 is explained by the one-instanton correction. The numerical experiment allows us to find log P inst Dirac (n, φ) = − 30n 2 + 72n + 11 72 √ 3 φ − 34n 3 + 96n 2 + 49n + 16 432 φ 2 − 4470n 4 + 17280n 3 + 14910n 2 + 12960n + 1081 58320 (3.26) where P inst Dirac = P inst Dirac,+ = P inst Dirac,− . To find this result, we first compute the expansion for various Landau levels n = 1, 2, . . . , and then determine each coefficient by assuming that the coefficient at order φ k is a polynomial of n with degree k + 1. This assumption is validated by checking whether the obtained result reproduces the ones for higher Landau levels or not. Note that the expansion (3.26) is obtained by the analysis for n ≥ 1, but we can extrapolate it to n = 0.
Combining these observations, one can evaluate both sides of (3.24) for n = 0, independently. The left hand side is The right hand side is Indeed these expansions perfectly agree! To check the relation for n ≥ 1, we need the bion correction. This correction can be extracted from the large order behavior of the perturbative expansion. Let us write the perturbative and the bion expansions as where our normalization is a (n) = 1. The resurgent analysis (see [23][24][25] for example) tells us that the information on a (1,1) (n) is encoded in the large order behavior of a (0) (n). The large order behavior has the following form in the limit → ∞. Here S = A/5 and C n is an n-dependent constant. We find C 1 = −81 √ 3 and C 2 = −2187 √ 3/2. We have the perturbative data a (0) (n) (n = 1, 2) up to = 100, and use them to estimate a (1,1) (n) (n = 1, 2) numerically by the method in appendix A. We then obtain the following numerical values of the bion corrections: (2) ≈ −1.4032762609, (3.31) Now, we compare these values with the prediction from the relation (3.24). Using the expansions of λ pert Dirac and P inst Dirac , one easily obtains the fluctuation around the bion: (3.32) The coefficients are actually in agreement with the numerical estimation (3.31). For the nonperturbative corrections around the top, it is not easy to check the relation because the bion fluctuation is hard to be evaluated. 6 Alternatively, we use (3.24) to predict the bion correction. For the top, we find the following one-instanton fluctuation (3.33) Our threesome relationship predicts the bion fluctuation: (3.34) It would be nice to check whether this prediction is correct or not.

Remarks on SUSY sine-Gordon QM
In this section, we give some remarks on nonperturbative effects of the supersymmetric sine-Gordon quantum mechanics. In [12,13,26], the instanton-anti-instanton correction was mainly studied in the context of resurgence theory. Here we look into the bandwidth of this model, which is not discussed in these references. The bandwidth is caused by one-instanton effect.

Band structure
We start with the following supersymmetric quantum mechanics: where σ i (i = 1, 2, 3) are the Pauli matrices, and W (x) is a superpotential. The first component corresponds to the fermionic sector, while the second component to the bosonic sector. We here consider the sine-Gordon potential: For our purpose, it is convenient to rescale the variable by x = √ gq. Then, the bosonic sector of the Hamiltonian is This Hamiltonian admits the expansion around g = 0, We can use the Bender-Wu method to compute the perturbative expansion of the energy.
Using the Mathematica package in [19], we obtain E pert (n, g) = n − n 2 4 g − n(2n 2 + 1) 32 g 2 − n 2 (5n 2 + 7) 128 g 3 − n(66n 4 + 182n 2 + 25) 2048 g 4 − n 2 (63n 4 + 290n 2 + 127) 2048 g 5 + O(g 6 ). (4.5) As was mentioned before, this perturbative expansion is very similar to (3.14). In particular, the expansion exactly vanishes for n = 0: Let us see the band structure. Since the supersymmetric sine-Gordon system has a 2π-periodic potential, the Floquet-Bloch theorem states that the wave function takes the form of where k is a continuous parameter, ranging from k = −1/2 to k = 1/2. Plugging this into the Schrödinger equation, we obtain an infinite set of relations By diagonalizing the left hand side of this equation, we obtain the discrete eigenvalues E k (m, g) (m = 0, 1, 2, . . . ) for a given wave number k and coupling g. Since the wave number moves in the first Brillouin zone (−1/2 ≤ k ≤ 1/2) continuously, the spectrum forms energy bands. We show the k-dependence of E k (g) for g = 1/2 in figure 7 (left). One characteristic feature of this model is that the upper band edge with band index 2m − 1 and the lower band edge with band index 2m coincide for all m ≥ 1. This means that there are an infinite number of zero-gap energies: E k=0 (2m − 1, g) = E k=0 (2m, g). For convenience we denote these bands as Another important fact is that E = 0 is always the lowest eigenvalue at k = 0: E k=0 (0, g) = 0. This is a reflection that the system has supersymmetry. The ground state does not receive any quantum corrections. In figure 7 (right), we plot the energy as a function of the coupling g.
Let us proceed to the computation of the bandwidth. We do it numerically. The bandwidth is given by ∆E band ± (n, g) := |E k=0,± (n, g) − E k=1/2,± (n, g)| For the lowest band n = 0, we also define ∆E band (0, g) := E k=1/2 (0, g). We first solve the eigenvalue equation (4.8) for various g numerically. We then fit these values under the ansatz The numerical fitting allows us to guess exact values of A, B and C. After that, we use the method in appendix A to determine c j 's. For the lowest band, we find ∆E band (0, g) ≈ 2 π e − 4 g P bion (0, g), (4.12) where P bion (0, g) is expected to come from the bion saddle in the path integral formulation. We observe that it has the following perturbative expansion (4.13) For n ≥ 1, we also find the leading correction where the fluctuation P inst (n, g) does not depend on two ±-bands, and it should be explained by the one-instanton saddle (and the one-anti-instanton saddle). We can guess a few terms log P inst (n, g) = − 6n 2 + 8n + 1 16 g − 10n 3 + 20n 2 + 7n + 2 64 g 2 − 330n 4 + 896n 3 + 546n 2 + 384n + 25 3072 g 3 + O(g 4 ). (4.15) Note that the bandwidth for n = 0 comes from the bion contribution O(e −4/g ), while for n ≥ 1 from the one-instanton contribution O(e −2/g ). This property is the same as the bandwidth of λ for the honeycomb lattice in section 3. The similar behavior is also found in the supersymmetric double well potential (or the Fokker-Planck model in the literature) [27].

PNP relations
There is a remarkable relationship between the perturbative sector and the nonperturbative sectors. We refer to it as the perturbative/nonperturbative (PNP) relation [14][15][16][17] (see also [28][29][30]). In our case, this relation results in the form where S inst = 2 is the instanton action. We note that the similar formula was found long time ago [31]. Since the right hand side in this equation contains only the perturbative quantity, the one-instanton correction is completely controlled by the perturbative sector! Using the perturbative expansion (4.5), one can check that it reproduces our guess (4.15) based on the numerical analysis. The authors in [12,13] conjectured that the fluctuation around the bion saddle also has a similar formula: This conjecture indeed reproduces (4.13) for n = 0. Also, the conjecture was confirmed by comparing it with the large order prediction [13]. Obviously, these two equations are almost the same form. After a simple computation, we arrive at P bion (n, g) P inst (n, g) 2 = ∂E pert (n, g) ∂n −1 .

(4.18)
This is just the formula we saw in the previous section.
In well-studied examples such as the double well potential and the cosine potential, this relation is not so surprising, because in these cases, we already know that the Zinn-Justin-Jentschura (ZJJ) exact quantization conditions [32,33] contain only the two non-trivial functions, a "perturbative" function B(E, g) and a "nonperturbative" function A(E, g). The complete trans-series expansion of the spectrum is expressed in terms of this two functions in principle. Moreover, these two functions A(E, g) and B(E, g) are exactly related by the PNP relation [14][15][16][17]. 7 On one hand, the formulae (4.16) and (4.17) are understood as a consequence of the PNP relation. On the other hand, the relation (4.18) is considered to be a direct consequence of the ZJJ exact quantization conditions. A remarkable point is that the relation (4.18) seems to hold universally. As will be reported in [18], it also works in the Hofstadter model [1]. In particular, we are not able to find the one-instanton formula like (4.16) for the Hofstadter model nor for the honeycomb lattice model so far. It seems to be problematic because in these cases the instanton action S inst = A is a complicated irrational number but the coefficients of the fluctuation are simple rational numbers. Nevertheless we have observed in the previous section that the threesome relation (3.24) is very likely true even for the honeycomb lattice.

Concluding remarks
In this paper, we analyzed the Bloch electrons in the honeycomb lattice under the magnetic flux. This system is physically important since graphene has the honeycomb structure. The system has a singular behavior at the Dirac points. We presented a systematic way to push the computation of the weak magnetic flux expansion perturbatively. We also looked into the nonperturbative bandwidth of the spectrum, and found a similarity to the supersymmetric sine-Gordon system. This fact suggests that physics near the Dirac points in graphene has a hidden supersymmetric structure (see also [21]). 8 The nonperturbative effects there is probably related to "Cheshire Cat Resurgence" [13,37].
Based on the detailed numerical analysis of the nonperturbative corrections, we finally found a simple threesome relation among the perturbative, one-instanton and instantonanti-instanton sectors. For simple quantum mechanical systems, this should be understood as a consequence of the ZJJ exact quantization conditions [32,33]. The remarkable point is that this relation seems to be true even for the 2d Bloch electron systems, in which neither the ZJJ quantization conditions nor the PNP relation have been understood. These observations suggest that the instanton fluctuation and the bion fluctuation are universally written as P inst fluc (n, g) = ∂P vac fluc (n, g) ∂n e − A(n,g) , P bion fluc (n, g) = ∂P vac fluc (n, g) ∂n e −2 A(n,g) , where A(n, g) is a function that is essentially same as the function A(E(n), g) appearing in the ZJJ quantization conditions. It is not necessary to impose the PNP relation between P vac fluc (n, g) and A(n, g) to derive (1.5). Our relation might give a clue to show up the nonperturbative structure in the Bloch electron systems. We gave a piece of evidence of the relation (1.5) for the honeycomb lattice.
We propose several issues. First, it is desirable to understand (1.5) more deeply in Bloch electron systems. We will investigate it for the Hofstadter model in more detail [18]. Second, it would be interesting to ask what quantum geometry of toric Calabi-Yau corresponds to this honeycomb system. We will report this issue in [7]. In this direction, it is also interesting to unifying the resurgent analysis with the topological string analysis, along the line in [38][39][40][41][42]. Finally, in our analysis, the similarity of the honeycomb lattice and the SUSY sine-Gordon QM was found. As in [11], if introducing an anisotropic parameter in the honeycomb lattice, it leads to a SUSY breaking. Does it has a relation to the deformation in [13,26]? An interesting point is that for particular choices of the parameter in [13,26], the system reduces to a quasi-exact solvable system. It would be interesting to look for such deformed electron systems.   The next-to-leading coefficient is also obtained by the sequence The behaviors of x (2) Q and its fifth Richardson transform are shown in figure 8 (right). By the similar way above, we obtain the numerical value R 35 [x (2) 215 ] = −0.0370370370370370370370370370373283 · · · . (A.8) It is easy to find x (2) ∞ = −1/27. In the last step to guess analytic values of the coefficients, the function Rationalize in Mathematica is useful. Alternatively, there is a very helpful website [44] to find complicated irrational numbers. Repeating this method, we found many analytic values, as in the main text.