{'fraction_non_alphanumeric': 0.03569941715237302, 'fraction_numerical': 0.009020260893699695, 'mean_word_length': 4.453083617101778, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Absorption and optical selection rules of tunable excitons in biased bilayer graphene

Biased bilayer graphene, with its easily tunable band gap, presents itself as the ideal system to explore the excitonic effect in graphene based systems. In this paper we study the excitonic optical response of such a system by combining a tight binding model with the solution of the Bethe-Salpeter equation, the latter being solved in a semi-analytical manner, requiring a single numerical quadrature, thus allowing for a transparent calculation. With our approach we start by analytically obtaining the optical selection rules, followed by the computation of the absorption spectrum for the case of a biased bilayer encapsulated in hexagonal boron nitride, a system which has been the subject of a recent experimental study. An excellent agreement is seen when we compare our theoretical prediction with the experimental data.


I. INTRODUCTION
Since graphene was first isolated [1] many other two dimensional materials have be discovered. Examples of these are the semiconductor monolayer transition metal dichalcogenides (TMDs) [2], the highly anisotropic phosphorene [3], or the insulator hexagonal boron nitride (hBN) [4]. Despite their own unique properties, these three materials have something in common, all of them have an optical response dominated by excitonic effects.
In the simplest possible picture, an exciton is formed upon the excitation of an electron to the conduction band, leaving a hole in the valence band. These two particles, which have opposite charges, interact via an electrostatic potential [5]. This situation is somewhat similar to what is found in the Hydrogen atom, and leads to the formation of bound states inside the band gap of the material. Due to the reduced screening in the out-of-plane direction, two dimensional materials host excitons with large binding energies , easily surpassing the energies usually found in their three dimensional counterparts. Furthermore, excitons in these 2D systems couple efficiently with light resulting in huge oscillator strengths. This has made the exploration of excitonic effects a rich research field, with many possible applications, such as deep-UV optoelectronics [4,6,7] in hBN, or the exploration of valleytronics [8] and single photon emitters [9,10] in TMDs, which are of great importance in the field of quantum information. Furthermore, it has been recently shown that the combination of a TMD with a Van der Waals heterostructure cavity enhances the light-mater interaction, thus allowing for unitary excitonic absorption [11].
Although present in many 2D materials, low energy excitonic bound states are noticeably absent in pristine monolayer graphene, something easily understood if one recalls that this material is a semi-metal, thus lacking the necessary band gap for the formation of such states. This, however, does not stop us from studying excitonic physics in graphene based systems, since through material engineering many systems can be designed to obtain the desired properties. A clear example of this is the case of biased bilayer graphene (BBLG), where two graphene monolayers are stacked and a displacement field is applied to them. Even though bilayer graphene (BLG) inherits the lack of a band gap from its monolayer form, the presence of the displacement field opens a gap in the system, which can be finely tuned from the mid to farinfrared [12][13][14]. This finite gap is the key ingredient to combine the desirable properties of excitons with the unique physics of graphene based systems.
The first theoretical study about excitons in BBLG was carried out in Ref. [15], where a pronounced resonance, sensitive to the displacement field, was predicted to appear inside the band gap in absorption measurements. More recently, in Ref. [16], an experimental study was carried out to characterize the excitonic optical response of BBLG encapsulated in hBN. Our goal with this paper is to produce a theoretical description of the experimental data measured in that work. To achieve this, we will describe the electronic properties using a tight binding model, while solving the Bethe-Salpeter equation (BSE) to obtain the energies and wave functions of the excitons. While oftentimes the BSE is coupled to DFT calculations , and solved in a fully numerical manner [15,[17][18][19][20], requiring huge computational effort, our approach uses a simple analytical treatment which requires a single numerical quadrature, thus greatly reducing the numerical cost of the computation. This paper is organized as follows. In Sec. II we introduce the tight binding model which describes the biased bilayer, and diagonalize it. In Sec. III, we start the discussion on the excitonic properties of the system, which we obtain by solving the Bethe-Salpeter equation. In Sec. IV, we discuss the optical properties of the biased bilayer. First, we discuss the optical selection rules; then we compute the absorption spectrum and compare it with the experimental results of Ref. [16]. Finally, in Sec. V, we give an overview of the paper and our closing remarks. An appendix describing how to efficiently solve the Bethe-Salpeter equation closes the paper.

II. TIGHT BINDING MODEL
To start our discussion let us study the electronic properties of a graphene bilayer subject to an external bias, which we schematically depict in Figure 1. To character- ize this system we shall employ a tight binding Hamiltonian written directly in momentum space. To obtain this Hamiltonian we account for nearest neighbors intralayer hoppings (γ 0 ), as well as nearest neighbors (γ 1 ) and next nearest neighbors (γ 3 , γ 4 and γ 5 ) interlayer hoppings. Doing so, one finds the Hamiltonian to be: when written in the basis {|1, b , |2, b , |1, t , |2, t } (with b and t denoting the bottom and top layers, respectively, and the labels 1 and 2 are defined as in Fig. 1). Here, V quantifies the bias and φ(k) = e iaky + e ia(kx √ 3−ky)/2 + e −ia(kx √ 3+ky)/2 is a phase factor determined by the honeycomb geometry of the individual layers. Since we will be mostly interested in the low energy response of the system, we shall restrict our study to the vicinity of the Dirac points of the the first Brillouin zone. Mathematically, this amounts to approximating the phase factor φ(k) as φ(k) ≈ 3 2 a(τ k x + ik y ), with τ = ±1 labeling the two Dirac points. Note that henceforth, the momentum k = (k x , k y ) is measured relatively to the chosen Dirac point. With this approximation one finds the low energy Hamiltonian where we introduced the Fermi velocity defined as v F = 3aγ 0 /2, and θ = arctan k y /k x . For simplicity, in what follows, we will be mainly concerned with the case where γ 3 = γ 4 = γ 5 = 0, i.e. accounting only for nearest neighbors hoppings. However, later in the text, we shall see how these hopping parameters affect the response of the system.
Diagonalizing the low energy Hamiltonian of Eq. (2), we find the electronic bands depicted in Fig. 2, which were obtained using γ 0 = 3 eV and γ 1 = 0.4 eV. In this On the left panel we show the case without bias, and on the right one the bands when a bias of V = 60 meV is considered. For convenience we label the valence and conduction bands with an additional index η, which is equal to −1 for the bands closer to the middle of the gap and 1 for the other two.
figure, we see that, for zero bias, one finds the usual bands of bilayer graphene, which host no gap since the lowest conduction band and the highest valence band touch at k = 0. When the bias is increased to V = 60 meV, the two aforementioned bands are deformed, taking the shape of a "Mexican hat". This deformation of the valence and conduction bands leads to the opening of a gap, something crucial to exploit excitonic effects. If the next nearest neighbors hoppings had been considered only minor changes would appear in the band structure. Even though small, these modifications in the bands have an important effect on the optical selection rules (which we discuss ahead), since they introduce the trigonal warping effect, which reduces the symmetry of our calculation from C ∞ to C 3 .
From the diagonalization of Eq. (2), we also find the eigenvectors associated with each band. We label these vectors as |c, η, k and |v, η, k , where the index η is defined as in Fig. 2, i.e. the bands closer to the middle of the gap have η = −1, while the other two have η = 1; they read: The specific definition of each entry is not of particular interest to the current analysis, however, the phase of each spinor was chosen carefully. As a consequence of the definition of θ, as k → 0 the complex exponential e iτ θ is discontinuous. To avoid this discontinuity in the eigenvectors, the phase of each one was chosen such that the complex exponentials appear multiplied by terms that vanish in the limit of small k. It is by now known that the phase choice of the Bloch factors plays a crucial role in determining the optical selection rules, and the criteria we established here is the one that directly leads to Hydrogen-like selection rules in the monolayer. Although other phase choices can be employed, they require the introduction of an additional angular quantum number associated with the pseudospin texture of the excitonic states [15,21,22]. Even though we do not shown the explicit form of the entries of each eigenvector, in light of the discussion regarding where the complex exponentials are place, it should be clear These relations between the amplitude of the different vector components are enough for us to gain intuition about the coupling strength of light with different excitonic states, which we introduce next.

III. BETHE-SALPETER EQUATION
Now that the single particle regime was studied, we move on to the excitonic part of the problem. To obtain the response of the system due to excitons we must first determine their energies and wave functions. These quantities can be obtained from the solution of the Bethe-Salpeter equation. This integral equation in momentum space may be written for a multi-band system as [23,24]: where, to simplify the notation, we omitted the index η (which is now included in the indexes c and v). We note that this equation takes the single particle energies E c/v,k and eigenvectors |c/v, k as the input, and couples the different bands through an electrostatic potential V (k), thus capturing the many body nature intrinsic to excitons. By solving this equation, one finds the exciton energies, E, and the associated wave functions, ψ cv (k). In our approach we consider the electrostatic potential V (k − q) to be the Rytova-Keldysh potential [25,26] (which is usually employed to describe excitonic effects in monolayers [5]). This potential can be obtained from the solution of the Poisson equation for a charge embedded in a thin film, by taking the limit of vanishing thickness. Although this potential presents a complex form in real space, its representation in momentum space is rather simple where α ∼ 1/137, is the mean dielectric constant of the media above and below the bilayer and r 0 corresponds to an in plane screening length which is related to the 2D polarizability of the material; its value can be found from the single particle bands and Bloch factors. The numerical value of r 0 is of great importance to accurately describe the excitonic properties of a given system, since it directly affects the screening experienced by the carriers. Although the values of r 0 are well documented for TMDs, hBN and other 2D materials, the same is not true for biased bilayer graphene, where, in principle, r 0 varies with the bias V . In Ref. [27] the value of r 0 as a function of V was determined, however, it was done using a 4 band tight binding model, which may not be sufficiently accurate. Ideally the value of r 0 should be determined from ab initio calculations taking into account several valence and conduction bands, because, as discussed in Ref. [28], depending on the material, a small number of bands of a minimal tight binding model may or may not be enough to accurately determine r 0 . Since a detailed ab initio study of r 0 in BBLG is lacking in the literature, we will use the results of [27] as a reference. There one finds that, as expected, r 0 decreases with increasing band gap, or alternative, it decreases with increasing bias in a non-trivial way.
To solve the BSE we start by noting that Eq. (7) corresponds in fact to four equations, one for each possible pair (v, c), forming an eigenvalue problem. To more easily solve the BSE, we assume that the excitons have a well defined angular quantum number m, such that ψ cv (k) = f cv (k)e imθ . Using this, and manipulating the phase choices of the spinors in Eq. (3), one can show that the BSE can be reduced to a 1D integral equation, which can then be efficiently solved with a single numerical quadrature (the details on how to achieve this are discussed in the Appendix).
After solving the BSE for a wide range of values for the bias (whilst using the results and parameters of the previous section, as well as the r 0 of Ref. [27] and taking = 6.9 [29], corresponding to the case of hBN encapsulated BLG), we observed that of the four possible sets of ψ cv (k), the contribution of the one with η c = η v = −1 was by far the dominant one, so much so, that the other wave functions appeared to vanish in comparison. This is a reasonable result to find, since intuition tells us that the bands with η = −1, i.e. the ones closer to the middle of the gap, should dominate the low energy response. Hence, to further optimize the calculation, one may restrict the BSE to a single pair of bands, thus greatly reducing the computational cost of the problem, without hindering the quality of the results.

IV. OPTICAL ABSORPTION
In this section our goal is to determine the optical selection rules of the biased bilayer and to compute its optical absorption. To achieve both these goals we shall start by looking at the optical conductivity of the system. In the dipole approximation, and considering normal incidence, the optical conductivity follows as [24] where the sum over n accounts for the different excitonic contributions with energy E n , and the second term corresponds to the non-resonant part of the optical response. The quantity Ω n is defined as where v, k|r|c, k is the interband matrix element of the dipole operator. Also, we introduced a phenomenological linewidth, Γ n , which we assume to be n−dependent, to more easily compare our theoretical predictions with experimental data. As we have discussed before, we consider only the contribution of the bands with η = −1, since the ψ cv associated with the other ones are essentially zero. To evaluate the dipole matrix element we use the relation Although when solving the BSE we considered only the contribution of the nearest neighbor hoppings γ 0 and γ 1 , to understand the role played by the other hoppings in the optical response, we shall take them into account when the commutator [H, r] is evaluated.
Let us now consider that our system is excited by a circularly polarized electric field. The relevant dipole matrix elements for such a situation are v, k|x ± iy|c, k , corresponding to positive and negative circular polarization. Using Eq. (11) to evaluate the matrix element, we find that the terms proportional to γ 0 give the contribution It is then clear from Eq. (10), that such a result implies that only states with m = −τ and m = −3τ , i.e. p− and f −states, can be excited. From this result, one also sees that positive polarization excites states with m = 1 in the valley with τ = −1, and states with m = −3 in the opposite valley; an analogous results holds negative polarization. Regarding the coupling strength between the photons and the allowed exciton states, and recalling the comments made on the magnitude of eigenvector entries after Eq. (3), one can expect the resonances associated with the p−states to be far more intense than the ones associated with the f −states, which should present a vanishingly small oscillator strength.
The contributions to the matrix element v, k|x ± iy|c, k proportional to γ 3 and γ 4 , lead to the same selection rules as the ones we have just discussed. Furthermore, since γ 1 γ 3 , γ 4 , one finds that these hoppings can be safely ignored when studying the linear optical response of the bilayer. In stark contrast with this observation, the part of the dipole matrix element which is proportional to γ 5 introduces new selection rules. Explicitly: where we see that by accounting for the contribution of  switches the roles of the two circular polarizations. Moreover, since linear polarization can be viewed as the linear combination of the two circular polarization components, the analysis we performed is easily applied to the case of linearly polarized electric fields.
Having analyzed the optical selection rules of the biased bilayer, and making use of the solutions of the BSE, the optical conductivity is readily obtained. A quantity which is determined by the optical conductivity and is experimentally more relevant is the optical absorption. Here, we define the absorption as A = 1 − |r| 2 − |t| 2 , where r and t are the reflection and transmission coefficients, respectively. These two coefficients are found from the longitudinal conductivity as [30]: with α ∼ 1/137, σ( ω) the optical conductivity and σ 0 = e 2 /4 the conductivity of a graphene monolayer.
Considering an in-plane linearly polarized electric field (for example along the x direction), using once again = 6.9 (which describes a system encapsulated in hBN), and considering γ 5 = 0.1γ 1 [27], we find the absorption spectra depicted in Fig. 4 (a). Analyzing these results, we see that for each bias two resonances are well resolved, the smaller one at lower energies being due to the 1s exciton, and the larger one corresponding to the excitation of the 2p states. The linewidths of these resonances were chosen in order to match the experimentally measured values of Ref . [16], which read Γ 1s = 0.4meV and Γ 2p = 1.3meV. Although optically allowed, the f − and g−states do not originate any noticeable resonance due to their small oscillator strength; the same is true for the more excited s−states. The more energetic p−states present a significant oscillator strength, however, since these are too close together, and too close to the band gap, they can hardly be resolved; instead, these states form the beginning of a plateau. This structure with an approximately constant magnitude is then extended beyond the band gap, due to the contribution of the interband excitations. Superimposed on the excitonic response, we also depict the interband conductivity obtained in the single particle (SP) picture using Fermi's golden rule. This allows us to see that the excitonic effect is not only responsible for the resonances inside the gap, but also leads to a different response at higher energies, as a result of the shift in oscillator strength from the interband response to the excitonic resonances. The results of Fig. 4 (a) agree well with those of [15], which were obtained through more complex numerical approaches. Those results, however, do not account for the s-resonance which we include in our calculation.
Focusing now on the dependence of the absorption with the the bias, we observe that as the bias is increased the main features of the absorption remain qualitatively the same, but the optical response is shifted to higher energies. Further evidence of this behavior is given in panel (b) of the same figure where we depict the position of the 1s and 2s resonances as a function of the bias. This allows us to confirm that the resonances are in fact shifted to higher energies as the bias increases, while also becoming more separated from each other. The shift of the resonances to higher energies with increasing bias was to be expected, since, as we also depict, the band gap grows as the bias is turned up. However, we must note that the trend of the peak's position differs slightly from that of the band gap. This is effect is due to dependence of r 0 with V , since as the bias increases, r 0 decreases, leading to more tightly bound excitons, appearing deeper into the gap. These two competing effects result in the almost linear dependence of the location of the resonances with V .
Finally, in panel (c) of Fig. 4, we compare our theoretical prediction with the experimental points of Long Ju and coworkers [16]. To facilitate the comparison between theory and experiment, we chose a bias of V = 52 meV since it produced the best description of the experimental results; the remaining parameters retain the same values we considered so far (with r 0 taken from [27]). Comparing the theoretical prediction with the experimental data an excellent agreement is seen, in both the position and the magnitude of the two resonances. Besides this, we also see that the theoretical plateau at higher energies matches well with the experimental one. Although our plateau appears to be further away from the 2p resonance than in the experimental data, it seems possible to assign the structures that are experimentally captured on the high energy side of the largest peak to other excitonic state, which, due to their proximity, merge together, creating a single broader resonance.
At last, in Fig. 5, we depict a density plot of the 1s and 2p exciton wave functions for different bias. As expected, for every case the 1s wave function is finite at k = 0 while it vanishes for the 2p states. Moreover, as the bias increases the wave functions become more delocalized in momentum space. This increased delocalization in momentum space, translates to a higher localization in real space, in agreement with the aforementioned idea that the excitons become more tightly bound as the bias rises. This relation between real and momentum space also allows us to estimate the size of an exciton. Taking the case of V = 70 meV as an example, we see that its 2p wave function presents a spread in momentum space of about ∆k = 0.2 nm −1 . Using ∆r = 2π/∆k, we estimate that the same exciton should be extended over approximately 40 nm in real space.

V. CONCLUSION
In this paper we studied the excitonic response of biased bilayer graphene. To achieve this, we started by describing the single particle electronic properties of the system using a tight binding Hamiltonian. Afterwards, the single particle solutions were used as the input to the Bethe-Salpeter equation, whose solution defines the excitonic properties, i.e. their energies and wave functions. To solve the BSE we transformed the 2D integral equation into a 1D problem by carefully choosing the phase of the Bloch factors. This reduction in the dimension of the problem allows us to more efficiently solve the BSE, which now requires a single numerical quadrature.
To study the optical properties of the biased bilayer, we computed its conductivity, which allowed us to extract the optical selection rules and gain intuition on the strength of the photon-exciton coupling. We found that, in agreement with the literature, when trigonal warping is ignored, only p−states (with m = ±1) are bright; even though f −states are optically allowed they present a tiny oscillator strength. When trigonal warping is considered, and thus the symmetry of the calculation is reduced from C ∞ to C 3 , we find that s−states (with m = 0) become optically allowed, and present a noticeable oscillator strength (at least for the 1s state).
Making use of the optical conductivity, we computed the absorption spectrum for the case where a biased bilayer graphene is encapsulated in hBN, which corresponds to the system where an experimental study was carried out by Long Ju and coworkers in [16]. Comparing our theoretical prediction with the experimental data, an excellent agreement was observed. Our model accurately described the position and magnitude of the 1s and 2p resonances, while also capturing the transition from excitonic resonances to the continuum of states, which originates the plateau that appears in the experimentally measured absorption. The dependence of the resonance's location with the bias also agreed with the experimental data.
Note added: After the completion of this work we became aware of a similar unpublished work [31] on the same system we treated in this paper.
Appendix A: On the solution of the Bethe-Salpeter equation In this appendix we shall give a more in depth description on how to numerically solve the Bethe-Salpeter equation (BSE). This approach is an extension of the method presented in Ref. Ref. [32] to treat excitons in 2D semiconductors using the Coulomb potential. We take Eq. (7) of the main text as our starting point. Using the ansatz ψ cv (k) = f cv (k)e imθ , and taking the thermodynamic limit, we find: This equation corresponds to a two dimensional integral equation whose numerical solution is computationally demanding. To reduce the numerical weight of the calculation we wish to reduce the problem to that of a one dimensional integral equation. This can be achieved if the product of Bloch factors has the form since this allows us to perform a variable change from to ϑ = θ q −θ k , thus eliminating one of the momentum dependent integrals; here λ is some integer, and A cc vv λ (k, q) are coefficients determined by the explicit computation of the spinor product. In principle this is always possible for the materials we are considering, given that the phases of the Bloch factors are appropriately chosen. Inserting this into the previous equation, and noting that V (k−q) ≡ V (k, q, θ q −θ k ), one finds where we explicitly introduced the variable change dθ q → dϑ. Now, recalling the definition of Rytova-Keldysh potential V (k, q, ϑ), we introduce a new function, I ν (k, q) corresponding to the integral over dϑ: with κ(k, q, ϑ) = k 2 + q 2 − 2kq cos ϑ. Notice how only cos(νθ) enters the integral, since the analogous term in sin(νϑ) vanishes due to parity. From inspection, it should be clear that when q = k the function I ν (k, q) is numerically illbehaved, and as such must be treated carefully. For convenience we shall express I ν (k, q) in terms of partial fractions as where from these two terms only the first one, J ν (k, q), is problematic when k = q. Before we explain how to avoid this numerical problem, let us first write the BSE in a more useful way. First, we write (A7) Then, we define B cc vv m (k, q) = λ J m+λ (k, q)A cc vv λ (k, q) and C cc vv m (k, q) = λ K m+λ (k, q)A cc vv λ (k, q). Using these, one finds Let us now focus on the numerical problem associated with B cc vv m (k, q). To treat the divergence that appears when k = q, we introduce an auxiliary function g m (k, q) and introduce the modification with g m defined such that lim q→k B cc vv m (k, q) − g m (k, q) = 0. Following Ref. [32], we define g m as g m = B cc vv m (k, q) 2k 2 k 2 + q 2 (A10) With the analytical part of the calculation taken care of, we shall now discuss how to numerically solve the equation we have arrived to. First, we introduce a variable change which transforms the improper integral over [0, ∞), into one with finite integration limits, such as [0, 1]. One possibility is to define q = tan[πx/2]. Afterwards, we discretize the variables k and x (and consequently q): where N is the number of points and w j is the weight function of the chosen numerical quadrature; also, q j ≡ q(x j ) and dq/dx j ≡ [dq/dx] x=xj . Furthermore, we note that ∞ 0 g m (k i , p)pdp is numerically well behaved as opposed to the original integral, ∞ 0 B cc vv m (k i , p)pdp. Depending on the problem, and on how localized the excitonic wave functions are in momentum space, it may be useful to split the original improper integral into regions, in order to allow for a thinner mesh in the relevant portion of k-space, and a coarser one in the region where the wave functions are already close to zero.
Regarding the choice of quadrature, we employ a Gauss-Legendre quadrature, which is defined as [33] b a f (y)dy ≈ where with ξ i the i−th zero of the Legendre polynomial P N (y), and with P N (ξ i ) ≡ [dP N (y)/dy] y=ξi . At last, the only thing left to do is to realize that this equation can be expressed as an eigenvalue problem of a 4N × 4N matrix. This matrix can be thought of as a 4 × 4 matrix of matrices, each one with dimensions N × N . The 16 blocks come from the different combinations of the indexes c, c , v and v , with each block corresponding to a N × N matrix stemming from the numerical discretization of the integral. Solving the eigenvalue problem one finds the exciton energies and wave functions.