Low-energy theory for strained graphene: an approach up to second-order in the strain tensor
Maurice Oliva-Leyva Chumin Wang










moliva@iim.unam.mx

chumin@unam.mx

Instituto de Investigaciones en Materiales, Universidad Nacional Autónoma de
México, Apartado Postal 70-360, 04510 Mexico City, Mexico.






An analytical study of low-energy electronic excited states in an uniformly strained graphene is carried out up to second-order in the strain tensor. We report an new effective Dirac Hamiltonian with an anisotropic Fermi velocity tensor, which reveals the graphene trigonal symmetry being absent in low-energy theories to first-order in the strain tensor. In particular, we demonstrate the dependence of the Dirac-cone elliptical deformation on the stretching direction respect to graphene lattice orientation. We further analytically calculate the optical conductivity tensor of strained graphene and its transmittance for a linearly polarized light with normal incidence. Finally, the obtained analytical expression of the Dirac point shift allows a better determination and understanding of pseudomagnetic fields induced by nonuniform strains.





Low-energy theory for strained graphene: an approach up to second-order in the strain tensor
    Chumin Wang
    
============================================================================================





§ INTRODUCTION


Given the striking interval of elastic response of graphene <cit.>, can withstand a reversible stretching up to 25 %, strain engineering has been widely used to improve and/or to tune its electronic, thermal, chemical and optical properties <cit.>. For instance, theoretical predictions have been made of a band-gap opening by large uniaxial strains from both tight-binging approach <cit.> and density functional theory <cit.>, whenever the strain produces such a Hamiltonian modification beyond the inequalities obtained by Hasegawa, et al. <cit.> The emergence of the pseudomagnetic field caused by a nonuniform strain is possibly the most interesting strain-induced electronic effect, due to the possibility of observing a pseudoquantum Hall effect under zero external magnetic fields <cit.>. Nowadays, the transport signatures of the such fictitious fields are actively investigated <cit.>. Moreover, from a view point of basic research, strained graphene opens an opportunity to explore mixed Dirac–Schrödinger Hamtiltonian <cit.>, fractal spectrum <cit.>, superconducting states <cit.>, magnetic phase transitions <cit.>, metal-insulator transition <cit.>, among others exotic behaviours.


The concept of strain engineering has been also extended to the optical context <cit.>. The optical properties of graphene are ultimately provided by its electronic structure, which can be modified by strain. For example, pristine graphene presents a transparency defined by fundamental constants, around 97.7%, over a broad band of frequencies <cit.>. This remarkable feature is essentially a consequence of its unusual low-energy electronic band structure around the Dirac points. Under uniform strain, such conical bands are deformed which produces anisotropy in the electronic dynamics <cit.>. Accordingly, this effect gives rise an anisotropic optical conductivity of strained graphene <cit.> and, therefore, a modulation of its transmittance as a function of the polarization of the incident light, as experimentally observed <cit.>. From a theoretical viewpoint, this optoelectronic behaviour of strained graphene has been quantified by continuum approaches up to first-order in the strain tensor <cit.>. However, nowadays there are novel methods for applying uniaxial strain larger than 10% in a nondestructive and controlled manner <cit.>. So, a low-energy continuum theory for the electronic and optical properties of strained graphene, up to second-order in the strain tensor, seems to be needed <cit.>.


In this paper, we derive the effective Dirac Hamiltonian for graphene under uniform strain up to to second-order in the strain tensor. For this purpose, we start from a nearest-neighbor tight-binding model and carry out an expansion around the real Dirac point. Unlike previous approaches to the first-order in strain, we show how the obtained low-energy Hamiltonian reveals the trigonal symmetry of the graphene. Also, we calculate the optical conductivity of strained graphene and characterize its transmittance for a uniaxial strain up to second-order in the stretching magnitude.
These findings describe in a more accurate form the electronic and optical properties of strained graphene and, hence, can be potentially utilized towards novel optical characterizations of the strain state of graphene.






§ TIGHT-BINDING MODEL AS STARTING POINT


Strain effects on electronic properties of graphene are usually captured by using a nearest-neighbor tight-binding model <cit.>. Within this approach, one can demonstrate that the Hamiltonian in momentum space for graphene under a uniform strain is given by <cit.>

    H(k)=-∑_n=1^3 t_n(
    [           0 e^-ik·δ_n^';  e^ik·δ_n^'           0 ]),

where the strained nearest-neighbor vectors are obtained by δ_n^'=(I̅+ϵ̅)·δ_n, being I̅ the (2×2) identity matrix and ϵ̅ the rank-two strain tensor, whose components are independent on the position. Here, we choose the unstrained nearest-neighbor vectors as

    δ_1=a_0/2(√(3),1),  δ_2=a_0/2(-√(3),1), δ_3=a_0(0,-1),

where a_0 is the intercarbon distance for pristine graphene. Thus, the x (y) axis of the Cartesian coordinate system is along the zigzag (armchair) direction of the honeycomb lattice.
Owing to the changes in the intercabon distance, the nearest-neighbor hopping parameters are modified. Here we consider this effect by means of the commonly used model <cit.>

    t_n=t_0e^-β(|δ_n^'|/a_0-1),

where t_0=2.7 is the hopping parameter for pristine graphene and β≈3.


From equation (<ref>) follows that the dispersion relation near the Fermi energy of graphene under uniform strain is given by two bands,

    E(k)=±| t_1e^ik·δ_1^' + t_2e^ik·δ_2^' + t_3e^ik·δ_3^'|,

which remains gapless as long as the triangular inequalities, | t_1-t_2|≤| t_3|≤| t_1+t_2|, are satisfied <cit.>. Evaluating equation (<ref>) for uniaxial strains, V. Pereira, et al., found the minimum uniaxial deformation that leads to the gap opening is about 23% <cit.>. This result is confirmed by the ab initio calculations, finding that this gap in strained graphene requires deformations larger than 20% <cit.>. Therefore, the use of an effective Dirac Hamiltonian obtained from equation (<ref>) is justified for uniform deformations up to the order of 10%.

For this purpose, it is important to take into account a crucial detail: the strain-induced shift of the Dirac points in momentum space. In absence of deformation, the Dirac points K_D (determined by condition E(K_D)=0) coincide with the corners of the first Brillouin zone. Then, to obtain the effective Dirac Hamiltonian in this case, one simply expand the Hamiltonian (<ref>) around such corners, e.g., K_0=(4π/3√(3)a_0,0). However, in presence of deformations, the Dirac points do not coincide even with the corners of the strained first Brillouin zone <cit.>. Thus, to obtain the effective Dirac Hamiltonian, one should no longer expand the Hamiltonian (<ref>) around K_0. As demonstrated <cit.>, such expansion around K_0 yields an incorrect derivation of the anisotropic Fermi velocity. The appropriate procedure is to find first the new positions of the Dirac points and then carry out the expansion around them <cit.>.




§ EFFECTIVE DIRAC HAMILTONIAN



As first step, we determine the new positions of Dirac points from the condition E(K_D)=0, up to second order in the strain tensor, which is the leading order used throughout the rest of the paper. Essentially, we calculate the strain-induced shift of the Dirac point K_D from the corner K_0 of the first Brillouin zone by using equation E(K_D)=0, which leads to

    ∑_n=1^3t_ne^iK_D·δ_n^'=
     ∑_n=1^3t_ne^iK_D· (I̅+ϵ̅)·δ_n=
     ∑_n=1^3t_ne^iG·δ_n=0,

where G≡(I̅+ϵ̅)·K_D is the effective Dirac point. As demonstrated in Appendix <ref>, G can be expressed as

    G=K_0 + A^(1) + A^(2) + 𝒪(ϵ̅^3),

where

    A_x^(1) + i A_y^(1)=β/2a_0(ϵ_xx-ϵ_yy - 2iϵ_xy),

and

    A_x^(2) + i A_y^(2)=β(4β+1)/16a_0(ϵ_xx-ϵ_yy + 2iϵ_xy)^2.





Notice that the correction up to first order, A^(1), coincides with the value previously reported <cit.>, which is interpreted as a gauge field for nonuniform deformations <cit.>. On the other hand, the expression (<ref>) for the second-order correction A^(2) is one of the main contributions of this work. To demonstrate its relevance, we numerically calculate the positions of G for two deformations and compare them with the analytical results given by (<ref>-<ref>). As illustrated in Fig. <ref>(a) for a uniaxial strain along zigzag direction, the values of G_x estimated up to first order in the strain magnitude ϵ (blue solid circles) clearly differ from the exact numerical values of G_x (gray line) as ϵ increases, while the values of G_x estimated up to second order (red open circles) show a significantly better approximation. The case of a shear strain is an even more illustrative example of the relevance of A^(2). According to the first-order correction, G_x does not change under a shear strain, which is at variance with the exact numerical result displayed in Fig. <ref>(b). In contrast, the values of G_x estimated up to second order present a good agreement with the numerical values over the studied range of ϵ. Beyond the present work, the second-order correction A^(2) for nonuniform strain could be relevant to a more complete analysis of the strain-induced pseudomagnetic fields. For example, in presence of a deformation field given by u=(u(y),0), for which ϵ̅_xx=ϵ̅_yy=0 and ϵ̅_xy=∂_yu(y)/2, the pseudomagnetic field B_ps, derived from the standard expression B_ps=∇×A^(1), results equal to zero. However, if A^(2) is taken into account by means of the possible generalized expression B_ps=∇×A^(1) + ∇×A^(2), one can demonstrate that the resulting pseudomagnetic field B_ps is not zero. The implications of this issue will be discussed with details in an upcoming work.

Knowing the position of the Dirac point K_D, through equation (<ref>), one can now proceed to the expansion of Hamiltonian (<ref>) around K_D, by means of k=K_D+q, to obtain the effective Dirac Hamiltonian. Following this approach up to second order in the strain tensor ϵ̅, the effective Dirac Hamiltonian can be written as (see Appendix <ref>)

    H=ħ v_0τ·(I̅ + ϵ̅ - βϵ̅ -βϵ̅^2 - βκ̅_1 + β^2κ̅_2)·q,

where v_0=3t_0a_0/2ħ is the Fermi velocity for pristine graphene, τ=(τ_x,τ_y) is a vector of (2×2) Pauli matrices describing the pseudospin degree of freedom,

    κ̅_1=1/8([     (ϵ_xx-ϵ_yy)^2 -2ϵ_xy(ϵ_xx-ϵ_yy); -2ϵ_xy(ϵ_xx-ϵ_yy)           4ϵ_xy^2 ]),

and

    κ̅_2=1/4([ ϵ_xx^2-ϵ_yy^2+2ϵ_xxϵ_yy+2ϵ_xy^2                       4ϵ_xxϵ_xy;                       4ϵ_xxϵ_xy                2(ϵ_yy^2-ϵ_xy^2) ]).


It is important to emphasize that the explicit form of equations (<ref>) and (<ref>) is a consequence of the Cartesian coordinate system xy chosen. For an arbitrary coordinate system x̃ỹ, rotated by an angle ϑ respect to the system xy, the new expressions for κ̅_1 and κ̅_2 should be found by means of the transformation rules of a second order Cartesian tensor <cit.>.

From  equation (<ref>) one can recognize the Fermi velocity tensor as

    v̅=v_0(I̅ + ϵ̅ - βϵ̅ -βϵ̅^2 - βκ̅_1 + β^2κ̅_2),

which generalizes the expression, v_0(I̅ + ϵ̅ - βϵ̅), for the Fermi velocity tensor up to first-order in the strain tensor reported  in Refs. [Oliva13,Volovik14a].

As a consistency test, let us consider an isotropic uniform strain of the graphene lattice, which is simply given by ϵ̅=ϵI̅. Under this deformation, the new intercarbon distance a is rescaled as a=a_0(1+ϵ), whereas the new hopping parameter t, expanding equation (<ref>) up to second order in strain, results t=t_0(1-βϵ+β^2ϵ^2/2). Therefore, the new Fermi velocity, v=3t a/2ħ, obtained straight away from the nearest-neighbor tight-binding Hamiltonian, takes the value v=v_0(1-βϵ+ϵ-βϵ^2+β^2ϵ^2/2). This result can be alternatively obtained by evaluating our tensor (<ref>) for ϵ̅=ϵI̅.

The tensorial character of v̅ is due to the elliptic shape of the isoenergetic curves around K_D. Notice that the principal axes of the Fermi velocity tensor up to first-order in the strain tensor, v_0(I̅ + ϵ̅ - βϵ̅), are collinear with the principal axes of ϵ̅. Therefore, within the effective low-energy Hamiltonian up to first-order in the strain tensor, the anisotropic electronic behaviour is only originated from the strain-induced anisotropy. Nevertheless, the terms κ̅_1 and κ̅_2 in equation (<ref>) suggest that the second-order deformation theory might reveal the anisotropy (trigonal symmetry) of the underlying honeycomb lattice.



To clarify this issue, let us to consider graphene subjected a uniaxial strain such that the stretching direction is rotated by an arbitrary angle θ respect to the Cartesian coordinate system xy (see Fig. <ref>). In this case, the strain tensor (ϵ̅) in the reference system xy reads

    ϵ̅(θ)=ϵ(
    [ cos^2θ - νsin^2θ    (1+ν)cosθsinθ;    (1+ν)cosθsinθ sin^2θ - νcos^2θ ]),

where ϵ is the strain magnitude. Note that both ϵ̅(θ) and ϵ̅(θ+180^∘) represent physically the same uniaxial strain, which can be confirmed in equation (<ref>). It is important to mention that for θ=n 60^∘ (θ=90^∘+n 60^∘), being n an integer, the stretching is along a zigzag (armchair) direction of graphene lattice.


As discussed above, under the strain (<ref>), the Fermi velocity tensor up to first-order in the strain tensor, v_0(I̅ + ϵ̅ - βϵ̅), is diagonal in the coordinate system x'y', rotated by the angle θ respect to the coordinate system xy. However, the Fermi velocity tensor (<ref>), up to second-order in the strain tensor, is diagonal in a coordinate system x”y”, rotated by an angle θ_v such that

    tan2θ_v=2v_xy/v_xx-v_yy,

which determines the direction of lower electronic velocity. In the reciprocal space, the angle θ_v characterizes the pulling direction of isoenergetic curves, i.e., the principal axis of the isoenergetic ellipses, as illustrated in Fig. <ref>(d).

In Fig. (<ref>), we show the difference θ = θ_v-θ, numerically calculated from equation (<ref>), as a function of the stretching direction θ for two different strain magnitudes ϵ=5% and 10%. The observed six-fold behaviour of θ can be analytically evaluated by

    θ   ≈    -β(2β+1)(1+ν)/16(β-1)ϵsin(6θ) 
          ×(1-β(1-ν)/2ϵcos(6θ))

in good agreement with the numerical values, as shown in Fig. (<ref>). From the last expression, it follows that the principal axes of the Fermi velocity tensor (<ref>) are only collinear with the principal axes of ϵ̅(θ) for θ=n 30^∘, i.e., when the stretching is along the zigzag or armchair crystallographic directions. This result demonstrates that our Hamiltonian (<ref>), a second-order deformation theory, reveals the trigonal symmetry of underlying honeycomb lattice.






§ OPTICAL PROPERTIES


An anisotropic Dirac system described by the effective Hamiltonian

    H=ħ
    v_0τ·(I̅+Δ̅)·q,

being Δ̅ a symmetric (2×2) matrix such that Δ_ij≪1, presents an anisotropic optical response captured by the conductivity tensor (see Appendix <ref>):

    σ̅(ω)    ≈   σ_0(ω){I̅ - (Δ̅)I̅ + 2Δ̅ + Δ̅^2
       +   1/2[((Δ̅))^2
    + (Δ̅^2)]I̅-2(Δ̅)Δ̅},

where ω is the frequency of the external electric field and σ_0(ω) is the optical conductivity of the unperturbed Dirac system, i.e., the optical conductivity of unstrained graphene. Equation (<ref>) is a generalization up to second-order in Δ̅ of previous expression until first-order in Δ̅ for the optical conductivity of an anisotropic Dirac system, as it can be seen in equation (17) of Ref. [Oliva2016].


Now, comparing equations (<ref>) and (<ref>), the optical conductivity tensor σ̅(ω) of strained graphene is  straightforward obtained by making the replacement:

    Δ̅= ϵ̅ - βϵ̅ -βϵ̅^2 - βκ̅_1 + β^2κ̅_2,

into equation (<ref>). Regarding terms up to second-order in the strain tensor, it results

    σ̅(ω)   =   σ_0(ω) [ I̅ + β̃(ϵ̅)I̅ - 2β̃ϵ̅
       +   (5β+2β̃^2/4(ϵ̅^2) + 4β̃^2-2β^2-β/8((ϵ̅))^2)I̅
       +    (β̃^2-2β)ϵ̅^2 -2β̃^2(ϵ̅)ϵ̅  - 2βκ̅_1 +2β^2κ̅_2],

where β̃=β-1. This equation generalizes previous works <cit.>, in which the optical conductivity of graphene under uniform strain was reported up to first-order in the strain tensor. 

Let us make a proof about the consistency of equation (<ref>). When graphene is at half filling, i.e., the chemical potential
equals to zero, the optical conductivity σ_0(ω) is frequency-independent and is given by the universal value e^2/(4ħ) <cit.>. It is important to emphasize that this result is independent on the value v_0 of the Fermi velocity <cit.>. Therefore, under an isotropic uniform strain ϵ̅=ϵI̅, which only leads to a new isotropic Fermi velocity v=v_0(1-βϵ+ϵ-βϵ^2+β^2ϵ^2/2), the optical conductivity does not change and remains equal to σ_0=e^2/(4ħ), at least within the Dirac cone approximation <cit.>. In other words, any expression reported as optical conductivity tensor for uniformly strained graphene, as a function on the strain tensor, to be evaluated for ϵ̅=ϵI̅ must give rise σ_0I̅, as occurred when one evaluates the tensor (<ref>).

The optical conductivity up to first-order in the strain tensor, σ_0[I̅ + β̃(ϵ̅)I̅ - 2β̃ϵ̅], under a uniaxial strain (<ref>) can be characterized by σ_∥=σ_0[1-β̃ϵ(1+ν)] and σ_⊥=σ_0[1+β̃ϵ(1+ν)], where σ_∥ (σ_⊥) is the optical conductivity parallel (perpendicular) to the stretching direction (see blue lines in Fig. <ref>). Within the first-order approximation, the optical conductivity along the stretching direction decreases by the same amount that the transverse conductivity increases, independently of θ. This behaviour is modified when second-order terms are taken into account.

In Fig. <ref>(a), we plot the components of the optical conductivity tensor (<ref>) versus the stretching magnitude ϵ for a uniaxial strain along the armchair direction. The perpendicular conductivity to the stretching direction, σ_xx(ϵ), does not have appreciable difference respect to the lineal approximation σ_0[1+β̃ϵ(1+ν)] whereas the parallel conductivity, σ_yy(ϵ), noticeably differs from σ_0[1-β̃ϵ(1+ν)] with increasing strain. On the other hand, Fig. <ref>(b) displays a contrary behaviour of the optical conductivity for a uniaxial strain along the zigzag direction. For this case, the parallel conductivity, σ_xx(ϵ), looks slight different from σ_0[1-β̃ϵ(1+ν)] whereas the perpendicular conductivity, σ_yy(ϵ), is noticeably greater than σ_0[1+β̃ϵ(1+ν)] with increasing strain. This increase of σ_yy(ϵ) respect to σ_0[1+β̃ϵ(1+ν)] might help to give a better understanding of the change in the transmission of hybrid graphene integrated microfibers elongated along their axial direction <cit.>. For example, in Figure 2(b) of Ref. [Chen2016], it is possible to appreciate that the experimental data of this change gradually differ, with increasing strain, from the theoretical calculation using the first-order linear approximation σ_0[1+β̃ϵ(1+ν)], which can be improved by considering the second-order contribution as shown in Fig. <ref>(b). 



To complete our discussion about the emergence of the trigonal symmetry of graphene in the continuum approach presented here, we now study the transmittance of linearly polarized light on strained graphene. Considering graphene as a two-dimensional sheet with conductivity σ̅ and from the boundary conditions, vacuum-graphene-vacuum, for the electromagnetic field on the interfaces, the transmittance for normal incidence reads as <cit.>

    T=(1 + [σ_xxcos^2θ_i + σ_yysin^2θ_i + σ_xysin 2θ_i]/2ε_0c)^-2

where ε_0 is the vacuum permittivity, c is the speed of light in vacuum and θ_i is the incident polarization angle. Note that for a pristine graphene with σ̅=σ_0I̅, equation (<ref>) reproduces the experimentally observed constant transmittance T=(1+πα/2)^-2≈1-πα over visible and infrared spectrum <cit.>, being α≈1/137 the fine-structure constant.
From equation (<ref>) it can be seen that an anisotropic absorbance yields a periodic modulation of the transmittance as a function of the polarization direction θ_i <cit.>.


For the case of a uniaxial strain, and assuming the chemical potential
equal to zero, from equations (<ref>), (<ref>) and (<ref>) it follows that the transmittance up to second-order in the strain magnitude ϵ is given by

    T    =    1 - πα + παβ̃(1+ν)ϵcos2(θ_i-θ_v) - πα/2β̃^2(1+ν)^2ϵ^2
           + πα/2(1+γcos6θ)(1+ν)^2ϵ^2cos2(θ_i-θ_v),

where γ=β(2β+1)/4. Expression (<ref>) reveals two new remarkable features in comparison with the first-order theory. As illustrated in Fig. <ref>, the transmittance mean value, ⟨ T⟩= 1 - πα - παβ̃^2(1+ν)^2ϵ^2/2, has a negative shift with respect to the first-order average value T_0 = 1 - πα. Second, the transmittance oscillation amplitude (T) is determined by

    T (θ) = 2παβ̃(1+ν)ϵ + πα(1+γcos6θ)(1+ν)^2ϵ^2.


While the first-order expression for the transmittance oscillation amplitude, 2παβ̃(1+ν)ϵ, is independent on the stretching direction θ, T of equation (<ref>) depends on θ. For example, for a uniaxial strain along the zigzag (armchair) direction with θ=n 60^∘ (θ=90^∘ + n 60^∘), T takes its highest (lowest) value, as displayed in Fig. <ref>. This strectching direction dependent T might be used to confirm experimentally the present theory up to second-order in the strain tensor, as done for small strain less than 1% <cit.>.






§ CONCLUSION


We have analytically deduced a new effective Dirac Hamiltonian of graphene under a uniform deformation up to second-order in the strain tensor, including new Dirac-point positions that are qualitatively different from those predicted by first-order approaches, as occurred for the shear strain. Moreover, based on a detailed analysis about the anisotropic Fermi velocity tensor, we demonstrated how our second-order deformation theory reveals the trigonal symmetry of graphene unlike the previous first-order results.

We further derived, for the first time, analytical expressions for the high-frequency electric conductivity and light transmittance of a strained graphene up to second-order in the strain tensor. The magnitude of this transmittance oscillates according to the incident light polarization and the oscillation amplitude depends on the stretching direction, in contrast to the first-order prediction. In fact, within the first-order theory, the maximal transmittance occurs when the light polarization coincides to the stretching direction. However, the second-order theory predicts such coincidences only for stretching along zigzag and armchair directions. Therefore, the obtained light transmittance results can be experimentally verified by optical absorption measurements and they would be used for characterizing the deformation states of strained graphene. In general, the analytical study presented in this article has the advantage of being concise and establishes a reference point for upcoming numerical and experimental investigations.

It would be important to stress that the observed absence of lattice symmetry in the optical properties of strained graphene is due to the combination of the low-energy effective Dirac model and first-order approximation in the strain tensor. Such absence can be overcome by carrying out the study within the first-neighbour tight-binding model as occurred for high-energy electron excitations <cit.> or by introducing second-order effects in the strain tensor even within the simplest Dirac model, as done in this article. This finding of trigonal symmetry in optical response reveals the capability of low-energy effective Dirac theory to describe properly anisotropic electron behaviour in graphene under strong uniform deformations. However, a tight-binding model beyond nearest-neoghbour interactions would be required to analyze both the gap opening and the electron-hole spectrum symmetry induced by lattice strain <cit.>. Finally, the present work can be extended to perform an analytical study of the pseudomagnetic fields induced by nonuniform strains.




This work has been partially supported by CONACyT of Mexico through Project 252943, and by PAPIIT of Universidad Nacional Autónoma de México (UNAM) through Projects IN113714 and IN106317. Computations were performed at Miztli of UNAM. M.O.L. acknowledges the postdoctoral fellowship from DGAPA-UNAM.








§ DIRAC POINT POSITION


Here we provide the derivation of expressions (<ref>–<ref>) of main text.
Equation, E(K_D)=0, can be rewritten as

    ∑_n=1^3t_ne^iK_D·δ_n^'=
     ∑_n=1^3t_ne^iK_D· (I̅+ϵ̅)·δ_n=
     ∑_n=1^3t_ne^iG·δ_n=0,

where G≡(I̅+ϵ̅)·K_D is the effective Dirac point associated to a pristine honeycomb lattice with strained nearest-neighbor hopping integrals t_n. To solve equation (<ref>) in a perturbative manner, we cast the position of G as

    G=K_0 + A^(1) + A^(2) + 𝒪(ϵ̅^3),

where A^(1) ( A^(2)) is the correction from first (second) order in the strain tensor. Similarly, we consider Taylor expansions of t_n, up to second order in strain tensor, in the form

    t_n=t_0[1 + Δ_n^(1) +  Δ_n^(2) + 𝒪(ϵ̅^3)],

where Δ_n^(1) (Δ_n^(2)) are terms of the first (second) order in the strain tensor.

Substituting equations (<ref>) and (<ref>) into equation (<ref>), the coefficient of the first-order strain tensor should be equal to zero, which leads to

    ∑_n=1^3[Δ_n^(1)+ iA^(1)·δ_n]e^iK_0·δ_n=0.


Analogously, the coefficient of the second-order strain tensor should also be zero, yields

    ∑_n=1^3      [ Δ_n^(2) + iΔ_n^(1)A^(1)·δ_n - (A^(1)·δ_n)^2/2. 
          . + iA^(2)·δ_n] e^iK_0·δ_n=0.


From equation (<ref>), A^(1) can be determined and it is used as input of equation (<ref>) to obtain A^(2). To carry out this procedure, it is necessary to explicitly know Δ_n^(1) and Δ_n^(2)  as functions of the strain tensor.


Expanding t_n, up to second order in strain tensor, gives

    t_n/t_0   =   exp[-β(|δ_n^'|/a_0-1)]
       =   exp[-β( 1/a_0^2δ_n·ϵ̅·δ_n + 1/2a_0^4(ϵ̅·δ_n)^2..
          .. - 1/2a_0^4(δ_n·ϵ̅·δ_n)^2 + 𝒪(ϵ̅^3))]
       =    1 - β/a_0^2δ_n·ϵ̅·δ_n - β/2a_0^4(ϵ̅·δ_n)^2
           + β(β+1)/2a_0^4(δ_n·ϵ̅·δ_n)^2 + 𝒪(ϵ̅^3).

Then, by comparing equations (<ref>) and (<ref>) one obtains

    Δ_n^(1)= - β/a_0^2δ_n·ϵ̅·δ_n,

and 

    Δ_n^(2)= - β/2a_0^4(ϵ̅·δ_n)^2 + β(β+1)/2a_0^4(δ_n·ϵ̅·δ_n)^2.


Finally, substituting Δ_n^(1) into equation (<ref>), we get

    A_x^(1) + i A_y^(1)=β/2a_0(ϵ_xx-ϵ_yy - 2iϵ_xy),

and consequently, using this result and the expression of Δ_n^(2), equation (<ref>) can be rewritten as

    A_x^(2) + i A_y^(2)=β(4β+1)/16a_0(ϵ_xx-ϵ_yy + 2iϵ_xy)^2.


Note that equations (<ref>) and (<ref>) are the first- and second-order corrections to the Dirac point position given in equation (<ref>) of the main text.







§ EFFECTIVE DIRAC HAMILTONIAN


In order to derive the effective Dirac Hamiltonian given by equation (<ref>) in the main text, we start from the tight-binding model in momentum space for graphene under a uniform strain,

    H=-∑_n=1^3 t_n(
    [                0 e^-ik·(I+ϵ̅)·δ_n;  e^ik·(I+ϵ̅)·δ_n                0 ]),

and we consider momenta close to the Dirac point K_D, by means of the substitution k=K_D+q. Then, expression (<ref>) transforms as

    H= (
    [   0 h^∗;   h   0 ]),

where h=-∑_n=1^3 t_n e^i(K_D+q)·(I+ϵ̅)·δ_n.
Now, using equation (<ref>), h can be expanded up to first-order in q and second-order in ϵ̅ as


    h   =    -∑_n=1^3 t_n e^i[K_D·(I+ϵ̅)·δ_n+ q·(I+ϵ̅)·δ_n]
       ≈     -∑_n=1^3 t_n e^iK_0·δ_n e^i(A^(1)·δ_n + A^(2)·δ_n)e^i q·(I+ϵ̅)·δ_n
       ≈    -∑_n=1^3 t_n e^iK_0·δ_n[1 + i A^(1)·δ_n + i A^(2)·δ_n - (A^(1)·δ_n)^2/2][1 + iq·(I+ϵ̅)·δ_n] 
       ≈    -∑_n=1^3 t_n e^iK_0·δ_n[1 + i A^(1)·δ_n + i A^(2)·δ_n - (A^(1)·δ_n)^2/2 + iq·(I+ϵ̅)·δ_n. 
                               - . (A^(1)·δ_n)(q·(I+ϵ̅)·δ_n) - (A^(2)·δ_n)(q·δ_n) - i(A^(1)·δ_n)^2(q·δ_n)/2 ]

and substituting expression (<ref>) for t_n in equation (<ref>), the expansion of h results

    h   ≈    -t_0∑_n=1^3 e^iK_0·δ_n[ 1 + i A^(1)·δ_n + i A^(2)·δ_n - (A^(1)·δ_n)^2/2 + Δ_n^(1) + iΔ_n^(1)A^(1)·δ_n +  Δ_n^(2)_suming over n equal to zero. 
          + iq·δ_n_h_0+ iq·ϵ̅·δ_n_h_1,a- (A^(1)·δ_n)(q·δ_n) + iΔ_n^(1)q·δ_n_h_1,b-(A^(1)·δ_n)(q·ϵ̅·δ_n) + iΔ_n^(1)(q·ϵ̅·δ_n)_h_2,a
          . - i(A^(1)·δ_n)^2(q·δ_n)/2 - Δ_n^(1)(A^(1)·δ_n)(q·δ_n) - (A^(2)·δ_n)(q·δ_n) + iΔ_n^(2)(q·δ_n)_h_2,b].




By taking into account equation (<ref>) and ∑_n=1^3e^iK_0·δ_n=0, the q-independent terms in the last expression are cancelled. Thus, h can be rewritten as

    h=h_0+h_1,a+h_1,b+h_2,a,h_2,b,

where

    h_0= - t_0∑_n=1^3 e^iK_0·δ_n[iq·δ_n]=3t_0a_0/2(q_x+iq_y),


    h_1,a   =    - t_0∑_n=1^3 e^iK_0·δ_n[ iq·ϵ̅·δ_n]
       =    - t_0∑_n=1^3 e^iK_0·δ_n[ iQ·δ_n]
       =   3t_0a_0/2(Q_x+iQ_y)
       =   3t_0a_0/2[ϵ_xxq_x + ϵ_xyq_y + i(ϵ_xyq_x + ϵ_yyq_y)],


    h_1,b   =    - t_0∑_n=1^3 e^iK_0·δ_n (q·δ_n) [- (A^(1)·δ_n) + iΔ_n^(1)]
       =   -3t_0a_0/2β[ϵ_xxq_x + ϵ_xyq_y + i(ϵ_xyq_x + ϵ_yyq_y)],



    h_2,a   =    - t_0∑_n=1^3 e^iK_0·δ_n (q·ϵ̅·δ_n) [  -(A^(1)·δ_n) + iΔ_n^(1)]
       =    - t_0∑_n=1^3 e^iK_0·δ_n (Q·δ_n)[  -(A^(1)·δ_n) + iΔ_n^(1)]
       =   -3t_0a_0/2β[ϵ_xxQ_x + ϵ_xyQ_y + i(ϵ_xyQ_x + ϵ_yyQ_y)]
       =    -3t_0a_0/2β{(ϵ_xx^2+ϵ_xy^2)q_x + ϵ_xy(ϵ_xx+ϵ_xy)q_y
           + i[ϵ_xy(ϵ_xx+ϵ_xy)q_x + (ϵ_xx^2+ϵ_xy^2)q_y]}

and



    h_2,b   =    - t_0∑_n=1^3 e^iK_0·δ_n (q·δ_n) [ - i(A^(1)·δ_n)^2/2 - Δ_n^(1)(A^(1)·δ_n) - (A^(2)·δ_n) + iΔ_n^(2)]
       =    -3t_0a_0/2β1/8{ (ϵ_xx-ϵ_yy)^2q_x + 2ϵ_xy(ϵ_xx-ϵ_yy)q_y + i [2ϵ_xy(ϵ_xx-ϵ_yy)q_x - 4ϵ_xy^2q_y] }
           + 3t_0a_0/2β^21/4{ (ϵ_xx^2-ϵ_yy^2+2ϵ_xxϵ_yy+2ϵ_xy^2)q_x +  4ϵ_xxϵ_xyq_y + i [4ϵ_xxϵ_xyq_x + 2(ϵ_yy^2-ϵ_xy^2)q_y] },


being Q=q·ϵ̅. To simplify each term of h in equation (<ref>), we have used of K_0=(4π/3√(3)a_0,0), equations (2) of the main text for δ_n and the expressions obtained in the previous section for A^(1), A^(2), Δ^(1) and Δ^(2). In addition, note the same algebraic form between the initial expression of equation (<ref>) and equation(<ref>) if one defines  Q=q·ϵ̅. This similarity is also observed between equations (<ref>) and (<ref>).




In consequence, using equations (<ref>–<ref>) we obtain the contribution of each term of h to equation (<ref>) as

    (
    [     0 h_0^∗;   h_0     0 ])=
    ħ v_0τ·q,


    (
    [       0 h_1,a^∗;   h_1,a       0 ])=ħ v_0τ·ϵ̅·q,




    (
    [       0 h_1,b^∗;   h_1,b       0 ])=-ħ v_0βτ·ϵ̅·q,


    (
    [       0 h_2,a^∗;   h_2,a       0 ])=-ħ v_0βτ·ϵ̅^2·q,



and



    (
    [       0 h_2,b^∗;   h_2,b       0 ])=ħ v_0τ·(- βκ̅_1 + β^2κ̅_2)·q,

where v_0=3t_0a_0/2ħ is the Fermi velocity for pristine graphene, τ=(τ_x,τ_y) is a vector of (2×2) Pauli matrices, κ̅_1 and κ̅_2 are respectively the matrices (10) and (11) of the main text. To obtain the expressions (<ref>–<ref>) in terms of Pauli matrices, we used the identity


    ([                                    0 χ_xxq_x+χ_xyq_y - i[χ_xyq_x+χ_yyq_y]; χ_xxq_x+χ_xyq_y + i[χ_xyq_x+χ_yyq_y]                                    0 ])=
    τ·χ̅·q,


where χ_ij are the elements of an arbitrary (2×2) matrix χ̅, being χ̅=I, χ̅=ϵ̅,  χ̅=-βϵ̅,  χ̅=-βϵ̅^2 and χ̅=-βκ̅_1 + β^2κ̅_2 for equations (<ref>) to (<ref>), respectively. It is worth mentioning that expression (<ref>) is the effective Dirac Hamiltonian for pristine graphene, while (<ref>) and (<ref>) are corrections to first-order in ϵ̅, previously derived in Ref. [Oliva13]. The second-order corrections in ϵ̅, equations (<ref>) and (<ref>), are among the principal contributions of our work.

Finally, combining equations (<ref>), (<ref>) and (<ref>–<ref>), we obtain the effective Dirac Hamiltonian for graphene under a uniaxial strain, up to second-order in the strain tensor ϵ̅, given by

    H=ħ v_0τ·(I̅ + ϵ̅ - βϵ̅ -βϵ̅^2 - βκ̅_1 + β^2κ̅_2)·q,

which is the equation (<ref>) reported in the main text.




§ OPTICAL CONDUCTIVITY OF AN ANISOTROPIC DIRAC SYSTEM


In this section, we derive the optical conductivity tensor σ̅_ij(ω) of an anisotropic Dirac system described by the effective Hamiltonian

    H=ħ
    v_0τ·(I̅+Δ̅)·q,

where the anisotropic behaviour is expressed through the perturbation Δ̅, which is a symmetric (2×2) matrix such that Δ̅_ij≪1. Essentially, we now extend, up to second-order in Δ̅, a previous calculation of σ̅_ij(ω) up to first-order in Δ̅ reported in Ref. [Oliva2016].

Assuming that the considered system has linear response to an external electric
field of frequency ω, its optical conductivity σ̅_ij(w) can be calculated by combining the Hamiltonian (<ref>) and the Kubo formula. Following the approach used in
Refs. [Ziegler06,Ziegler07], σ̅_ij(ω) can be expressed as a double integral with respect to two energies E, E':

    σ̅_ij(ω)    =    ie^2/ħ∫∫Tr{ v_iδ(H-E')v_jδ(H-E)}
          ×1/E-E'+ħω - iαf(E)-f(E')/E-E'd
    Ed E',

where f(E)=(1+exp[E/(k_BT)])^-1 is the Fermi function at
temperature T,  is the trace operator including the summation over the q-space (as defined in equation (7) of Ref. [Ziegler07]) and v_l=i [H,r_l] is the velocity operator in the l-direction, with l=x,y.

To calculate the integral (<ref>) it is convenient to make the change of variables

    q=(I̅+Δ̅)^-1·q^*,

which yields that the Hamiltonian (<ref>) becomes
H=ħ v_0τ·q^*, corresponding to the case of a unperturbed and isotropic Dirac system, as unstrained graphene. At the same time, the velocity operator components transform as

    v_x   =   i [H,r_x]=∂ H/∂ q_x, 
       =   (∂ H/∂ q_x^*∂
    q_x^*/∂ q_x + ∂ H/∂ q_y^*∂ q_y^*/∂ q_x), 
       =   (1 + Δ̅_xx)v_x^*
    + Δ̅_xyv_y^*,

and analogously

    v_y=(1 + Δ̅_yy)v_y^*
    + Δ̅_xyv_x^*,

where v_x^*=(∂ H/∂ q_x^*) and
v_y^*=(∂ H/∂ q_y^*) are the velocity
operator components for the unperturbed Dirac system.

Then, substituting equations (<ref>) and (<ref>) into
equation (<ref>) we find

    σ̅_xx(ω)   =   [(1+Δ̅_xx)^2+Δ̅_xy^2]Jσ_0(ω),
    σ̅_yy(ω)   =   [(1+Δ̅_yy)^2+Δ̅_xy^2]Jσ_0(ω),

and

    σ̅_xy(ω)   =   σ̅_yx(ω)
       =   [2Δ̅_xy+Δ̅_xy(Δ̅_xx+Δ̅_yy)]Jσ_0(ω),

where J is the Jacobian determinant of the transformation (<ref>) originated by expressing the trace operator  of equation (<ref>) in the new variables q^* and σ_0(ω)
is the optical conductivity of the unperturbed Dirac system, i.e., the reported optical conductivity of unstrained graphene <cit.>. Note that equations (<ref>-<ref>) can be written in a compact manner as

    σ̅(ω)=(I̅ + 2Δ̅ + Δ̅^2)Jσ_0(ω).


Now, if J is expressed up to second-order in Δ̅, results

    J   =   [(I̅+Δ̅)^-1]
    ≈(I̅-Δ̅+Δ̅^2)
       ≈    1-(Δ̅) + [(Δ̅)]^2-(Δ̅)
       ≈    1-(Δ̅) + [(Δ̅)]^2/2 + (Δ̅^2)/2,

where (Δ̅)=Δ̅_xx+Δ̅_yy.



Finally, substituting equation (<ref>) into (<ref>) we obtain that the optical conductivity tensor of the anisotropic Dirac system, described by the Hamiltonian (<ref>), is given by

    σ̅(ω)    ≈   σ_0(ω){I̅ - (Δ̅)I̅ + 2Δ̅ + Δ̅^2
       +   1/2[((Δ̅))^2
    + (Δ̅^2)]I̅-2(Δ̅)Δ̅},

where the second term is the contribution to second-order in the perturbation Δ̅ to the conductivity.









