Optical Transparency in an effective model for Graphene

Motivated by experiments confirming that the optical transparency of graphene is defined through the fine structure constant and that it could be fully explained within the relativistic Dirac fermions in 2D picture, in this article we investigate how this property is affected by next-to-nearest neighbor coupling in the low-energy continuum description of graphene. A detailed calculation within the linear response regime allows us to conclude that, somewhat surprisingly, the zero-frequency limit of the optical conductivity that determines the transparency remains robust up to this correction.


Introduction
Graphene is a two-dimensional allotrope of carbon, arranged as a honeycomb lattice with a C 3v ⊗ Z 2 symmetry [1] that determines its remarkable physical properties [2,3,4,5]. In particular, the electronic spectrum arising from an atomistic tightbinding model displays two non-equivalent points K + , K − where the conduction and valence bands touch, and in whose vicinity the dispersion relation is approximately linear. This leads to an effective, low-energy continuum model where the electronic properties of the material are well captured by those of relativistic Dirac fermions in 2D. Among the plethora of physical consequences of this fact that have been already predicted and measured [2,3,4,5,6,7,8], we noticed an interesting experiment that measures the optical transparency of single and few-layer graphene [9]. The transparency is a physical property that is determined by the optical conductivity, i.e. the linear response to an electromagnetic field, in the zero-frequency limit. A variety of experiments confirm [9,10,11,12,13] that the measured transmittance is indeed compatible with the effective single-particle model of relativistic Dirac fermions in graphene. A number of different theoretical works have exploited this fact to calculate the light absorption rate in graphene from a "relativistic" quantum electrodynamics perspective [14,15,16,17,18,19,20,21,22]. An interesting question that remains open is up to what extent this effective model is valid in the representation of this optical property, since it arises from a tight-binding microscopic atomistic model that involves only the nearest neighbor hopping. In this article, we decided to explore what is the contribution to the optical conductivity arising from the next-to-nearest neighbors coupling in the atomistic Hamiltonian, included as a quadratic correction to the kinetic energy operator within the continuum effective model for graphene. Such a model has been considered in Ref. [23] to fully account for the Anomalous Integer Quantum Hall Effect in this material and the underlying wave equation is referred to in literature as Second Order Dirac Equation [24]. For our purposes, let us recall that within the linear response theory, general Kubo relations allow to express the transport coefficients in terms of retarded correlators [25], that for a pair of observablesÔ 1 ,Ô 2 are defined by (ζ = ± for Bosons and Fermions, respectively) These retarded correlators differ from the usual time-ordered ones that, by construction, are obtained via functional differentiation of the standard generating functional constructed form a path-integral formulation in quantum field theory. This rather technical inconvenience can be overcome by connecting the different propagators using a Lehmann representation, or alternatively to work in the Matsubara formalism at finite temperature and use analytic continuation a posteriori [25]. There is however a third, and more direct alternative, which is to express the generating functional in the contour time path (CTP), also known as Keldysh formalism in the condensed matter literature [26,27]. In this work, we choose the CTP formalism to explicitly calculate the polarization tensor as a retarded correlator of the current operators, which provides the correct definition of the optical conductivity within linear response theory. With these ideas in mind, we have organized the remaining of this article as follows: In Sect. 2, we present the details of the model. In Sect. 3 we present the Keldysh formalism to calculate the current-current correlator and in Sect. 4 we obtain the optical conductivity from the vacuum polarization tensor. We discuss our findings in Sect. 5. Some calculational details are presented in an Appendix.

Lagrangian, conserved current and generating functional
Graphene consist in one atom thick membrane of tightly packed carbon atoms in a honeycomb array. Its crystal structure, sketched in Fig. 1, is described in terms of two overlapping triangular (Bravais) sublattices so that for a given atom belonging to any of these sublattices, its nearest neighbors belong to the second sublattice, the nextto-nearest neighbors to the original sublattice and so on. The band structure at the next-to-nearest approximation is of the form where t and t are the nearest and next-to-nearest hopping parameters and where a 1.42Å is the interatomic distance. The points K + and K − at which f (K ± ) = 0 define the so-called Dirac points. Around K + , with tan(ϑ) = k y /k x . Around K − one merely has to replace ϑ → −ϑ in Eq. (2.3). The isotropic portion of this model was first considered in Ref. [23] as a natural framework to explain the Anomalous Integer Quantum Hall Effect in graphene. The anisotropic term in that work was treated perturbatively and shown not to contribute to the energy spectrum at first order. In the presence of electromagnetic interactions, the model is described by the Lagrangian [23] where θ = mv F . Here, ψ † and ψ are regarded as independent fields whose equations of motion are derived from the variation of the action with respect to these fields, namely, and similarly for ψ.
The Lagrangian in Eq.(2.4) remains invariant against the local change in the dynamical variables and the external electromagnetic field that is, it has a U (1) gauge symmetry. Noether's Theorem leads to the existence of the locally conserved current δψ . (2.7) The corresponding charge density is Figure 2. (Color online) The contour γ = γ + ⊕ γ − is depicted in the figure. The double folding of the time axis is displayed, by showing that always two points t − and t + , located in the time ordered t − ∈ γ − and and anti-time ordered t + ∈ γ + branches of the contour correspond to the same chronological time instant t. and the current density It is straightforward to verify, from the equations of motion, that j µ is conserved, Notice also that we can write With these ingredients, we can formulate the corresponding current-current correlator.

Generating functional in the Contour Time Path.
We seek to calculate the polarization tensor, defined as a retarded current-current correlator that, in linear response, determines the optical conductivity. For that purpose, we choose to represent the field-theory described in the previous section on the Contour Time Path (CTP) [26,27]. Let us define the contour γ = γ − ⊕ γ + , where γ − represents the time-ordered branch while γ + the anti-time-ordered branch, as depicted in Fig.2. Therefore, we define a contour evolution parameter τ ∈ γ, such that Also notice that, as depicted in Fig.2, both t + and t − have a unique correspondence to a given chronological instant of time t ∈ IR. Correspondingly, for operators and fields defined with their time arguments along the CTP, we have the definitions Then, the generating functional of (current) Green's functions of this two-dimensional system, defined on the CTP reads is the effective contribution to the action for the electromagnetic field. The path-integral on the CTP induces by construction the contour-ordering between the fields, defined by the operation T between two operatorsÔ 1 (τ ) andÔ 2 (τ ) in the Heisenberg picture (ζ = ± for Bosons/Fermions, respectively) Here, we have defined the contour Heaviside function as with the symbol > c indicating the relation "later than in the contour". In general physical situations where the sources and external fields do not break time-reversal invariance, ψ − (x, t) = ψ + (x, t), and the CTP becomes just a useful trick to express at once all the different correlators. Consider for instance the contour-ordered correlator between two fields, This single definition, depending on the location of the parameters τ 1 , τ 2 ∈ γ, generates at once four different propagators: (3.10) Here, we have defined the usual time-orderT and anti-time-orderT operators. Notice that not all correlators are independent, since they satisfy It is customary to organize the correlators above in the matrix form Using the definitions above, the retarded and advances correlators can be expressed as linear combinations of the previous ones From the CTP generating functional defined in Eq. (3.3), it is possible to generate the average current components while the second functional derivative gives the current-current correlation function, where the first term is the diamagnetic contribution [28] δj 17) and the others are the paramagnetic ones. We take the currents in normal order with respect to the fermionic field, so that J µ [A = 0] = 0. The linear response of the system to the external electromagnetic field is described by the second derivative in Eq. (3.16) evaluated at A µ = 0 [28], Then, the density response is The spatial components of the current are given by Here, we have defined the differential operatorŝ Applying Wick's theorem on the CTP for the definition of the current-correlator (correlators associated to disconnected diagrams vanish): The previous relation allows us to define the corresponding components of the polarization tensor in the CTP contour indices α, β = ±, The retarded component of the polarization tensor is obtained following the general prescription explained in Eq.(3.14), In terms of Fourier transforms, we have Here, the different propagators for the Hamiltonian model considered are, in Fourier space (F: Feynman, R: Retarded, A: Advanced), (3.29) In particular, for the linear response theory, we need the retarded component of the polarization tensor Here, the Fourier transform of the retarded component is given by where the different terms are defined by and a similar expression for Γ l cd (p + 2q). Below we obtain the polarization tensor explicitly.

The polarization tensor
The polarization tensor Π kl (p) contains the information about the conductivity on the plane of this two-dimensional system and also about its properties of transmission of light through it [16,28]. We are interested in the consequences of the application of harmonic homogeneous electric fields which, in the temporal gauge, are related with the vector potential by E k = −∂A k /∂t = iωA k . Since the conductivity is determined by the linear relation between the current and the applied electric field, J k = σ kl E l , from Eqs. (3.15), (3.18) and (3.31), we can write for the conductivity as a function of the frequency [16,28] σ kl (ω) = Π R kl (p) p 0 p→(ω,0) . (4.1) So, in the following we evaluate Π R kl (ω, 0) from Eq.(3.31), that hence requires the evaluation of the three integrals defined in Eq.(3.32). Let us start with Π F A kl (p), Specializing this expression to the case p = (ω, 0), we write By writing q 1 = Q cos ϕ, q 2 = Q sin ϕ, and noticing that the denominator is independent of ϕ, it is straightforward to get for the trace in the numerator integrated over ϕ, ,(4.5) for k, l = 1, 1 or 2, 2, and a vanishing result for k, l = 1, 2 or 2, 1.
Since the previous result is a quadratic polynomial in q 0 , and the denominator in Eq. (4.4) is a quartic expression in the integration variable, the integral over q 0 can be done on the complex plane, taking into account the position of the simple poles of the integrand with respect to the real axis. With a dimensional regularization of the resulting integral (dimension d = 2 − s), as described in detail in Appendix, one finds A similar procedure, as described in Appendix, leads to the corresponding expressions for the other two pieces of the retarded polarisation tensor We notice that the three separate parts above, which together yield the retarded polarization tensor, display a pole at s = 0. However, when added together according to Eq.(3.31), the poles exactly cancel to yield a finite result The result above must be multiplied by a factor of 2 due to the spin degeneracy, and another factor of 2 due to valley degeneracy. Thus, the final result of the optical conductivity is (4.10) Remarkably, this is the same result that is obtained for the usual linear dispersion approximation. Therefore, we conclude that the optical conductivity, and hence the transparency in graphene are not affected by next-to-nearest neighbor contributions to the tight-binding microscopic model, that translate into a quadratic correction to the kinetic energy, as considered in this work.

Conclusions
Among the many outstanding properties of graphene which can be described within the Dirac limit, its optical transparency is entirely explained in terms of the fine structure constant. A natural question is to ask the extent at which such picture deviates from the experimental measurements. In this regard, in this article we considered the next-tonearest neighbors contribution which in the continuum corrects the kinetic term with a quadratic contribution. Introducing the CTP formalism, we calculate the linear response current-current correlator from which the optical conductivity is derived. Within this formalism, it is straightforward to obtain the retarded part of the polarization tensor after a dimensional regularization of the involved integrals. Remarkably and somehow unexpectedly, we found the conductivity of the Dirac limit to be robust against quadratic corrections. This encouraging results opens the possibility of testing deviations of the Dirac limit in graphene in other physical phenomena. These results are currently under scrutiny and results will be reported elsewhere.

Appendix: Regularization of the momentum integrals
In this Appendix, we present in detail the dimensional regularization method used to calculate the momentum integrals defined in the main text. Let us consider the term Eq. (4.3). After taking the trace and performing the angular integral as shown in Eq. (4.5), we have to evaluate with B F A given in Eq. (4.4) with |q| = Q. Clearly, on the q 0 -plane, the integrand has three poles on the positive imaginary plane at q and a single pole on the negative imaginary plane at q We evaluate the q 0 integral by means of the residue theorem, closing the contour on the lower plane. The result of this procedure can be expressed as Here, we have defined the numerator as the polynomial function 3) By simply counting powers in numerator and denominator, it is clear that the remaining integral is divergent and needs regularization. For this purpose, we first perform a partial fraction decomposition of the denominator as follows After this, the integral splits into two contributions For each integral (i.e. Q j = Q 1 , Q 2 respectively), let us analyze the asymptotic behavior of the integrand at large momentum values, say for Q > Q * , with Q * an arbitrary but large momentum scale. In this regime, where we have defined P Asymp F A (Q, ω, Q j ) as the polynomial obtained by truncating the asymptotic expansion above up to O[Q −3 ], for Q > Q * . Therefore, using this expansion, we regularize each of the integrals using the prescription (d = 2 − s) (5.7) After lengthy but straightforward algebra, we obtain in the limit → 0 + Let us now consider the expression for Π RF 11 (ω, 0), as obtained after calculating the trace and angular integration according to Eq.(4.5) In this case, on the q 0 -plane the integrand has three poles on the negative imaginary plane, q Therefore, we calculate the integral over q 0 using the residue theorem, by choosing a contour that closes on the upper complex plane. Thus, .
The numerator of the resulting integrand is defined by the quartic polynomial function ) and then the integral is clearly divergent. A consistent regularization procedure is applied in this case as well. By performing the same partial fraction expansion of the denominator, as in Eq.(5.4), we find that the integral splits into two pieces ( For each integral (i.e. Q j = Q 1 , Q 2 respectively), we analyze the asymptotic behavior of the integrand at large momentum values, say for Q > Q * . In this regime, where we have defined P Asymp RF (Q, ω, Q j ) as the polynomial obtained by truncating the asymptotic expansion above up to O[Q −3 ], for Q > Q * . Therefore, using this expansion, we regularize each of the integrals using the prescription (d = 2 − s) After straightforward manipulations, we obtain in the limit → 0 + Finally, let us consider the term Eq.(3.32). After tracing and performing the angular integral, Clearly, on the q 0 -plane, the integrand has two poles on the positive imaginary plane at q ± v F Q, and two poles on the negative imaginary plane at q We evaluate the q 0 integral by the residue theorem, closing the contour on the upper plane. Thus, Π RA 11 (ω, 0) = e 2 4m 2 ∞ 0 dQ (2π) 3 Q P RA (Q, ω) (ω + 2i ) (ω − 2v F Q + 2i ) (ω + 2v F Q + 2i ) .
The numerator of the resulting integrand is defined by the quartic polynomial function and hence the diverging integral needs also a regularization. As in the former two cases, we first do a partial fraction decomposition of the denominator, to obtain where we have defined Q 3 = (ω + 2i )/(2v F ), Q 4 = −Q 3 . For each integral (i.e. Q j = Q 3 , Q 4 respectively), we analyze the asymptotic behavior of the integrand at large momentum values, say for Q > Q * . In this regime, where we have defined P Asymp RA (Q, ω, Q j ) as the polynomial obtained by truncating the asymptotic expansion above up to O[Q −3 ], for Q > Q * . Therefore, using this expansion, we regularize each of the integrals using the prescription (d = 2 − s)