Anomalous Nernst effect in strained graphene coupled to a substrate inducing a time-reversal symmetry breaking

The usual Nernst effect is a thermomagnetic phenomenon of a triggered electric current transverse to a temperature gradient and an external magnetic field applied perpendicular to both the temperature gradient and the charge current [1]. This effect is currently discussed as a sensitive probe in the context of high-T_c superconductivity [2] and also for graphene-based systems, e.g. monolayer graphene [3, 4] and graphite [5]. Owing to the Dirac fermion character in graphene, the Nernst coefficient deviates markedly from the Mott relation for the n = 0 Landau level (LL) around the Dirac point [3, 4]. Theoretically this deviation is yet to be fully understood [3]. On the other hand, the Nernst effect goes along with the quantum Hall effect which is inherent for a topologically protected edge current. The number of edge states is just the Chern–Simon number (or the winding number) which is given by an integral over the Berry curvature on the momentum torus for the two-dimensional (2D) Brillouin zone (BZ) with periodic boundary. The difference between the Hall conductance and the Nernst coefficient is that the latter is determined not only by the Berry curvature but also by the entropy generation around the Fermi surface. Therefore, the Nernst effect is sensitive to the change of the Fermi energy, the temperature and other quantities. The common feature of these two effects is that in an external electric field the charge carrier attains an extra anomalous velocity with ℏv_k = eE × Ω_νK(K')z, where E is the electric field, Ω_νK(K')z is the Berry curvature, ν = 1 (2) stands for the conduction (valence) band and K (K') indicates the Dirac cone located at K (K') (see figure 1(b)). Hence, a finite Berry curvature may lead to an anomalous Nernst effect in the absence of an external magnetic field which is the case discussed here.

Zoom In Zoom Out Reset image size

For gapless graphene [6], the Berry curvature shows a singularity at the Dirac points so that the Berry phases ±π are associated with the two Dirac cones. A substrate-induced gap (hereafter referred as the substrate gap) was observed experimentally when a SiC substrate is used [7] and in an ab initio calculation with a BN substrate [8]. For the graphene layer, such a substrate gap yields a finite Berry curvature for each cone. However, these Berry curvatures have opposite signs [9, 10]. In a recent study, we showed that in a time-reversal symmetry (TRS)-breaking environment, strain can lead to an anomalous Hall effect [11].

Since both the quantum Hall effect and the Nernst effect are intrinsically linked to the topology of the momentum space, it is interesting to investigate how to engineer the topology of the Dirac cones of graphene. In this respect, the strain is currently intensely studied [12]. It is expected, however, that the strain-induced static electric potential does not break the TRS. To study whether the TRS is still maintained in an essentially transient non-equilibrium transport (e.g. Nernst effect when strain builds up) is out of the scope of this paper [13]⁴. Here, we consider a specific model that provides an explicit TRS breaking. We study a graphene layer on a 'Haldane' substrate under a uniaxial strain. The substrate effect is described within a Haldane-type model [14, 15], i.e. a zero flux through one unit cell but with a vector potential structure. The Haldane model takes the next-nearest-neighbor hopping into account and leads to massive Dirac fermions. The nearest-neighbor hopping leads to coupling to the planar motion of Dirac fermions. In this work, we will show that a novel anomalous Nernst effect emerges and can be engineered by the strain. In contrast to the case without strain [10], the total Nernst coefficient (TNC) does not vanish and undergoes a sign change when increasing an applied gate voltage (for a small flux). We found a robust six-leaf geometry of the TNC as a function of the angle of the strain with respect to the zigzag direction.

The Pauli matrices σ₁ = σ_x and σ₂ = σ_y represent the E representation of C_3v group at the Dirac points of graphene and are related to the A and B sublattice character. σ_x and σ_y indicate that there is a coupling between A and B sublattices which can be in the lowest order provided by the nearest-neighbor hopping, while σ₃ = σ_z essentially describes the breaking of this sublattice degeneracy and gives rise to mass generation of the Dirac fermion. This mass can be generated by a stagger field due to the static electric interaction with the substrate. For example, if the substrate has a honeycomb structure (or from a top view is felt by the Dirac fermions in the graphene layer), but the two atoms in the unit cell are not equivalent, then a staggered field acts on the two carbon atoms in the unit cell of the graphene. This coupling generates a mass of the Dirac fermion and enters into the coupling to σ₃. Another contribution entering the coupling to σ₃ is the TRS breaking term in the Haldane model. The Haldane TRS breaking term is a consequence of an extra phase due to the next-nearest-neighbor hopping. The low-energy effective Hamiltonian is written as H = σ·d(k). From the theory of invariants [16], we obtain the d vector in the K and K' cones without including explicitly the external electric field

$$\begin{eqnarray} &&\begin{array}{llll} \displaystyle d_{0} = t_{2}\bar{\eta}\cos\phi, \\ \displaystyle d_{1} = (\tilde{v}_{F}k_{x}+\tilde{\mathbf{v}}\cdot\mathbf{k})\tau_{z}-\mathcal{A}_{x}, \\ \displaystyle d_{2} = \tilde{v}_{F}k_{y}+(\mathbf{k}\times\tilde{\mathbf{v}})_{z}+\mathcal{A}_{y}\tau_{z},\\ \displaystyle d_{3} = \frac{\Delta}{2}-\sqrt{3}\tau_{z}t_{2}\bar{\eta}\sin\phi, \end{array} \end{eqnarray} \tag{ 1 }$$

where $\vartriangle$ is the energy gap induced by the stagger field of the substrate (substrate gap), τ_z denotes the valley degree, $\tilde {v}_{F}=v_{F}+b_{62}(u_{xx}+u_{yy})$, $\tilde {v}_{x}=b_{63}(u_{xx}-u_{yy})$, $\tilde {v}_{y}=2b_{63}u_{xy}$, $\bar {\eta }=3+\sum _{i}\frac {\delta t_{2i}}{t_{2}}$ is a dimensionless function varying with strain, δt_2i is the variation induced by strain to the next-nearest-neighbor hopping labeled by i = 1,2,3 and ϕ is the flux of the effective magnetic field generated by the substrate (see the appendices). We should point out that the d₀ and d₃ terms are derived by the next-nearest-neighbor hopping (hopping in the same sublattice) and the d₁ and d₂ terms are obtained from the nearest-neighbor hopping (hopping between different sublattices). Furthermore, we introduced the strain tensor as defined in the continuous elastic theory as $u_{ij}=\frac {1}{2}(\frac {\partial u_{i}}{\partial r_{j}}+\frac {\partial u_{j}}{\partial r_{i}}+\frac {\partial u_{z}}{\partial r_{i}}\frac {\partial u_{z}}{\partial r_{j}})$, where u is a polar vector indicating the displacements of atoms. In principle, all the b_ij parameters can be determined experimentally or from ab initio calculations [16] in a similar manner as done for semiconductors. Another way is to link these parameters to those derived from tight-binding formalism. This is obviously a transparent and effective approach to connect with microscopic method. However, this way is not explored thoroughly for graphene. For example, the case to study the deformed next-nearest-neighbor hopping is still missing, to our knowledge. For the effect of the deformed nearest-neighbor hopping, research is being undertaken extensively. If only the deformed nearest-neighbor hopping integrals are taken into account [17, 18], one finds an effective gauge field that is independent of the momentum. Recently, it was pointed out that the strain influence should be envisaged by including the effects of the deformed lattice in the real space, the deformed momentum vectors and also the deformed hopping integrals [19]. After expanding the Hamiltonian around the Dirac points, extra terms can be obtained which may be momentum-dependent. This changes definitely the topology of the Dirac cones. A corresponding effect is, however, still unexplored for the next-nearest-neighbor case. To that end and for the purpose of our study, we worked out an effective Hamiltonian in tight-binding formalism for this case and obtained a concrete model with a TRS breaking under strain in appendix B.

Thus, every term in equation (1) can be determined by tight-binding parameters and finally determined by solely the strain tensor. In this way, the model is well founded and its parameters are linked in a transparent way to microscopic quantities. $\mathcal {A}_{x}=b_{61}(u_{xx}-u_{yy})$ and $\mathcal {A}_{y}=b_{61}(2u_{xy})$ are to be viewed as effective vector potentials in a minimal coupling scheme for the strain-induced transport [17, 18]. It is instructive to recall that the vector potential is expressible in terms of the hopping parameters in a tight-binding formalism $\mathcal {A}_{x}=\frac {1}{2}\left (\triangle t_{02}+\triangle t_{03}-2\triangle t_{01}\right )$ and $\mathcal {A}_{y}=\frac {\sqrt {3}}{2}\left (\triangle t_{02}-\triangle t_{03}\right )$, where $\triangle t_{0i}$, i = 1,2,3, are the variation of the nearest-neighbor hopping integrals induced by the strain with respect to the change of the nearest-neighbor vectors δ_i shown in appendix D (a similar expression was given in [20]). Thus, the essential feature is that vector d mediates a map from a 2D momentum space to a three-dimensional vector space. We will study the consequences of this mapping beyond the minimal coupling.

Without strain, the cross section of the Fermi level and the Dirac cone is a circle (for E_f is not exactly at the Dirac point). Upon including the strain, this circular topology of the cones is changed into an ellipsoidal ring. Gaps are opened at the original Dirac points at K and K'. However the positions of the minimal energy gap are no longer located at these original Dirac points, they are floating elsewhere when changing the strain since the strain serves in effect as a magnetic field in the momentum space and the Dirac points can be regarded as monopoles in the same space. To capture this behavior, we use the linearized momentum with respect to the original Dirac points. In the absence of strain, the Hamiltonian reduces to the linear Hamiltonian around these points. To study the Nernst effect carefully, we have to use all the components in the first BZ. Thus, we use a rhombic BZ marked by the shaded regions in figure 1(b). In this way, the number of states in the BZ is conserved after linearization of the spectrum around the Dirac points, which is as an alternative to the usual Debye prescription [21].

To be specific, we consider throughout the work a uniaxial strain only. The x direction is selected as the zigzag direction; the momentum space x direction (k_x) is shown in figure 1. The angle of the applied tension with respect to x direction is θ. We first calculate the energy contour for the conduction band in the presence of the strain, as shown in figures 2(a) and (b). First of all, we observe that the circular Dirac cone is deformed to an elliptical cone for both K and K' cones. The original Dirac points are shifted away by the strain since the strain effect serves partly as a gauge field in the momentum space. For zero gap graphene, the Dirac points are monopoles in the BZ. Although there is a gap considered here, they should be affected by the induced gauge field. Although the TRS breaking in the model studied here is based on the Haldane model, it is instructive to think about some other situations with a possible TRS breaking. The TRS breaking environment could be created by applying a weak magnetic field in the graphene sheet so that the orbital motion and the Zeeman effect can be ignored. Graphene easily wrinkles into the third dimension with nanometer-scale ripples that have been observed in microscopic studies [22]. The existence of the rippling has been demonstrated to be not only extrinsic, e.g. by surface roughness of the substrate, but also intrinsic in graphene [23]. Therefore, an in-plane magnetic field induces additionally a random magnetic field both in the xy direction and z direction [24–26]. TRS is expected not to be present at all. Alternatively, the graphene layer may be put on a magnetic substrate to break the TRS. In this case, a non-vanishing effect enters into the third component of the d vector. This is reflected by the asymmetric K and K' in figures 2(a) and (b) since different gaps are opened in addition to the uniform gap $\vartriangle$. It is useful to recall that another mechanism to break the TRS might be having spontaneous microscopic currents due to the order parameter of a density wave as that in the context of underdoped cuprate superconductor where a square lattice is considered [27, 28].

Zoom In Zoom Out Reset image size

What we intend to study is an anomalous Nernst effect that emerges for a non-vanishing Berry curvature. This Berry curvature manifests itself in an additional anomalous velocity term in the presence of an external electric field ℏv_k = eE × Ω. Thus, the velocity multiplied by the entropy density gives rise to the Nernst coefficient [10, 27]

$$\begin{equation} \alpha^{n}_{xy}=4\pi\alpha_{0}\sum_{\nu}\int\frac{\mathrm{d}^{2}\mathbf{k}}{(2\pi)^{2}}\Omega_{\nu nz}(\mathbf{k})S_{\nu}(\mathbf{k}), \end{equation} \tag{ 2 }$$

where α₀ = ek_B/h, e is the electron charge, k_B is the Boltzmann constant, ℏ is the Planck constant and n = K or K', ν = + 1 (−1) for the conduction (valence) band, S_ν(k) = −f_νk ln f_νk − (1 − f_νk) ln (1 − f_νk) is the entropy density of the ν band and f_νk = f[E_ν(k)] is the Fermi distribution function. Please note that equation (2) yields the Nernst coefficient for one cone. The integral over k should be taken in the respective triangle BZ shown in figure 1(b).

In figures 2(c) and (d), the entropy density is calculated for the K and K' cones for E_f = 1.2 eV. Different from the energy contour, the entropy density displays a maximum around the Fermi energy. The ellipsoidal topology is clearly seen from the maximum of the entropy density, as well. Another observation is that TRS breaking is observed in the symmetry breaking of the two Dirac cones. This can be inferred from the values of the entropy density and the different positions of the shifted Dirac points in the momentum space.

It has been shown in [11] that the Berry curvatures for the two cones have the same signs rather than opposite signs in the absence of the gap Δ, but in a TRS breaking environment. Therefore, the anomalous transverse velocities will be the same for the two cones leading to a net Hall current driven by an external electric field. It is natural to expect that an anomalous Nernst coefficient will be associated with this anomalous transverse velocity. In figures 3(a) and (b), TNC which is the sum of the contributions from the two Dirac cones is calculated while varying the Fermi level for the different magnetic fluxes and the strain magnitudes in figures 3(a) and (b), respectively.

Zoom In Zoom Out Reset image size

The peak in the TNC when the Fermi level lies in the conduction band can be understood as follows. The Nernst coefficient is affected by two factors. (a) The Berry curvature which exhibits a maximum value at the Dirac points or the bottom of the conduction band. (b) The maximum of the entropy density is, however, around the Fermi level leading to a finite-width Fermi ring in the momentum space. Thus, when the Fermi level increases, the area of the entropy density in the momentum space increases; while the Berry curvature decreases. The overall effect manifests itself as a peak of the Nernst coefficient. The peak shift induced by the flux and the strain is a consequence of increasing the asymmetry of the energy and the gaps in the two cones (see figure 2(a)). It is interesting to note that the magnitude of the peak of the TNC is first enhanced slightly by increasing the strain from a small distortion and then is decreased by a strong strain (see figure 3(b)).

When the gap Δ is present, the behavior of the Nernst coefficient is determined by the relative magnitudes of this gap and the gap induced by the TRS breaking which is affected by the strain. Hereafter we call the gap Δ the stagger field gap (SFG) and the latter the time-reversal-breaking (TRB) gap. When SFG is dominant, the Berry curvatures for the two Dirac cones are of the opposite signs. Since the total gaps at the two Dirac points are asymmetric, one cone yields at first a dominant contribution to the Nernst signal when the Fermi level is touching the band edge of this cone. When further increasing the Fermi level, the other cone will affect the Nernst signal by reducing it due to the opposite signs of the Berry curvatures of the two cones. This explains the marked difference and even the sign change shown in figure 3(c). However, when the TRB gap is enlarged with larger flux, the Berry curvature of one cone (K) is gradually converted by its sign. This makes the negative TNC turn to a positive TNC at K point. When the TRB gap is dominant, the signs of the Berry curvatures for the two cones are the same, leading to a two-peak structure of the TNC varying with the Fermi level. In figure 3(d), the peaks of the TNC are shifted to a larger Fermi energy with increasing the strain in accordance with the enlarged total gap.

We shall compare the behavior of the Nernst signal in the whole gate voltage range, scanning the case when the Fermi level both in the conduction band and the valence band. Since both SFG and TRB gaps do not change the particle–hole nature of the conduction and valence bands, the Berry curvatures take opposite signs for these two bands. Therefore, the sign of the Nernst coefficient is anti-symmetric with respect to the neutral Dirac point. This is in complete contrast to the case when applying a magnetic field. In the presence of an external magnetic field, Landau quantized orbits are formed. The Nernst signal is [3]

$$\begin{equation} \alpha_{xy}(E_{\mathrm{f}})=\frac{k_{\mathrm{B}}e}{h}\beta\sum_{n}\int_{E_{n}}^{\infty}\mathrm{d}\epsilon(\epsilon-E_{\mathrm{f}})\left(-\frac{\partial f}{\partial\epsilon}\right), \end{equation} \tag{ 3 }$$

where β = (k_BT)⁻¹ and E_n indicates the energy of nth LL. The Nernst signal shows peaks near the LLs, while it tends to zero in between the LLs [29, 30]. This gives rise to a symmetric profile of the Nernst signal when crossing the Dirac point due to the symmetric distribution of the LLs with respect to the Dirac point. This distinctive symmetry difference indicates the different nature of the case studied here to the physics under an external magnetic field.

In figure 4, the variation of the TNC with the strain angle is shown for zero Δ, finite Δ and different strain cases. The graph exhibits a six-leaf geometry as we vary θ. The maxima are all located at θ = 0°, 60°, 120°, 180°, 240° and 300° (for the selected E_f = 0.7 eV). In between the maxima the Nernst coefficient turns six minima at θ = 30°, 90°, 150°, 210°, 270° and 330° (a few directions, i.e. θ = π/6, π/3 and π/2, are shown explicitly in figure 1). The positions of the maxima and minima of the Nernst coefficients are robust against varying the parameters of the system, such as ε, the SFG Δ and the Fermi levels. It is interesting to note that the directions where the Nernst coefficients attain maxima pass through the two Dirac points, while the TNC is minimal when the strain is along the primitive vectors in the momentum space and σ_d1, σ_d2 and σ_d3 in the C_3v group for the Dirac cones [31]. Estimating the order of the obtained TNC we note that α₀ ≈ 3.33 nA K⁻¹ corresponds to 100 μV K⁻¹ of a quantum resistivity $\frac {h}{e^{2}}$. Thus, 10⁻²α₀ corresponds to 1000 nV K⁻¹, which might be measurable.

Zoom In Zoom Out Reset image size

We demonstrated that a novel anomalous Nernst effect driven by an external electric field rather a magnetic field can be induced in strained graphene coupled to a TRB environment. The Haldane model under strain is studied as a specific illustration. We investigated the action of the Berry curvatures of the two Dirac cones and how it is changed by strain. The interplay of the stagger-field-induced gap of the substrate and the TRS-breaking-induced gap is studied. The former preserves the opposite signs of the Berry curvatures, whereas the latter renders the same signs of the Berry curvatures. For weak TRS breaking, the total Nernst signal experiences a sign change when the Fermi level is elevated. When the TRB gap is large, TNC shows a two-peak structure with varying Fermi level. This sign change does not occur if only the SFG is missing. Hence, this effect may serve as a thermal-electric probe for the existence of such a gap. We found a robust six-leaf configuration of the TNC as the angle of the strain varies with respect to the zigzag direction.

We thank C L Jia, K-H Ding and V K Dugaev for useful discussions. JB was supported by DFG. ZGZ was supported by DFG-SFB 762, and is now supported by Hundred Talents Program of The Chinese Academy of Sciences.

The Hamiltonian including the second nearest-neighbor hopping reads

$$\begin{equation} H_{2}=-\sum_{\langle\langle i,j\rangle\rangle}\left(t^{a}_{ij}\mathrm{e}^{+\mathrm{}\frac{e}{\hbar}\int_{R_{a_{i}}}^{R_{a_{j}}}\mathbf{A}\cdot \mathrm{d}\mathbf{r}}a^{\dagger}_{i}a_{j} +t^{b}_{ij}\mathrm{e}^{+\mathrm{}\frac{e}{\hbar}\int_{R_{b_{i}}}^{R_{b_{j}}}\mathbf{A}\cdot \mathrm{d}\mathbf{r}}b^{\dagger}_{i}b_{j}+\textrm{H.c.}\right), \end{equation} \tag{ A.1 }$$

A is the vector potential in the Haldane model, the integrals in the exponents are to be performed along the paths connected by the atoms located at R^a(b)_i and R^a(b)_j for A (B) sublattice. A method to construct a Wigner–Seitz (WS) unit cell is shown in figure A.1(a) in which a hollow position is selected as the starting ends of primitive vectors. Therefore, each atom among the six atoms located at the corners of the honeycomb belongs to three WS unit cells and there are two atoms effectively in one WS unit cell. The vectors δ_i are the three vectors connecting the nearest neighbors. a₁ = δ₁ − δ₃ and a₂ = δ₁ − δ₂ are two primitive vectors; in addition a₃ = δ₃ − δ₂. Since the hollow position is labeled by R_i for the ith WS unit cell, the locations of the six atoms are

$$\begin{eqnarray} &&\begin{array}{llllll} \displaystyle \mathbf{A}_{1}=\mathbf{R}_{i}-\boldsymbol{\delta}_{1}, \\ \displaystyle \mathbf{A}_{2}=\mathbf{R}_{i}-\boldsymbol{\delta}_{3}, \\ \displaystyle \mathbf{A}_{3}=\mathbf{R}_{i}-\boldsymbol{\delta}_{2}, \\ \displaystyle \mathbf{B}_{1}=\mathbf{R}_{i}+\boldsymbol{\delta}_{3}, \\ \displaystyle \mathbf{B}_{2}=\mathbf{R}_{i}+\boldsymbol{\delta}_{2}, \\ \displaystyle \mathbf{B}_{3}=\mathbf{R}_{i}+\boldsymbol{\delta}_{1}. \end{array} \end{eqnarray} \tag{ A.2 }$$

Zoom In Zoom Out Reset image size

If we use a vector to indicate the hopping between the carbon atoms (for example ta^†b), the reversed vector means a Hermitian conjugate hopping (for example tb^†a). For the nearest-neighbor hopping, we only need to write the terms constituting a clockwise loop in one WS unit cell, and the Hermitian conjugate terms will be recovered by doing the same thing in another WS unit cell. From figure A.1(b), we clearly see this and there are two inverse vectors for each nearest 'bond'. For the next-nearest hopping, we need to include both the clockwise and anti-clockwise loops along the two triangles since the next-nearest 'bonds' are not the common edges of the WS unit cells (see figure A.1(c)).

For each sublattice, we can group three terms out of the six terms, and the three terms left are then the Hermitian conjugate of the selected three ones. We select A₃ → A₂ → A₁ → A₃ and B₁ → B₃ → B₂ → B₁. In the presence of the vector potential, we have $({\rm sum})_{1}=\int _{A_{1}}^{A_{3}}\mathbf{A}\cdot \mathrm{d}\mathbf{r}+\int_{A_{3}}^{A_{2}}\mathbf{A}\cdot \mathrm{d}\mathbf{r}+\int _{A_{2}}^{A_{1}}\mathbf {A}\cdot \mathrm {d}\mathbf {r}= -\oint _{\circlearrowleft }\mathbf {A}\cdot \mathrm {d}\mathbf {r}=-\int\int _{\vartriangle A_{1}A_{2}A_{3}}\mathbf {B}\cdot \mathrm {d}\mathbf {r}=\int\int _{S-\vartriangle A_{1}A_{2}A_{3}}\mathbf {B}\cdot \mathrm {d}\mathbf{r}$, where S is the area of the honeycomb, $\vartriangle A_{1}A_{2}A_{3}$ is the triangle formed by three points located at A₁, A₂ and A₃. The zero flux through the honeycomb has been used for the last equation. Since the field B(r) is supposed to have full symmetry of the crystal, we have $\frac {e}{\hbar }(\textrm {sum})_{1}=2\pi (\textrm {sum})_{1}/\phi _{0}=\frac {2\pi }{\phi _{0}}\int\int _{S-\vartriangle A_{1}A_{2}A_{3}}\mathbf {B}\cdot \mathrm {d}\mathbf {r}=3\phi$, where $\phi =\frac {2\pi }{\phi _{0}}\left [2\Phi _{a}+\Phi _{b}\right ]$, Φ_a(b) is the flux through the region labeled by a and b in figure A.1(a).

The items for A sublattice in the nearest-neighbor hopping are now

$$\begin{eqnarray} &&\fl ({a}) =-\sum_{i}\left\{\mathrm{e}^{\mathrm{i}\phi}\left(t^{a}_{13}a^{\dagger}_{1}a_{3}+t^{a}_{32}a^{\dagger}_{3}a_{2}+t^{a}_{21}a^{\dagger}_{2}a_{1}\right)+\mathrm{H.c.}\right\} \nonumber\\ &&= -\sum_{i}\left\{\mathrm{e}^{\mathrm{i}\phi}t^{a}_{13}\frac{1}{N}\sum_{\mathbf{kk'}}\mathrm{e}^{-\mathrm{i}\mathbf{k}\cdot(\mathbf{R}_{i}-\boldsymbol{\delta}_{1})} \mathrm{e}^{\mathrm{i}\mathbf{k'}\cdot(\mathbf{R}_{i}-\boldsymbol{\delta}_{2})}a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}'}+(\textrm{term for } t^{a}_{32})+(\textrm{term for } t^{a}_{21})+\mathrm{H.c.}\right\}. \nonumber\\ \end{eqnarray} \tag{ A.3 }$$

We use a relation $\frac {1}{N}\sum _{i}\mathrm {e}^{\mathrm {i}\mathbf {R}_{i}\cdot (\mathbf {k}-\mathbf {k}')}=\delta _{\mathbf {k}\mathbf {k}'}$ in the above equation, leading to

$$\begin{equation} \fl ({a})=-2\sum_{\mathbf{k}}\left[\cos(\phi+\mathbf{k}\cdot\mathbf{a}_{2})t^{a}_{13}+\cos(\phi-\mathbf{k}\cdot\mathbf{a}_{3})t^{a}_{32}+ \cos(\phi-\mathbf{k}\cdot\mathbf{a}_{1})t^{a}_{21}\right]a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}. \end{equation} \tag{ A.4 }$$

If we define l₁ = a₁, l₂ = −a₂ and l₃ = a₃, we have

$$\begin{equation} \fl ({a})=-2\sum_{\mathbf{k}}\left[\cos(\phi-\mathbf{k}\cdot\mathbf{l}_{2})t^{a}_{13}+\cos(\phi-\mathbf{k}\cdot\mathbf{1}_{3})t^{a}_{32}+ \cos(\phi-\mathbf{k}\cdot\mathbf{l}_{1})t^{a}_{21}\right]a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}. \end{equation} \tag{ A.5 }$$

For graphene layer without strain, we set all the second nearest hopping to be the same t₂. We have then

$$\begin{equation} ({a})=-2t_{2}\sum_{\mathbf{k}}\left[\sum_{i=1}^{3}\cos(\phi-\mathbf{k}\cdot\mathbf{l}_{i})\right]a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}. \end{equation} \tag{ A.6 }$$

A similar calculation can be performed for the B sublattice giving rise to

$$\begin{eqnarray} &&\fl (b) = -2\sum_{\mathbf{k}}\left[\cos(\phi+\mathbf{k}\cdot\mathbf{a}_{1})t^{b}_{13}+\cos(\phi-\mathbf{k}\cdot\mathbf{a}_{2})t^{b}_{32}+ \cos(\phi+\mathbf{k}\cdot\mathbf{a}_{3})t^{b}_{21}\right]b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}\nonumber\\ &&= -2t_{2}\sum_{\mathbf{k}}\left[\sum_{i=1}^{3}\cos(\phi+\mathbf{k}\cdot\mathbf{l}_{i})\right]b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}. \end{eqnarray} \tag{ A.7 }$$

Thus, the Hamiltonian for the next-nearest-neighbor hopping is

$$\begin{eqnarray} &&\fl H_{2}=(a)+(b) = -2t_{2}\sum_{\mathbf{k}}\left\{\cos\phi\left(\sum_{i=1}^{3}\cos\mathbf{k}\cdot\mathbf{l}_{i}\right)\left(a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}+ b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}\right)\right. \nonumber\\ &&+\left.\sin\phi\left(\sum_{i=1}^{3}\sin\mathbf{k}\cdot\mathbf{l}_{i}\right)\left(a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}- b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}\right)\right\}. \end{eqnarray} \tag{ A.8 }$$

In these appendices, we use K_± to indicate K and K'. In the presence of strain, the next-nearest-neighbor hopping is changed, and the vectors l_i in real space are also deformed. We write explicitly t'_2i = t₂ + δt_2i, where i = 1, 2 and 3 indicate the three next-nearest-neighbor 'bonds'. And we have l'_i = l_i + δl_i. Under the strain, the Haldane model for the next-nearest-neighbor hopping is

$$\begin{eqnarray} &&\fl H_{2} = -2\sum_{\mathbf{k}}\left\{\! \cos\phi\left(\sum_{i=1}^{3}t'_{2i}\cos\mathbf{k}\cdot\mathbf{l}'_{i}\right)\left(a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}+ b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}\right)\!+\! \sin\phi\left(\sum_{i=1}^{3}t'_{2i}\sin\mathbf{k}\cdot\mathbf{l}'_{i}\right)\left(a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}- b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}\right)\!\right\}\!,\nonumber\\ \end{eqnarray} \tag{ B.1 }$$

$$\begin{eqnarray} &&\fl (\mathrm{I}) = \sum_{i=1}^{3}t'_{2i}\cos\mathbf{k}\cdot\mathbf{l}'_{i}\nonumber\\ &&=\sum_{i}(t_{2}+\delta t_{2i})\cos\left[(\mathbf{K}_{\pm}+\delta\mathbf{k})\cdot(\mathbf{l}_{i} +\delta\mathbf{l}_{i})\right] \nonumber\\ &&\approx \sum_{i}(t_{2}+\delta t_{2i})\left\{\cos(\mathbf{K}_{\pm}\cdot\mathbf{l}_{i})-\sin(\mathbf{K}_{\pm}\cdot\mathbf{l}_{i}) \left(\mathbf{K}_{\pm}\cdot\delta\mathbf{l}_{i}+\delta\mathbf{k}\cdot\mathbf{l}_{i}\right)\right\} \nonumber\\ &&\approx \sum_{i}\left\{t_{2}\cos(\mathbf{K}_{\pm}\cdot\mathbf{l}_{i}) -t_{2}\sin(\mathbf{K}_{\pm}\cdot\mathbf{l}_{i})\left(\mathbf{K}_{\pm}\cdot\delta\mathbf{l}_{i}+\delta\mathbf{k}\cdot\mathbf{l}_{i}\right) +\delta t_{2i}\cos(\mathbf{K}_{\pm}\cdot\mathbf{l}_{i}) \right\}.\nonumber\\ \end{eqnarray} \tag{ B.2 }$$

To proceed further we list

$$\begin{eqnarray} \mathbf{l}_{1} &=& a\left(\frac{1}{2},\frac{\sqrt{3}}{2}\right), \nonumber\\ \mathbf{1}_{2} &=& a\left(\frac{1}{2},-\frac{\sqrt{3}}{2}\right), \\ \mathbf{1}_{3} &=& a\left(-1,0\right),\nonumber \end{eqnarray} \tag{ B.3 }$$

where $a=\sqrt {3}a_{\mathrm {cc}}$, a_cc is the bond length for carbon–carbon atoms. By using $\mathbf {K}_{\pm }=\frac {2\pi }{a}\left (\pm \frac {2}{3},0\right )$, we have

$$\begin{eqnarray} \mathbf{K}_{\pm}\cdot\mathbf{l}_{1} &=& \pm\frac{2\pi}{3}, \nonumber\\ \mathbf{K}_{\pm}\cdot\mathbf{l}_{2} &=& \pm\frac{2\pi}{3}, \\ \mathbf{K}_{\pm}\cdot\mathbf{l}_{3} &=& \mp\frac{4\pi}{3}.\nonumber \end{eqnarray} \tag{ B.4 }$$

Therefore, we obtain

$$\begin{equation} \sum_{i}\cos\mathbf{K}_{\pm}\cdot\mathbf{l}_{i}=-3/2, \end{equation} \tag{ B.5 }$$

$$\begin{equation} \sum_{i}\delta t_{2i}\cos\mathbf{K}_{\pm}\cdot\mathbf{l}_{i}=-\frac{1}{2}(\delta t_{21}+\delta t_{22}+\delta t_{23}), \end{equation} \tag{ B.6 }$$

$$\begin{equation} -t_{2}\sum_{i}\sin(\mathbf{K}_{\pm}\cdot\mathbf{l}_{i})\mathbf{K}_{\pm}\cdot\delta\mathbf{l}_{i}=\mp t_{2}\mathbf{K}_{\pm} \cdot(\delta\mathbf{l}_{1}+\delta\mathbf{l}_{2}+\delta\mathbf{l}_{3}), \end{equation} \tag{ B.7 }$$

$$\begin{equation} -t_{2}\sum_{i}\sin(\mathbf{K}_{\pm}\cdot\mathbf{l}_{i})\delta\mathbf{k}\cdot\mathbf{l}_{i}=\mp t_{2}\delta\mathbf{k} \cdot(\mathbf{l}_{1}+\mathbf{l}_{2}+\mathbf{l}_{3}). \end{equation} \tag{ B.8 }$$

Therefore, we find

$$\begin{equation} (\mathrm{I})=-\frac{3}{2}t_{2}-\frac{1}{2}\sum_{i}\delta t_{2i}\mp t_{2}\left[\mathbf{K}_{\pm}\cdot\sum_{i}\delta\mathbf{l}_{i} +\delta\mathbf{k}\cdot\sum_{i}\mathbf{l}_{i}\right]. \end{equation} \tag{ B.9 }$$

Similar calculations yield

$$\begin{equation} \fl (\mathrm{II})=\sum_{i}t'_{2i}\sin(\mathbf{k}\cdot\mathbf{l}'_{i})=-\frac{3\sqrt{3}}{2}t_{2}\pm\frac{\sqrt{3}}{2}\sum_{i}\delta t_{2i}-\frac{1}{2} t_{2}\left[\mathbf{K}_{\pm}\cdot\sum_{i}\delta\mathbf{l}_{i} +\delta\mathbf{k}\cdot\sum_{i}\mathbf{l}_{i}\right]. \end{equation} \tag{ B.10 }$$

Since $\sum _{i}\mathbf {l}_{i}=0$, we obtain

$$\begin{eqnarray} &&\begin{array}{ll} \displaystyle (\mathrm{I}) = -\frac{3}{2}t_{2}-\frac{1}{2}\sum_{i}\delta t_{2i}\mp t_{2}\left(\mathbf{K}_{\pm}\cdot\sum_{i}\delta\mathbf{l}_{i}\right), \\ \displaystyle (\mathrm{II}) = \pm\frac{3\sqrt{3}}{2}t_{2}\pm\frac{\sqrt{3}}{2}\sum_{i}\delta t_{2i}-\frac{1}{2} t_{2}\left(\mathbf{K}_{\pm}\cdot\sum_{i}\delta\mathbf{l}_{i}\right). \end{array} \end{eqnarray} \tag{ B.11 }$$

Therefore, we arrive at the Hamiltonian

$$\begin{eqnarray} &&\fl H_{2} = -2\sum_{\mathbf{k}}\left\{\cos\phi\left[-\frac{3}{2}t_{2}-\frac{1}{2}\sum_{i}\delta t_{2i}\mp t_{2}\left(\mathbf{K}_{\pm}\cdot\sum_{i}\delta\mathbf{l}_{i}\right)\right](a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}+b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}) \right.\nonumber\\ &&+ \left.\sin\phi\left[\pm\frac{3\sqrt{3}}{2}t_{2}\pm\frac{\sqrt{3}}{2}\sum_{i}\delta t_{2i}- \frac{t_{2}}{2}\left(\mathbf{K}_{\pm}\cdot\sum_{i}\delta\mathbf{l}_{i}\right)\right](a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}-b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}) \right\}. \end{eqnarray} \tag{ B.12 }$$

We consider explicitly uniaxial strain as a study case, and the displacement of the l vectors is given by (see the strain tensor in [17])

$$\begin{eqnarray} \sum_{i}\delta\mathbf{l}_{i} &=& \left(\begin{array}{@{}c@{}} \sum_{i}\delta l_{ix} \\ \sum_{i}\delta l_{iy} \end{array}\right) = \sum_{i}\varepsilon\left(\begin{array}{@{}cc@{}} \varepsilon_{11} & \varepsilon_{12} \\ \varepsilon_{12} & \varepsilon_{22} \end{array} \right)\left(\begin{array}{@{}c@{}} l_{ix} \\ l_{iy} \end{array} \right) \nonumber\\ &=&\varepsilon\left(\begin{array}{@{}cc@{}} \varepsilon_{11}\sum_{i}l_{ix} & \varepsilon_{12}\sum_{i}l_{iy} \\ \varepsilon_{12}\sum_{i}l_{ix} & \varepsilon_{22}\sum_{i}l_{iy} \end{array}\right) =0. \end{eqnarray} \tag{ B.13 }$$

With this we infer the Hamiltonian including the next-nearest-neighbor hopping

$$\begin{equation} H_{2} = -2t_{2}\sum_{\mathbf{k}}\left[-\frac{1}{2}\bar{\eta}\cos\phi(a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}+b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}) \pm \frac{\sqrt{3}}{2}\bar{\eta}\sin\phi(a^{\dagger}_{\mathbf{k}}a_{\mathbf{k}}-b^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}) \right], \end{equation} \tag{ B.14 }$$

where $\bar {\eta }=3+\sum _{i}\frac {\delta t_{2i}}{t_{2}}$ is a dimensionless function varying with strain. The single-particle Hamiltonian reads

$$\begin{equation} h_{2}(\mathbf{k})=t_{2}\bar{\eta}\cos\phi\,\sigma_{0}-\sqrt{3}\tau_{z}t_{2}\bar{\eta}\sin\phi\,\sigma_{3}. \end{equation} \tag{ B.15 }$$

We use the relation derived in [17]

$$\begin{equation} t_{2}=V_{pp\pi}(l=a=\sqrt{3}a_{\mathrm{cc}})=t_{0}\exp\left[\frac{-3.37}{(l/a_{\mathrm{cc}}-1)}\right], \end{equation} \tag{ C.1 }$$

where t₀ is the nearest-neighbor hopping parameter. The deformed next-nearest-neighbor hopping is given by the variation of the length l

$$\begin{eqnarray} t'_{2} &=& t_{0}\exp\left[\frac{-3.37}{\left(\frac{l-\delta l}{a_{\mathrm{cc}}}-1\right)}\right] \nonumber\\ &\approx & t_{2}+\frac{3.37}{(\sqrt{3}-1)^{2}}t_{2}\frac{\delta l}{a_{\mathrm{cc}}}. \end{eqnarray} \tag{ C.2 }$$

Thus, we have

$$\begin{equation} \frac{\delta t_{2i}}{t_{2}}=\frac{3.37}{(\sqrt{3}-1)^{2}}\frac{|\delta\mathbf{l}_{i}|}{a_{\mathrm{cc}}}, \end{equation} \tag{ C.3 }$$

$$\begin{eqnarray} \delta\mathbf{l}_{i} &=& \varepsilon\left(\begin{array}{@{}cc@{}} \varepsilon_{11} & \varepsilon_{12} \\ \varepsilon_{12} & \varepsilon_{22} \end{array} \right)\left(\begin{array}{@{}c@{}} l_{ix} \\ l_{iy} \end{array}\right) \nonumber\\ &=&\varepsilon\left(\begin{array}{@{}cc@{}} \varepsilon_{11}l_{ix}, & \varepsilon_{12}l_{iy} \\ \varepsilon_{12}l_{ix}, & \varepsilon_{22}l_{iy} \end{array}\right), \end{eqnarray} \tag{ C.4 }$$

where δt_2i is the difference between the deformed t'₂ and t₂ associated with the vector l_i. So we have

$$\begin{eqnarray} &&\fl \sum_{i}\frac{\delta t_{2i}}{t_{2}} = \frac{3.37\varepsilon}{(\sqrt{3}-1)^{2}}\sum_{i}\left[(\varepsilon^{2}_{11} +\varepsilon^{2}_{12})\left(\frac{l_{ix}}{a_{\mathrm{cc}}}\right)^{2} +2\varepsilon_{12}(\varepsilon_{11}+\varepsilon_{22})\left(\frac{l_{ix}}{a_{\mathrm{cc}}}\right)\left(\frac{l_{iy}}{a_{\mathrm{cc}}}\right) \right.\nonumber\\ &&+ \left.\left(\varepsilon^{2}_{11}+\varepsilon^{2}_{12}\right)\left(\frac{l_{iy}}{a_{\mathrm{cc}}}\right)^{2}\right]^{\frac{1}{2}}. \end{eqnarray} \tag{ C.5 }$$

In the case of a uniaxial strain, the strain tensor is given by

$$\begin{eqnarray} &&\varepsilon_{11}=\cos^{2}\theta-\sigma\sin^{2}\theta, \nonumber\\ &&\varepsilon_{12}=(1+\sigma)\cos\theta\sin\theta, \\ &&\varepsilon_{22}=\sin^{2}\theta-\sigma\cos^{2}\theta,\nonumber \end{eqnarray} \tag{ C.6 }$$

where σ = 0.165 is the Poisson ratio [32] and θ is the angle between the tension direction and the zigzag direction.

The Hamiltonian of the nearest-neighbor hopping under strain can be derived as

$$\begin{equation} H_{1}=-\sum_{\mathbf{k}}\left[\underbrace{\sum_{j}t'_{0j}\mathrm{e}^{\mathrm{i}\boldsymbol{\kappa}\cdot\boldsymbol{\delta}'_{j}}}_{(\mathrm{III})}a^{\dagger}_{\mathbf{k}}b_{\mathbf{k}}+\mathrm{H.c.}\right], \end{equation} \tag{ D.1 }$$

where κ is the momentum vector with respect to the Γ point, $t'_{0j}=t_{0}+\triangle t_{0j}$ and $\triangle t_{0j}$ is the variation of the hopping integral induced by the strain. In deriving this equation, we have assumed that the BZ is maintained, but the structure in real space and the hopping integrals are deformed. In principle, BZ and the vectors in momentum space are also deformed. In some sense, the distortion of the vectors in momentum space can be compensated by the deformed vectors in real space. Thus, we assume here that the structure in the BZ is not changed for simplicity

$$\begin{eqnarray} &&\fl (\mathrm{III}) = \underbrace{t_{0}\sum_{j}\mathrm{e}^{\mathrm{i}(\mathbf{K}_{\pm}+\mathbf{k})\cdot\boldsymbol{\delta}_{j}}}_{(\mathrm{III})_{1}} +\underbrace{t_{0}\sum_{j}\mathrm{i}(\mathbf{K}_{\pm}+\mathbf{k})\cdot\triangle\boldsymbol{\delta}_{j} \mathrm{e}^{\mathrm{i}(\mathbf{K}_{\pm}+\triangle\mathbf{k})\cdot\boldsymbol{\delta}_{j}}}_{(\mathrm{III})_{2}} \nonumber\\ &&+ \underbrace{\sum_{j}\triangle t_{0j}\mathrm{e}^{\mathrm{i}(\mathbf{K}_{\pm}+\mathbf{k})\cdot\boldsymbol{\delta}_{j}}}_{(\mathrm{III})_{3}} +\underbrace{\sum_{j}\triangle t_{0j}\mathrm{i}(\mathbf{K}_{\pm}+\mathbf{k})\cdot\triangle\boldsymbol{\delta}_{j} \mathrm{e}^{\mathrm{i}(\mathbf{K}_{\pm}+\triangle\mathbf{k})\cdot\boldsymbol{\delta}_{j}}}_{(\mathrm{III})_{4}}, \end{eqnarray} \tag{ D.2 }$$

where k is the momentum vector with respect to K_±. (III)₁ can be derived as usual

$$\begin{equation} (\mathrm{III})_{1}=v_{F}\left(\mp k_{x}+\mathrm{i}k_{y}\right), \end{equation} \tag{ D.3 }$$

where $v_{F}=\frac {3}{2}a_{0}t_{0}$, a₀ is the distance of the carbon–carbon atoms, and t₀ is the nearest-neighbor hopping without strain. We have

$$\begin{eqnarray} &&\fl (\mathrm{III})_{2} \approx \underbrace{\mathrm{i}t_{0}\sum_{j}\left(\mathbf{K}_{\pm}\cdot\triangle\boldsymbol{\delta}_{j}\right) \mathrm{e}^{\mathrm{i}\mathbf{K}_{\pm}\cdot\boldsymbol{\delta}_{j}}}_{(\mathrm{III})_{2a}}+ \underbrace{\mathrm{i}t_{0}\sum_{j}\mathrm{i}\left(\mathbf{K}_{\pm}\cdot\triangle\boldsymbol{\delta}_{j}\right)(\mathbf{k}\cdot\boldsymbol{\delta}_{j}) \mathrm{e}^{\mathrm{i}\mathbf{K}_{\pm}\cdot\boldsymbol{\delta}_{j}}}_{(\mathrm{III})_{2b}}\nonumber\\ &&+ \underbrace{\mathrm{i}t_{0}\sum_{j}(\mathbf{k}\cdot\triangle\boldsymbol{\delta}_{j})\,\mathrm{e}^{\mathrm{i}\mathbf{K}_{\pm}\cdot\boldsymbol{\delta}_{j}} }_{(\mathrm{III})_{2c}}. \end{eqnarray} \tag{ D.4 }$$

After some algebra, we have

$$\begin{eqnarray} (\mathrm{III})_{2a} &=& \mathrm{i}v_{F}\varepsilon(\pi_{1}+\mathrm{i}\tau_{z}\pi_{2}), \nonumber\\ (\mathrm{III})_{2b} &=& v_{F}\varepsilon\frac{a_{0}}{2}\left[\pi_{2}k_{x}-\pi_{1}k_{y}+\mathrm{i}\tau_{z}(\pi_{1}k_{x}+\pi_{2}k_{y})\right], \\ (\mathrm{III})_{2c} &=& v_{F}\varepsilon\left[-\tau_{z}(\varepsilon_{11}k_{x}+\varepsilon_{12}k_{y})+\mathrm{i}(\varepsilon_{12}k_{x}+\varepsilon_{22}k_{y})\right],\nonumber \end{eqnarray} \tag{ D.5 }$$

where τ_z = 1 (−1) for K₊ (K₋), π₁ = K_±,xε₁₂ + K_±,yε₂₂ and π₂ = K_±,xε₁₁ + K_±,yε₁₂. A straight calculation leads to

$$\begin{eqnarray} && \hspace{-12pt} (\mathrm{III})_{3} = (w_{1}+w_{2})+\mathrm{i}(w_{3}+w_{4}), \nonumber\\ && w_{1} = \triangle t_{01}-\triangle t_{+}/2, \nonumber\\ && w_{2} = \tau_{z}\frac{a_{0}}{2}\left(-\frac{3}{2}\triangle t_{+}k_{x}+\frac{\sqrt{3}}{2}\triangle t_{-}k_{y}\right), \nonumber\\ && w_{3} = \tau_{z}\frac{\sqrt{3}}{2}\triangle t_{-}, \\ && w_{4} = \frac{a_{0}}{2}\left[-\frac{\sqrt{3}}{2}\triangle t_{-}k_{x}+\left(2\triangle t_{01}+\triangle t_{+}/2\right)k_{y}\right], \nonumber\\ && \hspace{-6pt} \triangle t_{\pm} = \triangle t_{02}\pm\triangle t_{03}.\nonumber \end{eqnarray} \tag{ D.6 }$$

We ignore (III)₄ since the lowest order is the second order in a small change and independent of the momentum. If we link the formalism derived in the tight-binding method to that derived from the theory of invariants, we ignore the terms that are independent of k_x(y), i.e. (III)_2a, and the term including the deformed momentum k_x(y), i.e. (III)_2c. Then we re-derived the d₁ and d₂ vectors in (1), and link the parameters in this equation by

$$\begin{eqnarray} &&\begin{array}{llll} \displaystyle \frac{\tilde{v}_{x}}{v_{F}} = -\frac{2\pi}{3\sqrt{3}}\varepsilon\varepsilon_{11}+\frac{\triangle t_{02}+\triangle t_{03}}{2t_{0}}, \\ \displaystyle \frac{\tilde{v}_{y}}{v_{F}} = \frac{2\pi}{3\sqrt{3}}\varepsilon\varepsilon_{11}-\frac{1}{\sqrt{3}}\frac{\triangle t_{02}-\triangle t_{03}}{2t_{0}}, \\ \displaystyle \mathcal{A}_{x} = \frac{1}{2}\left(\triangle t_{02}+\triangle t_{03}-2\triangle t_{01}\right), \\ \mathcal{A}_{y} = \frac{\sqrt{3}}{2}\left(\triangle t_{02}-\triangle t_{03}\right). \end{array} \end{eqnarray} \tag{ D.7 }$$

$\triangle t_{0i}$ can be calculated as done in appendix C, and we list the equation here

$$\begin{equation} \fl \frac{\triangle t_{0i}}{t_{0}}=\frac{3.37}{(\sqrt{3}-1)^{2}}\varepsilon\left[(\varepsilon_{11}^{2}+\varepsilon_{12}^{2})\left(\frac{\delta_{ix}}{a_{0}}\right)^{2} +2\varepsilon_{12}(\varepsilon_{11}+\varepsilon_{22})\frac{\delta_{ix}}{a_{0}}\frac{\delta_{iy}}{a_{0}} +(\varepsilon_{12}^{2}+\varepsilon_{22}^{2})\left(\frac{\delta_{iy}}{a_{0}}\right)^{2}\right], \end{equation} \tag{ D.8 }$$

where δ_ix(y) are the x (y) component of the vector δ_i.