{'fraction_non_alphanumeric': 0.03348900331297649, 'fraction_numerical': 0.01069732167804255, 'mean_word_length': 4.5453740593182825, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coulomb effects in the absorbance spectra of two-dimensional Dirac materials

A wide range of materials like graphene, topological insulators and transition metal dichalcogenides (TMDs) share an interesting property: the low energy excitations behave as Dirac particles. This emergent behavior of Dirac quasiparticles defines a large class of media that are usually called Dirac materials. The linear and nonlinear optical properties of Dirac materials with a gap are still largely unexplored, and in this Letter we build the foundations of a novel way to study the linear optical properties of these two-dimensional media. Our approach is based on a new Dirac-like formulation of the standard semiconductor Bloch equations used in semiconductor physics. We provide an explicit expression of the linear absorbance -- which we call the relativistic Elliott formula -- and use this to quantify the variation of the continuum absorbance spectrum with the strength of the Coulomb interaction (the Sommerfeld factor). Our calculations also show how the Coulomb enhancements scales with the bandgap and vanishes for zero bandgap, shedding new light on the behaviour of graphene for low light intensities. The results presented are in good quantitative agreement with published experimental results. Our new theory will allow researchers to explore the nonlinear interactions of intense, ultrashort pulses with TMDs, and the framework is flexible enough to be adapted to different experimental situations, such as cavities, multilayers, heterostructures and microresonators.

Condensed matter physics is witnessing a rapid expansion in the fabrication of a wide variety of materials with Dirac fermion quasi-particle excitations. These seemingly diverse materials possess properties that are direct consequence of the Dirac spectrum of the quasiparticles and are universal. For example neutral superfluids d-wave superconductors and graphene, which are characterised by a massless Dirac fermion low energy spectrum, all exhibit the same power-law temperature dependence of the fermionic specific heat with the only differences arising from the dimensionality of the excitation phase space [1]. Amongst Dirac materials, gapless graphene is one of the most widely studied due to its unusual physical properties arising from the interplay of its reduced dimensionality and the nature of its excitation spectrum [2][3][4][5].
Perhaps the most intriguing characteristic of graphene is the so-called universal absorbance. The opacity of suspended graphene is defined solely by the fine structure constant (α = e 2 / c), the parameter that describes coupling between light and relativistic electrons which is traditionally associated with quantum electrodynamics (QED) rather than materials science [4]. This universal behavior is known to be broken by opening a gap, with the appearance of bound states (excitons) due to the electron-hole Coulomb interactions. A gap can be opened in graphene when the sample is deposited on a dielectric substrate [6] or it can be induced by impurities, lattice defects [7] and mechanical strain [8]. In terms of the low excitation Dirac Hamiltonian the effect of a gap opening is analogous to the chiral symmetry breaking in QED that generates a mass gap in the particle-antiparticle spectrum [9,10]. This suggests that a gapped Dirac material can be accurately described by a massive Dirac Hamiltonian. An example of two-dimensional ma-terials whose low energy excitations behave as massive Dirac particles is provided by monolayer transition metal dichalcogenides (TMDs) such as molybdenum disulfide (MoS 2 ) and tungsten diselenide (WSe 2 ) [11]. Akin to graphene, these materials display new physical properties, distinct from their bulk counterparts [12][13][14][15][16]. Here, emission is dominated by excitons and trions due to the strong Coulomb interactions arising from their low dimensionality and reduced dielectric screening. Remarkably, despite their tiny atomic thickness they can absorb up to 20% of the incident light [15,16]. The proper analysis of the measured spectra requires the identification of the optically active states. For this purpose, one derives dipole allowed selection rules that result from the conservation of the total angular momentum in the excitation and emission processes. In conventional direct gap GaAstype quantum wells the radial nature of the dipole moment imposes that only s-type excitonic states couple to light [11,17]. This leads to the well-known excitonic Rydberg series. It has been shown that in the case of monolayer TMDs in the regime of strong Coulomb interactions the system collapses into a excitonic insulator phase and optically bright p-excitons are allowed [11]. However recent works pointed out that for a gapped Dirac system the symmetry properties of the K-point result in a non radial dipole moment with a nontrivial angular dependence related to the appearance of a Berry phase [18]. This peculiar characteristic suggests that different dipole allowed transitions should be present and optically active p-states should appear even in the weak Coulomb interaction regime.
In this Letter we calculate the absorbance spectra of Dirac materials in the relatively new theoretical framework of the instantaneous eigenstates [19] and we apply this formalism to predict the spectra of realistic mate-rials. We derive the Elliott formula for a gapped Dirac system taking into account both the discrete and the continuum part of the Wannier eigenvalue problem. In particular we study how the absorbance spectrum evolves as the bandgap is reduced to zero. We show that even in the presence of relatively strong electron-hole Coulomb interaction the gap reduction leads to the universal absorbance of graphene and we show how this behavior is closely related to the structure of the relativistic hydrogenic problem.
Before detailing our calculations we show in figure 1 the absorbance spectrum calculated for the Molibdenum disulfide (MoS 2 ). We observe a high absorbance value (15%) around the energy gap in agreement with previous experimental and theoretical results [16,20,21]. Further agreement with experimental data can be observed in the strength of the 1s exciton peak (5%) and the value of the binding energy [15]. The intensity of the continuum absorbance peak is around three time stronger than the 1s exciton one. This is a peculiar characteristic shared by Dirac materials and in contrast to non-Dirac materials.
In what follows we show how this spectrum relates to the relativistic structure of the low-energy quasiparticle excitations. The starting point for our analysis is the low-energy Hamiltonian for the band structure in the vicinity of the Dirac points (k = 0). Governed by the symmetry properties of the hexagonal lattice, the lowest order k · p Hamiltonian including the light-matter interaction and neglecting the spin-orbit interaction has the form where v F is the Fermi velocity, A(t) is the vector potential, ∆ is the energy gap, I 2 is the 2 × 2 identity matrix, V (r) is the electrostatic potential and σ = (σ x , σ y ), σ z are the Pauli matrices. The usual approach used to study the optical properties of this Hamiltonian parallels that of the Semiconductor Bloch equations (SBEs) [22]. However recent works show that a full understanding of the nonlinear optics of Dirac materials requires one to go beyond the SBEs [18,19,23]. In 2010 Ishikawa derived an extended version of the SBEs for graphene using the formalism of instantaneous eigenstates [19] which has been recently generalised for a gapped material [18]. The theory developed in these works unveil novel and previously unexplained nonlinear optical properties of graphene and gapped graphene, however they do not include Coulomb interactions that are in general very strong in two-dimensional semiconductors, preventing meaningful predictions to be made. The inclusion of Coulomb interactions within this new formalism requires the second quantisation of Hamiltonian Eq. (15). To accomplish this, we need to quantise the instantaneous Dirac field; the procedure is simple but rather lengthy (see Supplementary Material) and it results in the second quantised Hamiltonian (2) where a k and b k are respectively the annihilation operators for electrons and holes, while ] are the instantaneous energy and Rabi frequency respectively. This Hamiltonian differs from the usual second quantised 2D semiconductor case because both k (t) and Ω k (t) explicitly depend on the interaction field A(t). The presence of the Coulomb potential V (r), leads to a renormalisation of both the instantaneous energy and the Rabi frequency. Introducing the microscopic polarisation p k = b −kâk and the electron-hole occupation numbers n e k = â † kâ k , n h k = b † −kb −k we can show that they satisfy the following system of renormalised Dirac-Bloch equations (RDBs): are the renormalised energy and Rabi frequency. The phase term (γ k ) is a gap induced Berry phase (see Supplementary information for details). The formulation of Eqns 2 and 27 are the first results of this Letter. As was the case for the SBEs [22,[24][25][26], the RDB equations framework can be flexibly applied to many optical pulsed experiments and provide a route to the microscopic understanding of linear and nonlinear properties of gapped TMDs. This approach will pave the way for a complete understanding of the excitonic structure and the interaction of ultrashort, intense light pulses with Dirac materials, an area that at present is completely unexplored and potentially full of groundbreaking discoveries.
As a first example of the practical application of our new equations, we now consider the RDBs in the low intensity and low density limit (e/cA(t) << ∆/(2v F ) and n e k − n h k = −1) where only the polarization equation survives. We will also assume a continuous wave (CW) radiation field. This will allow us to explore analytically excitonic effects on the linear optical properties of gapped TMDs. Eq. (27) then reduces to with the non-renormalised energy reducing to k = v 2 F p 2 + ∆ 2 /4. To this equation is associated the well known Wannier stationary eigenvalue problem, that, for a massive Dirac quasiparticle reads here Ψ ν (r) (ν = n, j) is a two component spinor eigenfunction of the electron-hole pair states which in polar coordinates has the form where j = m + 1/2 is the eigenvalue of the "isospinangular" momentumĴ z =L z + 1 2 σ z along the z axis and n is the principal quantum number. The corresponding eigenvalues are with γ = j 2 − α 2 c . The constant α c is the dimensionless Coulomb coupling strength and is determined by the background dielectric constant [27].
In what follows we use the Wannier states in Eq. (6) to derive the electric susceptibility and hence the linear absorbance; thereby recasting the Elliott formula for gapped Dirac materials. To solve Eq. (4) we expand the microscopic polarisation in term of these Wannier states, after a lengthy but straightforward calculation [17] we obtain, Here µ cv (r) = µ cv (r) 1 1 is the electric dipole in spinor form, E(ω) is the electric field in frequency domain, L 2 the area of the sample and δ is an energy broadening. We can now define the macroscopic polarization as follows and inserting Eq. (8) we get where the oscillator strength is given by (11) From the definition of polarisation in c.g.s units P (ω) = L 2 d χ(ω)E(ω), we finally obtain the electron-hole pair susceptibility as The nominal layer thickness d has to be introduced to compute the volume within which the dipoles are induced but as we shall see it falls out in the calculation of the absorbance A which is the absorption coefficient multiplied by the layer thickness d. The Influence of this effective layer thickness on the excitonic ground-state has been studied in a recent work [28]. From Eq. (12), using the Plemelj-Sokhosky theorem, we obtain the final expression of the linear absorbance: the relativistic Elliott formula for Dirac materials where Ω = ω/∆ is the scaled photon energy and Θ is the Heaviside step function. The two terms of Eq. (13) are the discrete and the continuum contributions of the absorbance spectrum. Despite its striking similarity to the non-relativistic case Eq. (13) is not trivial to evaluate since the integral given in Eq. (11) can not be solved in closed form and must be evaluated numerically with particular care (some detail about the structure of the dipole moment are given in the Supplementary material). The factor C(Ω) in the second term of Eq. (13) is what we may call the relativistic Sommerfeld factor: Where the sum here extend over the continuum states. This quantifies how the presence of electron-hole interactions enhances the continuum absorbance with respect to the free carrier limit [17]. It is maximum at the energy gap (Ω = 1) with a value that depends on the coupling constant α c and tends to unity as the photon energy increases. This behavior is closely related to the relativistic energy spectrum and to the structure of the electron-hole continuum wave-function in contrast with the non-relativistic case which has a value of 2 at the bandedge. We now use Eq. (13) to demonstrate the development of the Coulomb effects on the absorbance spectrum for varying dielectric screenings and energy gaps. In figure  2 we plot the absorbance spectrum while varying the strength of the Coulomb interaction α c . Experimentally it is possible to tune the strength of the Coulomb interaction by changing the substrate on which the sample is deposited since it changes the overall dielectric constant. As expected the bound-state exciton lines reduce their binding energy and oscillator strength as the Coulomb interaction weakens (i.e. reducing α c ). For the continuum states, the Coulomb enhancement (Sommerfeld factor) is strongest at the bandedge, see Eq. (14). Approaching the zero Coulomb limit (α c = 0) Eq. (13) reproduces exactly the free carrier absorbance spectrum that can be computed analytically as A(Ω) = πα (Ω 2 +1)/Ω 2 . As can be seen from the dashed-dotted line in figure 2 this noninteracting model is unable to reproduce the giant continuum absorbance seen in experiment. It is therefore clear that the giant absorption of Dirac media such as TMDs is due to Coulomb interactions. The absorbance reduces quite rapidly as the photon energy is increased and at high photon energy the absorbance approaches the universal (graphene) absorbance value of 2.3%. In this limit the high energy electron and hole scattering wavefunctions involved are essentially identical to the free particle wavefunctions. The inset in figure 2 shows how the absorbance at the band-gap varies with α c . We can clearly observe that it increases dramatically as the coupling constant approaches the critical value (α c = 0.5) where a phase transition to an excitonic insulator occurs [11,29,30]. This phase transition is a fundamental signature of purely Dirac excitons and it potentially represents a measure of the Diracness [31] of the two body electron-hole system since it disappears for a parabolic exciton dispersion [32]. One of the most remarkable features of graphene from a physics point of view is that there appears to be no Coulomb enhancement of the continuum absorbance, despite the existence of final state Coulomb interactions. Figure 3 shows the continuum absorbance calculated from Eq. (13) varying the energy gap (∆). We observe that reducing the bandgap the curves approach the universal absorbance value of graphene πα at lower and lower photon energies. This indicates that the relevance of Coulomb interactions is confined to a region close to the energy gap and once the gap is closed the system behaves as a nearly free electron gas as seen in experiments [4] and in the Das Sarma paper [33] by using a renormal-isation group approach. Within our theoretical approach this phenomenon can be explained by observing that in the relativistic hydrogenic model the energy gap is the exact analog of the Rydberg constant that is the only relevant energy scale. This means that tuning the gap changes the relevant energy scale of the system; when it is reduced to zero all photon mediated transitions can be considered as high energy transitions and in this regime the system it is known to reproduce the free absorbance limit.
In conclusion we studied the linear optical properties of two-dimensional gapped Dirac materials. We derived a new system of Coulomb-renormalised Dirac Bloch equations based on the quantisation of instantaneous eigenstates. We provided an explicit expression of the absorbance Elliott formula that to the best of our knowledge was still missing in previous theoretical studies. We demonstrated how the absorbance spectrum evolves as the energy gap and the strength of Coulomb interaction are reduced to zero. In particular we showed that even in the presence of relatively strong electron-hole Coulomb interactions the gap reduction leads to universal absorbance and we relate this to the structure of the relativistic hydrogenic problem. Our new theoretical framework will be potentially used to predict the behaviour of the interaction of ultrashort pulses with Dirac media in a variety of experimental situations, leading to a fresh understanding on new mechanisms of higher-harmonics generation and the study of cavities, multilayers and TMDenhanced microresonators.
LDMV acknowledges support from EPSRC under the auspices of the Scottish Centre for Doctoral Training in Condensed Matter Physics.

QUANTISATION OF THE INSTANTANEOUS DIRAC FIELD AND RENORMALISED DIRAC-BLOCH EQUATIONS
The starting point of our analysis is the low-energy Hamiltonian for the band structure in the vicinity of the Dirac points. Dictated by the symmetry properties of the hexagonal lattice, the lowest order k · p Hamiltonian including the light-matter interaction has the form where A µ (x, t) = (A(t), V (r)) is the three-vector potential, ξ is a valley index referring to the Dirac points (K, K ) and σ = (ξσ x , σ y ), σ z are the Pauli matrices. Since we are not considering the effect of the spin-orbit coupling we can restrict to a single valley, namely ξ = 1. The goal of this section is to derive a renormalised system of Dirac-Bloch (RDB) equations from the second quantised form of the Hamiltonian in eq. (15). In order to do so we consider the following interacting Dirac equation An ansatz solution to this equation can be given in terms of instantaneous eigenstates of the secular equation These eigenstates can be written in the normalised from here |π k | = p + e c A(t), θ k = arctan(π x /π y ) and λ = ±1 labels the conduction and valence bands. k (t) is the positive branch of the instantaneous energy.
The addition of the gap leads to an inequivalence of K and K sublattices and consequently to the appearance of a Berry phase γ λ,k (t) = t −∞ dτ u † λ,k (τ )· u λ,k (τ ). Consequently the spinor (18) evolves in time as follows Where θ k (t) = 1 t −∞ k (τ )dτ is the dynamical phase and γ k (t) the positive branch Berry phase. The Dirac field ψ(r, t) can be then expanded in terms of instantaneous eigenstates with a k,λ are band electron ladder operators. Substituting this field in the Hamiltonian we get with η(t) = θ k (t) + γ k (t) being the sum of the Berry phase and the dynamical phase. Introducing the hole operators b † −k = a k,−1 , after a lengthy but simple algebraic calculation we get the final form of the second quantised Hamiltonian where Ω k (t) is a generalised Rabi frequency and it is given by Introducing the Coulomb interactions in the Hamiltonian we get we can now derive the equations for the population and inversion variables n e k = â † kâ k , n h k = b † −kb −k and p k = b −k a k . This lengthy procedure is quite standard and can be found in any book of semiconductor theory; for example [1]. It results in the following set of RDB equations where R and Ω R (t) = Ω k (t) + 1/ k =k V |k−k | p k are the renormalised energy and Rabi frequency.

THE WANNIER-DIRAC EQUATION
As in the case of 2D semiconductors the Wannier equation can be derived from the RDB equations (27). We consider the system under the approximation e/cA(t) << g /(2v F ) and under the low excitation approximation w k ≈ −1 Thus Eq. (27) simplifies to with the non-renormalised energy reducing to v F p 2 + ∆ 2 /4. This equation represents a two dimensional relativistic two particle system with an inhomogeneous term given by the optical field. In order to define the differential operator in the square root we write Eq. (28) as follows whereα i ,β are matrices, γ a spinor and the mass parameter is M = ∆/(2v 2 F ). The dispersion in Eq. (29) must satisfy Thus we must have This algebra is satisfied by taking (α 1 ,α 2 ) = (σ x , σ y ) = σ T andβ = σ z . To determine the spinor γ we can switch off the Coulomb potential and impose that each component q µ satisfy a Klein-Gordon equation with an inhomogeneous term µ k E(t). To accomplish this it is sufficient to take γ = 1 1 .
Taking the Fourier transform of equation Eq (29) we get We shall now solve the homogeneous stationary eigenvalue problem Where ν here indicates a set of quantum numbers. The solution to this equation is known (see for example [27]) and we report it here for completeness. It is convenient to solve the problem in polar coordinate. Writing the spinor in the form We can write Eq. (33) in spinor components (we will drop for the moment the quantum numbers subscript) where we introduced the "Compton" energy The coupling constant α c is computed under the random phase approximation i.e α c = α0 εRPA with α 0 = e 2 /(v F ). In polar coordinates we have i∂ x ± i∂ y = e iφr (∂ r ± 1/rL z ) where L z is the third component of the angular momentum operator L = (0, 0, L z ). We can write now the spinor in Eq. (34) in the form Where j = m + 1/2 is the eigenvalue of the "isospin-angular" momentumĴ z =L z + 1 2 σ z along the z axis. Substituting the spinor (36) in Eq. (33) we get the radial equation In spite of its apparent simplicity the solution of this system is far from being trivial. One way to solve it is to decouple it in two independent second order differential equations, to do so we write the solution in the following form looking at this system at ρ = 0 we derive the useful relation Using (41) we can easily decouple the system, obtaining These equations are in the confluent hypergeometric form One solution is given by Kummer hypergeometric function the 1 F 1 (a, b; z) [3]. Thus Q 1 and Q 2 are given by Since 1 F 1 (a, b, 0) = 1 looking again at (40) for ρ = 0 we get Thus we can write the eigenfunctions for the discrete spectrum (ν = n, j) The bound states occur when the first argument of the hypergeometric function is a negative integer i.e.
γ − α c E ν Substituting this prescription in the eigenfunctions after lenghty but straightforward calculations [2] we get Both the discrete and continuum eigenstates derived here are used in the paper for the calculation of the Elliot formula.

THE OPTICAL MATRIX ELEMENT
It is useful to compute the dipole moment in the position space. Using the more manageable Compton units we get This integral converges and can be solved analytically. With Mathematica c we get µ(r) = e L 2 (2π) 2 k c π 2 cos(φ r )[I 1 (k c r) − S −1 (k c r)] − i √ πk c sin(φ r )G 0, − 1 2 , Where G(a n , b n |z) is the Mejer G function and S α (z) is the the Struve function of order α. In the case of zero gap (k c → 0) the integral simplifies considerably giving The dipole moment in eq (54) is essential to compute the oscillator strength |I ν | 2 defined in the integral (11) of the main text. This integral is not trivial to compute due to the special functions in the dipole moment. In particular even if both the real and imaginary part of µ(r) converge [3], the numerical implementation is problematic in particular in the region r >> 1. A way to solve this problem is to find a fit for these functions, in the large r region, involving just simple polynomials and exponentials. Using Mathematica c we find the following best fit where with µ R,I (r) we indicate the radial components of the real and imaginary part of the dipole moment. To compute the oscillator strength we thus use the exact dipole moment between 0 and k c r cut = 2 and the fit between k c r cut = 2 and k c r cut = 30. We verified that integrating up to k c r cut = 30 was enough for the numeric integral to converge.

COMPARISON WITH EXPERIMENTS
We shall now compare the results of our theoretical predictions with previous experimental results. To do so we compute the real and imaginary part of the dielectric function for MoSe 2 and MoS 2 on a fused silica substrate and we compare our calculations with the experimental results in [16]. The complex dielectric function is simply given by where b is the background dielectric constant of the substrate, here we assume a thickness d = 1 nm for the monolayers.  figure 3a,i,c,k of [4]