Deformation induced pseudo-magnetic fields in complex carbon architectures

TL;DR

This paper demonstrates that complex carbon architectures like graphyne exhibit deformation physics similar to graphene, with emergent Dirac-Weyl equations and strain-induced fields, extending known phenomena to new materials.

Contribution

The authors develop two continuum theories for deformation in complex graphyne structures, showing their low-energy physics closely resembles that of graphene.

Findings

Deformation induces pseudo-magnetic fields in complex graphyne structures.

The low-energy physics of graphyne is described by an emergent Dirac-Weyl equation.

Key phenomena like pseudospin polarization and valley currents are present in these materials.

Abstract

We show that the physics of deformation in $\alpha$-, $\beta$-, and $6,6,12$-graphyne is, despite their significantly more complex lattice structures, remarkably close to that of graphene, with inhomogeneously strained graphyne described at low energies by an emergent Dirac-Weyl equation augmented by strain induced electric and pseudo-magnetic fields. To show this we develop two continuum theories of deformation in these materials: one that describes the low energy degrees of freedom of the conical intersection, and is spinor valued as in graphene, and one describing the full sub-lattice space. The spinor valued continuum theory agrees very well with the full continuum theory at low energies, showing that the remarkable physics of deformation in graphene generalizes to these more complex carbon architectures. In particular, we find that deformation induced pseudospin polarization and valley current loops, key phenomena in the deformation physics of graphene, both have their counterpart in these more complex carbon materials.

Figures

Equations

$$H_{TB}$$

$${\left[H(\mathbf{r},\mathbf{p})\right]}_{\alpha \beta }=\frac{1}{V_{UC}}\sum _j{M}_{j\alpha \beta }{\eta }_{\alpha \beta }(\mathbf{r},{\mathbf{K}}_j+\mathbf{p}/\hslash ),$$

$${\eta }_{\alpha \beta }(\mathbf{r},\mathbf{q})=\int d\mathbf{\delta }{e}^{-i\mathbf{q}.\mathbf{\delta }}{t}_{\alpha \beta }(\mathbf{r},\mathbf{\delta })$$

$$M_{j\alpha \beta }=e^{i{\mathbf{G}}_j.({\mathbf{\nu }}_{\alpha }-{\mathbf{\nu }}_{\beta })}$$

$$\mathbf{j}({\mathbf{R}}_i)=\frac{1}{V_{UC}}\sum _{\alpha }\left[n({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha })\mathbf{v}+\mathbf{v}n({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha })\right]$$

$$\mathbf{v}=\frac{1}{i\hslash }\left[\mathbf{r},H\right],$$

$$\mathbf{j}(\mathbf{r})=\Im \Psi (\mathbf{r}{)}^{†}\left[{\nabla }_{\mathbf{p}}H(\mathbf{r},\mathbf{p})\Psi (\mathbf{r})\right]$$

$$|\Psi ⟩=\sum _{{\mathbf{k}}_1\alpha }{c}_{{\mathbf{k}}_1\alpha }|{\Phi }_{{\mathbf{k}}_1\alpha }⟩$$

$$|{\Phi }_{{\mathbf{k}}_1\alpha }⟩=\frac{1}{\sqrt{N}}\sum _{{\mathbf{R}}_i}{e}^{i{\mathbf{k}}_1.({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha })}|{\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha }⟩$$

$$\displaystyle\frac{1}{2V_{UC}}\sum_{{\bf k}_{1},{\bf k}_{2},\alpha^{\prime},\beta}\Big{[}c_{{\bf k}_{1}\alpha^{\prime}}^{\ast}c_{{\bf k}_{2}\beta}\left\langle\Phi_{{\bf k}_{1}\alpha^{\prime}}\left|\hat{n}({\bf R}_{i}+\nu_{\alpha})\hat{{\bf v}}\right|\Phi_{{\bf k}_{2}\beta}\right\rangle$$

$$\displaystyle+\text{h.c.}\Big{]}$$

$$\frac{-i}{N\hslash }\sum _{{\mathbf{R}}_j}{e}^{-i{\mathbf{k}}_1.({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha })}{e}^{i{\mathbf{k}}_2.({\mathbf{R}}_j+{\mathbf{\nu }}_{\beta })}\left[{\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha }-{\mathbf{R}}_j-{\mathbf{\nu }}_{\beta }\right]$$

$$\times t_{\alpha \beta }({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha },{\mathbf{R}}_j+{\mathbf{\nu }}_{\beta })$$

$$\sum _{{\mathbf{R}}_j}f({\mathbf{R}}_j+{\mathbf{\nu }}_{\beta })=\frac{1}{V_{UC}}\sum _{{\mathbf{G}}_j}\hat{f}({\mathbf{G}}_j)e^{i{\mathbf{G}}_j.{\mathbf{\nu }}_{\beta }}$$

$$f(\mathbf{r})=e^{i{\mathbf{k}}_2.\mathbf{r}}\left({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha }-\mathbf{r}\right)t_{\alpha \beta }({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha },\mathbf{r})$$

$$\hat{f}(\mathbf{q})=\int d\mathbf{r}{e}^{-i(\mathbf{q}-{\mathbf{k}}_2).\mathbf{r}}\left({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha }-\mathbf{r}\right)t_{\alpha \beta }({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha },\mathbf{r})$$

$$\frac{e^{i({\mathbf{k}}_2-{\mathbf{k}}_1).({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha })}}{V\hslash }\sum _{{\mathbf{G}}_j}{e}^{i{\mathbf{G}}_j.({\mathbf{\nu }}_{\alpha }-{\mathbf{\nu }}_{\beta })}{\nabla }_{\mathbf{q}}{t}_{\alpha \beta }({\mathbf{R}}_i+{\mathbf{\nu }}_{\alpha },{\mathbf{k}}_2+{\mathbf{G}}_j)$$

$$\frac{1}{\sqrt{V}}{e}^{-i{\mathbf{p}}_1.\mathbf{r}}\left\{{\nabla }_{\mathbf{p}}{H}_{\alpha \beta }(\mathbf{r},\mathbf{p})\right\}\frac{1}{\sqrt{V}}{e}^{i{\mathbf{p}}_2.\mathbf{r}}.$$

$${\varphi }_{{\mathbf{p}}_1\alpha }(\mathbf{r})$$

$${\varphi }_{{\mathbf{p}}_2\beta }(\mathbf{r})$$

$$⟨{\Phi }_{{\mathbf{k}}_I\alpha }\left|\hat{n}(\mathbf{r})\hat{\mathbf{v}}\right|{\Phi }_{{\mathbf{k}}_J\beta }⟩={\varphi }_{{\mathbf{p}}_1\alpha }^{†}(\mathbf{r})\left[{\nabla }_{\mathbf{p}}H(\mathbf{r},\mathbf{p}){\varphi }_{{\mathbf{p}}_2\beta }(\mathbf{r})\right]$$

$$\Psi (\mathbf{r})=\sum _{{\mathbf{p}}_1}{c}_{{\mathbf{p}}_1\alpha }{\varphi }_{{\mathbf{p}}_1\alpha }(\mathbf{r})$$

$$\Psi (\mathbf{r}{)}^{†}\mathbf{j}\Psi (\mathbf{r})=\frac{1}{2}\left\{\Psi (\mathbf{r}{)}^{†}\left[{\nabla }_{\mathbf{p}}H\Psi (\mathbf{r})\right]+{\left[{\nabla }_{\mathbf{p}}H\Psi (\mathbf{r})\right]}^{†}\Psi (\mathbf{r})\right\}$$

$$\mathbf{j}(\mathbf{r})=\frac{1}{2m}\left[\Psi (\mathbf{r}{)}^{\ast }\mathbf{p}\Psi (\mathbf{r})-\Psi (\mathbf{r})\mathbf{p}\Psi (\mathbf{r}{)}^{\ast }\right]$$

$$\mathbf{j}(\mathbf{r})=\Psi (\mathbf{r}{)}^{†}{v}_F\mathbf{\sigma }\Psi (\mathbf{r})$$

$$H$$

$$j_j=V_F^{ij}(\mathbf{r}){\sigma }_i+i{W}_F^{ij}(\mathbf{r}){\sigma }_i+M^{ijk}(\mathbf{r}){\sigma }_i{p}_j.$$

$$\nabla .\mathbf{j}(\mathbf{r})$$

$$p_j{M}^{ijk}{\sigma }_i{p}_k-2i{W}^{ij}{\sigma }_i{p}_j$$

$$p_j{V}^{ij}{\sigma }_i-2i{\Gamma }_i{\sigma }_i$$

$$\bm{\nabla}.{\bf j}({\bf r})=\frac{1}{i\hbar}\Im\Bigl{[}-2\Psi^{\dagger}H\Psi\Bigr{]}$$

$${\partial }_{t}n(\mathbf{r})=\frac{1}{i\hslash }\left[{\Psi }^{†}\left(H\Psi \right)-{\left(H\Psi \right)}^{†}\Psi \right]$$

Full text

I Introduction

One of the most remarkable features of graphene is the robustness of its ultra-relativistic low energy physics. The Dirac-Weyl Hamiltonian that describes the quasi-particles of pristine graphene remains a valid description even under substantial deformation, simply being augmented by effective pseudo-magnetic and electric fields that encode the deformation in a low energy description. These fields, except for the requirement that the pseudo-magnetic field change sign at conjugate valleys, behave exactly as physical electric and magnetic fieldsgui10 ; Electronic , resulting in a rich phenomenology of deformation induced physics in single layer graphene. For realistic strain the induced magnetic field can reach hundreds of Tesla, a remarkable effect observable as deformation induced Landau levels in grapheneScience-NB . This deep connection between structural deformation and an induced electromagnetic field promises a control over electronic properties unrivaled in any three dimensional material, and generates novel physical effects such as deformation induced valley filtersValley-filter ; Valley-filter2 ; Valley-filter6 ; Valley-filter7 ; Valley-filter8 , and psuedospin polarizationpseudo1 ; pseudo2 .

Following the experimental realization of graphene low energy conical intersections have been predicted for several all-carbon materials, each with a substantially more complex lattice structure than that of graphene6 ; 5 ; 7 ; 8 ; 13 ; s1 ; s2 ; s5 ; s7 ; s8 ; s9 ; s11 ; s12 ; s13 ; s14 ; 16 ; 17 ; 23 ; 24 ; 25 ; 26 ; 27 ; 42 ; 46 ; 53 ; 54 ; 61 ; 43 ; r1 ; r2 ; 1 . For example, $6,6,12$-graphyne and $\beta$-graphyne both possess 18 carbon atoms in their unit cell, as opposed to the 2 atom unit cell honeycomb lattice of graphene. This entails a much more difficult chemistry of their fabrication5 ; 9 ; 44 , but also a fundamental difference in the low energy physics. In graphene sub-lattice space is isomorphic to $SU(2)$ pseudo-spin space, and it is this that underpins the connection between deformation and effective electric and magnetic fields. However, this simple relation between the pseudospin and sub-lattice degrees of freedom is lost in these more complex materials. While the structural physics of non-uniform deformation in the graphynes has been quite intensively investigateds8 ; s9 ; s11 ; s13 , the corresponding attention has not been devoted to the electronic theory of general non-uniform deformations, with most electronic investigations focusing either on uniform uniaxial and biaxial strainss1 ; s2 ; s5 ; s12 ; s14 ; 49 ; 50 , “rotating” strains7 , or phononic excitationss4 . A natural question is therefore: how much of the rich electronic physics of non-uniform deformation in graphene finds a counterpart in these more complex carbon architectures?

The purpose of the present paper is to answer this question. To that end we generalize the continuum theory of deformation in graphene to materials with arbitrary numbers of atoms in the unit cell. As a minimal description of the electronic structure entails one $\pi$-orbital per basis atom, this theory now necessarily takes two forms: one involving all sub-lattice degrees of freedom, and a down-folded theory describing only the spinor degree of freedom of the low energy Dirac cone. For all three graphynes we find that this latter description is, at low energies, in very close agreement with the full continuum theory. Thus the intimate connection between structure and effective electromagnetic field is preserved in these more complex architectures, and the rich physics of deformation in graphene generalizes to the semi-metallic graphynes.

II Continuum theory of complex carbon materials

While the continuum theory of deformation in graphene is very well developed, the same attention has not been paid to the graphynes. The complex lattice structures of these materials (see Fig. 1) render either cumbersome (the expansion of a tight-binding Hamiltonian) or inapplicable (transport of the Dirac-Weyl equation to non-Minkowski metric) the methods used in deriving the continuum theory of deformation for graphene. To avoid this we employ a methodologykiss15 ; shall2016 ; shall17 ; Gupta2018 ; Rost2019 based on an exact map of the tight-binding Hamiltonian to a general continuum operator, recently used to treat stacking deformations in bilayer structurekiss15 ; shall17 , and acoustic and optical deformations in single layer grapheneGupta2018 . We first briefly review this theory, before extending it to include a treatment of currents, and subsequently describing both the complete and down-folded versions required to investigate deformation in graphyne.

II.1 An exact map of the tight-binding Hamiltonian

In this approach the “input” consists simply of a hopping function envelope $t_{\alpha\beta}({\bf r},\bm{\delta})$ describing the tight-binding hopping matrix amplitude from a position ${\bf r}$ on sub-lattice $\alpha$ to ${\bf r}+\bm{\delta}$ on sub-lattice $\beta$ in the material, with the change in electron hopping due to deformation encoded in the ${\bf r}$-dependence of this function. The tight-binding Hamiltonian is therefore

[TABLE]

where ${\bf R}_{i}$ and $\bm{\nu}_{\alpha}$ are the lattice and basis vectors of the underlying high symmetry lattice.

The corresponding continuum Hamiltonian $H({\bf r},{\bf p})$ that is exactly equivalent to this tight-binding Hamiltonian is given byRost2019

[TABLE]

where $V_{UC}$ is the unit cell volume, the sum is over all ${\bf G}_{j}$ vectors of the reciprocal lattice, ${\bf K}_{j}={\bf K}_{0}+{\bf G}_{j}$ with ${\bf K}_{0}$ a reference momentum in the Brillouin zone (for the graphynes we study this will be the momentum at which the conical intersection occurs),

[TABLE]

the so-called mixed space hopping function, and the $M_{j}$ matrices are given by

[TABLE]

As shown in Ref. Gupta2018, expansion of the ${\bf r}$-dependence of $\eta_{\alpha\beta}({\bf r},{\bf q})$ for slowly varying deformations and expansion of the ${\bf p}$-dependence for momenta near the conical intersection leads to a systematic treatment of deformation within the continuum picture.

II.2 An exact map of the current operator

In this work a focus will be on the valley current densities induced by deformation, and thus we require a current operator corresponding to the Hamiltonian Eq. 2. The current operator in tight-binding theory, the so-called bond current operator, is given bytod02 ; boy10

[TABLE]

where $\hat{n}({\bf R}_{i}+\bm{\nu}_{\alpha})=\left|{\bf R}_{i}+\bm{\nu}_{\alpha}\right\rangle\left\langle{\bf R}_{i}+\bm{\nu}_{\alpha}\right|$ is the density operator with the velocity operator

[TABLE]

and where we have included normalization by the unit cell volume $1/V_{UC}$.

One would expect that the method used to derive Eq. (2) from Eq. (1) should, if applied to Eq. (5), yield a current density given by

[TABLE]

with $H({\bf r},{\bf p})$ given by Eq. (2), i.e. in the continuum limit the relation between Hamiltonian and current operator should follow from Hamilton’s equations. However, given the very different forms of the tight-binding Hamiltonian Eq. (1) and the bond-current operator Eq. (5) in tight-binding theory it is not obvious that this is the case.

Indeed, the precise link between the tight-binding and continuum limits of the current operator has been the subject of recent discussion in the context of recovering the Schödinger current operator from the bond current formulaeboy10 for a quadratic band model. As noted in Ref. boy10, the fundamental difficultly involves taking the continuum limit of a discrete grid and the associated ambiguity in the definition of differential operators. Employing the same methodology involved in the derivation of Eq. 2, however, we find a general form for the continuum current operator precisely equivalent to the general tight-binding bond-current formula while avoiding all use of grid limits. This turns out to be just the intuitive result given by Eq. (7) from which we can then recover, as specific cases, both the Schrödinger and Dirac-Weyl current operators.

To derive this result we first consider the matrix element of Eq. 5 with a general state of the system

[TABLE]

where $\left|\Phi_{{\bf k}_{1}\alpha}\right\rangle$ denotes a Bloch function of sublattice $\alpha$ and crystal momentum ${\bf k}_{1}$:

[TABLE]

One finds for the matrix element $\left\langle\Psi\left|{\bf j}({\bf R}_{i})\right|\Psi\right\rangle$ the result

[TABLE]

We now work out the matrix element $\left\langle\Phi_{{\bf k}_{1}\alpha^{\prime}}\left|\hat{n}({\bf R}_{i}+\nu_{\alpha})\hat{{\bf v}}\right|\Phi_{{\bf k}_{2}\beta}\right\rangle$ in detail. Insertion of the Bloch functions, Eq. 9, yields

[TABLE]

where $t_{\alpha\beta}({\bf R}_{i}+\bm{\nu}_{\alpha},{\bf R}_{j}+\bm{\nu}_{\beta})=\left\langle{\bf R}_{i}+\bm{\nu}_{\alpha}\left|H\right|{\bf R}_{j}+\bm{\nu}_{\beta}\right\rangle$ is the usual tight-binding hopping matrix element. To derive a continuum limit we now employ the Poisson sum formula in the form

[TABLE]

for which the appropriate function $f({\bf r})$ is

[TABLE]

with the Fourier transform

[TABLE]

Evaluation of the integral and employing the Poisson sum in Eq. (11) then yields

[TABLE]

which, by defining a reference momenta ${\bf K}_{0}$ through ${\bf k}_{2}+{\bf G}_{j}={\bf G}_{j}+{\bf K}_{0}+{\bf p}_{2}$ and promotion of ${\bf p}_{2}$ to an operator, can straightforwardly be recast into a form involving the continuum Hamiltonian $H({\bf r},{\bf p})$ Eq. (2):

[TABLE]

To make the connection with a continuum Hamiltonian we introduce the vector plane waves

[TABLE]

where $\left|\alpha\right\rangle$ and $\left|\beta\right\rangle$ represent unit vectors in a space of dimension equal to the number of atomic degrees of freedom (for graphene these would just be pseudospin up and pseudospin down states). Employing these functions we then arrive at the desired operator equivalence:

[TABLE]

Insertion of this result back into Eq. 10 and the obvious definition for the system wavefunction in the continuum representation as

[TABLE]

then leads to our final result

[TABLE]

which is evidently the sought for intuitive form given by Eq. (7).

This expression trivially reproduces both the well known current operators for the Dirac-Weyl and Schrödinger equation. For a Schrödinger form $H=\frac{1}{2m}p^{2}$ we find

[TABLE]

whereas for the Dirac-Weyl Hamiltonian $H=v_{F}\bm{\sigma}.{\bf p}$ we have

[TABLE]

There is, however, an important caveat to Eq. (21). While the bond current operator always satisfies a discrete form of the continuity equation appropriate for the tight-binding Hamiltoniantod02 the continuum version of the bond current, Eq. (21), is not guaranteed to satisfy the continuity equation. This follows as it is simply the expectation value of the velocity operator. Curiously, it turns out that Eq. (21) indeed violates the continuity equation, but only for Hamiltonians containing a higher than second power in momentumCE , a common occurrence in the effective Hamiltonians of condensed matter. Other less intuitive definitions have been providedCorrect-current-operator1 ; Correct-current-operator2 ; definition and these can encode non-classical current contributions. However the breakdown of a classical relationship between velocity and current should probably be viewed as a failure of effective Hamiltonian theory.

II.3 The continuity equation for deformation in Dirac-Weyl materials

As our focus here will be on deformation in materials with low energy conical intersections, and the valley currents that deformation induces, it is useful to demonstrate that the current operator defined in the previous section indeed satisfies the continuity equation for such systems. To this end we consider the most general form of the Dirac-Weyl equation augmented by deformation induced fields which, up to second order in momentum, is given by

[TABLE]

In this expression there are two effective gauge field terms: ${\bf A}$, a real valued gauge2002-paper ; ref13 ; PhysRevLett.108.227205 ; non-uniform-strain ; AMORIM20161 ; Peeters-revisited that transforms in a complex way under spatial rotations, and $\bm{\Gamma}$, an imaginary gaugePhysRevLett.108.227205 that transforms as a field term in the Dirac-Weyl operator. $v_{F}^{ij}$ and $w_{F}^{ij}$ are the real PhysRevLett.108.227205 ; non-uniform-strain ; AMORIM20161 ; Peeters-revisited and imaginary velocity tensors, in which the $i$ index runs over the three Pauli matrices $\sigma_{0}$, $\sigma_{1}$, and $\sigma_{2}$ and the $j$ index over the two degrees of freedom of space. Finally there is also a trigonal warping term $M^{ijk}$ in which again $i$ runs over the Pauli matrices and the indices $j,k$ over the two degrees of spatial freedom. Higher orders of momenta can be included in Eq. (24) but, if Hermiticity is to be preserved, only by restricting the spatial fields to be slowly varyingGupta2018 .

It is not immediately evident that Eq. (24) satisfies the continuity equation $\bm{\nabla}.{\bf j}({\bf r})+\partial_{t}n({\bf r})=0$ (as both the Hamiltonians $v_{F}\bm{\sigma}.{\bf p}$ and $v_{F}\bm{\sigma}.{\bf p}+\sigma_{i}A_{i}$ obviously do) due both to the presence of the coordinate dependent velocity and mass tensors, as well as the fact that both gauge and velocity terms have both real and imaginary parts. However, as we now show using conditions that guarantee hermiticity of the deformation Hamiltonian, the continuity equation is indeed satisfied.

The $j$’th component of the corresponding current operator is given by

[TABLE]

with the divergence of the current density then given by

[TABLE]

Using the obvious relations $\Im[(p_{j}\Psi)^{\dagger}M^{ijk}\sigma_{i}p_{j}\Psi]=0$ and $(\psi^{\dagger}V^{ij}\sigma_{i}p_{j}\Psi)^{\dagger}=(V^{ij}\sigma_{i}p_{j}\Psi)^{\dagger}\Psi$ in conjunction with the hermiticity conditions obeyed by the HamiltonianPhysRevLett.108.227205 ; Gupta2018

[TABLE]

we have

[TABLE]

and using

[TABLE]

we then find

[TABLE]

and so the the continuity equation is satisfied.

III Continuum theory of $\alpha$-,$\beta$-, and $6,6,12$-graphyne

III.1 Tight-binding parameterization

A continuum theory of the pristine lattices, i.e. without deformation, is easily obtained from Eq. (2). In this case the mixed space hopping function, Eq. (3), loses its ${\bf r}$-dependence becoming simply the Fourier transform of the electron hopping function defined between each sub-lattice: $\hat{t}_{\alpha\beta}({\bf q})=\int d\bm{\delta}\,e^{-i{\bf q}.\bm{\delta}}t_{\alpha\beta}(\bm{\delta})$. We therefore require a functional form of $t_{\alpha\beta}(\bm{\delta})$ describing electron hopping for the three graphynes we consider: $\alpha$-, $\beta$-, and $6,6,12$-graphyne. To obtain such a form we take the tight-binding parameters of Ref. liu12, , which are defined only for the nearest neighbour hopping, and fit to a Gaussian functional form

[TABLE]

such that the nearest neighbour hopping is reproduced, with negligible next nearest (and further) hopping.

This Gaussian form is useful as (i) it renders the Fourier transform straightforward and (ii) deformation, which modifies the hopping matrix element through changes of the hopping vector, is easily included. The tight-binding band structures using this parameterization are shown as the full lines in Fig. 2. The expected low energy conical intersections are found: (i) at the K-point for $\alpha$-graphyne, (ii) on the line connecting the $\Gamma$ and M-points for $\beta$-graphyne, and (iii) two low energy cones one at the X-point (we denote this cone I) and a second on the line connecting $\Gamma$- and Y-points (denoted cone II) for $6,6,12$-graphyne. In fact, cone I is shifted somewhat further off the X-point as compared to ab-initio calculations (a line fraction of 0.2 as opposed to 0.08), and cone II somewhat closer to the $\Gamma$ point (a line fraction of 0.5 as opposed to 0.62 in ab-initio). Furthermore, with only nearest neighbour hopping the cone tilting found in ab-initio for $\beta$-graphyne and cone II in $6,6,12$-graphyne is not reproduced, lowering the value of $B$ somewhat in Eq. (32) we find restores the cone tilting. We have checked that sensible variation of the tight-binding parameters does not significantly change the results we present in subsequent sections for deformation in these materials.

III.2 Continuum theory for pristine lattices

To extract a tractable continuum description from Eq. (2) requires a Taylor expansion in momentum about the Dirac point. Expanding the hopping function in Eq. (2) to first order

[TABLE]

generates an expression of the form

[TABLE]

where the Hamiltonian at the Dirac momenta is

[TABLE]

and the matrices $H_{i}^{(1)}$ given by

[TABLE]

where $K_{ji}$ is the $i$’th component of the vector ${\bf K}_{j}$. The matrices Eqs. (35)-(36) are labelled by sub-lattice indices, i.e. represent $8\times 8$ matrices for $\alpha$-graphyne and $18\times 18$ matrices for $\beta$-graphyne and $6,6,12$-graphyne. These therefore describe not only the low energy conical intersection (described generically by a spinor degree of freedom) but also all other band manifolds. In Fig. 2 is shown the band structure obtained from Eq. (34), with additionally second order in momentum terms included, revealing an excellent agreement with the underlying tight-binding method.

To obtain a description of the low energy conical intersection we must down-fold the full continuum theory. To that end we diagonalize the Hamiltonian at the Dirac point $H^{(0)}$, Eq. (35), and apply the resulting unitary transformation $U$ to the full Hamiltonian $H^{full}_{0}({\bf p})$. This yields

[TABLE]

where $\bm{\epsilon}$ is a diagonal matrix whose entries are the eigenvalues at the Dirac point momenta, with the matrices $H^{(1)}$ now encoding hybridization of these bands at a finite momentum away from the Dirac point. This now allows us to identify the sub-space corresponding to degenerate eigenvalues at the Dirac point, and by retaining only this sub-space one arrives at a spinor valued low energy Hamiltonian. This procedure yields a Dirac-Weyl form, but with an SU(2) rotation. Undoing this with a further unitary transform $U_{s}$ we then find a generic final form for all three materials16 ; 46

[TABLE]

The band structure calculated using Eq. (38) is shown as the dot-dashed lines in Fig. 2, showing good agreement with the full tight-binding calculation for the low energy conical intersection. Further improvement requires additional bands to be included in the continuum theory, not higher orders of momenta.

While Eq. (38) provides a good description of the low energy manifold it does not directly provide the physical wavefunction; for this a back transformation is required. Given a spinor eigenvector $\phi$ of Eq. (38) we must firstly transform back to the global SU(2) frame

[TABLE]

which then provides the coefficients for constructing the physical wavefunction from the Dirac point wave-functions:

[TABLE]

where $\Psi_{i}$ is the $i$’th component of the physical wavefunction (i.e., in sub-lattice space) while the sum $j$ is over the band indices of the low energy manifold (recall that $U$ is the unitary transform that diagonalizes $H_{0}$ the Hamiltonian at the Dirac point momenta). In graphene such a back transformation is, of course, unnecessary; a low energy expansion directly yields a Dirac-Weyl equation. For any more complex material, however, the coefficients of the Dirac-Weyl spinor wavefunction merely parameterize the low energy conical intersection in terms of the Dirac point wavefunctions.

III.3 Continuum theory of deformation

Inclusion of deformation into the low energy scheme requires describing how electron hopping changes throughout the material, i.e. the full function $t_{\alpha\beta}({\bf r},\bm{\delta})$. This can be obtained simply from the geometric information of how the deformation changes the hopping $\bm{\delta}$ at point ${\bf r}$ in the material, via substitution of $\bm{\delta}$ in the hopping function of the high symmetry material $t_{\alpha\beta}(\bm{\delta})$ by $t_{\alpha\beta}(\bm{\delta}({\bf r}))$. The change in the square of the hopping vector due to the applied deformation field ${\bf u}({\bf r})$ is just

[TABLE]

which, upon substitution into a hopping function $t_{\alpha\beta}(\bm{\delta}^{2})$ yields, via a Taylor expansion for slowly varying fields, the mixed space hopping function

[TABLE]

In this expression $t_{\alpha\beta}^{(n)}({\bf q}^{2})$ are the the Fourier transforms of derivatives of the high symmetry hopping function:

[TABLE]

Insertion of this result directly into Eq. (2) leads to a Hamiltonian of the form

[TABLE]

with the matrices $H_{nm}$ given by

[TABLE]

where again $K_{jn}$ is the $n$’th component of the vector ${\bf K}_{j}={\bf K}_{0}+{\bf G}_{j}$. This along with Eq. 34 provides a continuum description in which all bands are included, i.e., in terms of all sub-lattice degrees of freedom.

To identify the low energy sector of the Hamiltonian with deformation, we employ the same transformations $U$ and $U_{s}$ that yielded the low energy sector for the pristine lattice. This, for all three materials, then results in the following generic form for the low energy deformation Hamiltonian

[TABLE]

where we have suppressed the ${\bf r}$-dependence is the deformation tensor and fields, and where the $f_{ij}$ and $g_{ij}$ are constants. This is very close to the form of deformation in graphene, which can obtained by substituting $f_{xx}=f_{yy}=\alpha$ and $g_{xx}=-g_{yy}=-g_{xy}=\beta$ into this expression. This latter condition,

[TABLE]

will hold for all materials in which the Dirac cone is protected by symmetry, as this ensures that biaxial strain cannot displace the Dirac cone of the high symmetry point. This is true for $\alpha$-graphyne, and we find Eq. 47 to be exactly satisfied by our deformation expansion, but not for $\beta$- or $6,6,12$-graphyne. For these materials biaxial deformation shifts the Dirac cone in momentum space; for $\beta$-graphyne towards the X point while for $6,6,12$-graphyne cone I moves towards the $X$ point and cone II towards the $Y$ point50 .

IV Electronic structure of non-uniform deformation in graphyne

Having established the basic theory we now examine the electronic consequences of non-uniform deformation. After a brief description of the numerical methodology, we first examine the robustness of the low energy conical manifolds under deformation, before considering the changes in electron densities and current densities induced by deformation.

IV.1 Numerical details

Solution of either the 2-band or full band Hamiltonian with deformation, Eqs. 44 and 46, is performed using a basis of eigenstates from the pristine system. For a momentum ${\bf k}_{0}$ the basis set is formed from all eigenstates of the pristine system with momenta $|{\bf k}-{\bf k}_{0}|<k_{cut}$ and energy $|\epsilon^{(0)}_{{\bf k}i}|<e_{cut}$. The advantage of this basis lies in the efficiency of convergence: we find that typically to converge the electronic structure in an energy window $E$ requires $e_{cut}\approx 1.5E$. The basis size is determined by $e_{cut}$ with $k_{cut}$ chosen so as the restrict the calculation to a single valley. We find that $e_{cut}$ chosen so that the basis size is of between 1000 and 2000 states is usually sufficient for convergence. For further numerical details we refer the reader to Ref. Gupta2018, .

IV.2 Spectral weight changes due to deformation

How are the band manifolds of the pristine lattices modified by deformation? As large scale periodic deformation results in a Brillouin zone reduced by a factor of $10^{2}$-$10^{3}$ as compared to the pristine lattice, the multiply back-folded and hybridized band structure becomes very hard to interpret. A much more useful quantity is what could be called a “poor man’s spectral function”:

[TABLE]

where

[TABLE]

In this expression $\left|\Psi_{{\bf k}j}\right\rangle$ is the $j$’th eigenstate at crystal momentum ${\bf k}$ of the system with deformation and $\left|\phi_{{\bf k}i}\right\rangle$ an eigenstate, also at ${\bf k}$, for the system without deformation. In the absence of deformation $\rho_{{\bf k}j}=1$ and Eq. 48 is simply the band structure of the pristine material. However in the presence of deformation $\left|\Psi_{{\bf k}j}\right\rangle$ will involve the coupling together of many eigenstates of the pristine system and $\rho_{{\bf k}j}<1$. Thus Eq. 48 represents how the bands of the high symmetry system are broadened through scattering induced by the deformation. To explore this we consider the deformation field shown in Fig. 3, similar to those typically employed in the discussion of non-uniform deformation in graphene. In panel (a) of this figure is shown the deformation field, while in panel (b) is shown the magnitude of the deformation tensor. The maximum value of the strain tensor is $\sim 7$%; our ab-initio calculations indicate that strains of $<7$% are within the elastic regime of these materials. As can be seen in Fig. 4 for each of the three graphynes we consider the “spectral function” follows closely the band structures of the pristine systems (see Fig. 2) but with the expected deformation induced broadening. The non-uniform “speckled” nature of the spectral intensity along the band lines can be understood as arising from the complex multiple intersections and subsequent hybridization that occurs when the bands are folded back to the Brillouin zone of the deformed system; in the extended zone scheme this will be manifest as a non-uniform spectral weight along the band lines.

IV.3 Charge inhomogenity and current flow

Having established the robustness of the low energy manifold to deformation, we now consider a description of deformation within a continuum theory of the low energy conical intersection. In this case the physics is encoded in the effective electric and pseudo-magnetic fields that augment the Dirac-Weyl equation, and these fields are shown in Fig. 5 for $\alpha$-, $\beta$- and $6,6,12$-graphyne, each with the same circularly symmetric deformation field shown in Fig. 3. Strikingly, for all three materials the form of the effective pseudo-magnetic and scalar fields is very similar. This is remarkable when one considers the very different lattice structures of these three systems, with $\alpha$- and $\beta$-graphyne possessing hexagonal lattices and 6,6,12-graphyne a rectangular lattice. While the pseudo-magnetic fields are comparable in magnitude, the scalar field is almost an order of magnitude greater for cone II of 6,6,12-graphyne than for $\alpha$-graphyne suggesting that the interplay of gauge and scalar fields, known to be significant for describing nanobubbles in graphene, would be especially important in this materialScience-NB ; PhysRevB.96.241405 .

Is the 2-band Dirac-Weyl theory of deformation valid in these more complex carbon architectures? To probe this question we now examine the deformation induced changes in electron density and current density calculated using the 2-band and full-band continuum theories. If the Dirac-Weyl continuum theory provides a valid description, then results from these two distinct continuum theories should be in close agreement. In rows 1 and 2 of Fig. 6 are shown electron density using the 2-band and full-band continuum theory respectively, with in rows 3 and 4 displayed the current density, again calculated from the 2- and full-band continuum theory respectively. For both density and current density it can be seen that the results of the two continuum theories are in very good agreement: the Dirac-Weyl description of deformation thus remains valid for these much more complex lattice structures.

We now consider the structure of the deformation induced changes to electron density and current density. For $\alpha$-graphyne these are very similar to those reported in previous studies of graphene for circularly symmetric deformation fields, with “charge flowers”Gaussian of $C_{6}$ symmetry and an associated $C_{3}$ symmetry current density patterneddy_currents ; currents2 ; currents3 , see panels 1a-4a of Fig. 6. For $\beta$- and $6,6,12$-graphyne, however, the density exhibits a much lower symmetry, especially striking for cone I of 6,6,12-graphyne. This reflects the $C_{2}$ lattice symmetry for 6,6,12-graphyne. The corresponding current densities are, however, much closer in form to those found in $\alpha$-graphyne and graphene.

As the effective fields induced by deformation showed no great distinction of form the differences in electron density indicates the importance of the transformation from the pseudo-spinor of the Dirac-Weyl equation back to the physical wavefunction discussed in Sec. II. However, the deformation induced current densities do correspond to the pseudo-magnetic fields shown in Fig. 5, with regions of strong current flow occurring at the nodal lines of the pseudo-magnetic field. The current density induced by deformation in these graphynes can therefore, just as in graphene, be understood as due to snake states2000 ; currents2 ; currents3 at regions where the pseudo-magnetic field changes sign.

Finally we examine the question of sub-lattice polarization. In graphene the deformation induced charge density is strongly polarized on sub-lattice A or B, which can be viewed in the Dirac-Weyl picture as a local pseudospin polarization. The more complex lattice structures for $\alpha$-, $\beta$-, and $6,6,12$-graphyne, which have 8, 18, and 18 atoms in the unit cell respectively, suggests that this physics will not be transferable to these more complex carbon architectures. Remarkably, as can be seen in Fig. 7, this is not the case. For each of these materials the atom projected density takes on one of only two forms, exhibited in the first two columns, with the correspondence between projection type and atom position in the unit cell shown in column three. However, while each of the 8 atoms of $\alpha$-graphyne have a density given exactly by one of these two projection types, for $\beta$-graphyne and 6,6,12-graphyne there are slight deviations amongst the 9 atom projections of each type.

IV.4 The role of optical deformation

Thus far we have not considered the role of local relaxations which will undoubtedly be induced by application of a deformation field. This is known to be more significant in these materials than in graphene (where it can also qualitatively change the physicsmidgap ). In the case of graphene the two atom unit cell leads to a natural effective Hamiltonian theory in terms of acoustic and optical deformation fields, as recently discussed by Gupta et al.Gupta2018 . However, the more complex unit cells of the graphynes imply many more optical modes. To simplify this situation, and by analogy with graphene, we can define optical and acoustic modes in terms of the two groups of atoms on which charge is localized due to pseudospin polarization. In this way we can define the average displacement of each group of atoms off their ideal positions under strain, and so define single optical and acoustic modes.

To investigate this we have performed ab-initio calculations for $6,6,12$-graphyne using the VASP software suite, in which we allow the 18 atoms of the unit cell to relax under an applied biaxial strain. In Fig. 8 is shown the magnitude of the resulting optical deformation, due to local atomic relaxation, given as a percentage of the nearest neighbour separation for a range of biaxial strains. As can be seen, for biaxial strains of up to 6%, the optical deformation is of the order of 1%. Beyond 8% biaxial strain the lattice substantially reconstructs. Thus optical deformation will likely play some role in the physics of these materials, as they do with graphene, and may in principle be included in an effective Hamiltonian description following the scheme outlined in Ref. Gupta2018 .

V Discussion

The principal question we have addressed is whether the intuitive Dirac-Weyl description of deformation in graphene generalizes to the more complex carbon architectures of the graphynes. To answer this question we have developed two distinct continuum theories: a Dirac-Weyl type theory, formally identical to that of graphene, and an continuum approach describing the full band structure. For $\alpha$-, $\beta$-, and $6,6,12$-graphyne these lead to very similar results for the deformation induced changes to the density and current density close to the Dirac point, showing that the Dirac-Weyl description remains valid for these materials. Deformation in the graphynes thus retains the remarkable connection between structural change and pseudo-magnetic and scalar fields found in graphene, and the rich physics of that material – valley filtersValley-filter ; Valley-filter2 ; Valley-filter6 ; Valley-filter7 ; Valley-filter8 , pseudospin polarizationpseudo1 ; pseudo2 , and deformation induced Landau laddersScience-NB – can be expected to be found in the graphynes, if they can be synthesized.

Acknowledgement

This work was carried out in the framework of SFB 953 of the Deutsche Forschungsgemeinschaft (DFG).