In-plane magnetoelectric response in bilayer graphene

TL;DR

This paper investigates the in-plane magnetoelectric response of bilayer graphene, linking it to valleytronics, and derives a comprehensive low-energy Hamiltonian including key magnetoelectric terms.

Contribution

It provides a detailed derivation of the low-energy Hamiltonian for bilayer graphene with in-plane fields, revealing an axion-like term and anisotropic corrections related to interlayer hopping.

Findings

Confirmed the existence of an axion-type pseudoscalar magnetoelectric term.

Identified anisotropic corrections related to skew interlayer hopping parameter γ4.

Suggested optical conductivity features to detect magnetoelectric effects.

Abstract

A graphene bilayer shows an unusual magnetoelectric response whose magnitude is controlled by the valley-isospin density, making it possible to link magnetoelectric behavior to valleytronics. Complementary to previous studies, we consider the effect of static homogeneous electric and magnetic fields that are oriented parallel to the bilayer's plane. Starting from a tight-binding description and using quasi-degenerate perturbation theory, the low-energy Hamiltonian is derived including all relevant magnetoelectric terms whose prefactors are expressed in terms of tight-binding parameters. We confirm the existence of an expected axion-type pseudoscalar term, which turns out to have the same sign and about twice the magnitude of the previously obtained out-of-plane counterpart. Additionally, small anisotropic corrections to the magnetoelectric tensor are found that are fundamentally related to the skew interlayer hopping parameter $\gamma_4$. We discuss possible ways to identify magnetoelectric effects by distinctive features in the optical conductivity.

Tables

Table 1. Table 1: Phenomenological parameters for bilayer graphene employed for numerical calculations in this paper adopted from Refs.  Winkler and Zülicke, 2015 ; Konschuh  et al. , 2012 .

intralayer nearest-neighbor hopping	${\gamma }_0$	=	$3.0\mathrm{eV}$
intralayer second-nearest-neighbor hopping	${\gamma }_0^{\prime }$	=	$0.22\mathrm{eV}$
dimer-dimer hopping	${\gamma }_1$	=	$0.32\mathrm{eV}$
non-dimer-dimer hopping	${\gamma }_3$	=	$0.25\mathrm{eV}$
non-dimer-non-dimer hopping	${\gamma }_4$	=	$0.14\mathrm{eV}$
interlayer distance	$d$	=	$0.335\mathrm{nm}$
lattice constant	$a$	=	$0.245\mathrm{nm}$

Equations

$${\alpha }_{ij}=$$

$${\alpha }_{\theta }=$$

$${\mathcal{L}}_{\text{ME}}=-e{n}_{\mathrm{v}}({\xi }_{\parallel }{\mathcal{E}}_{\parallel }\cdot {\mathbf{B}}_{\parallel }+{\xi }_z{\mathcal{E}}_z\cdot {\mathbf{B}}_z)\delta (z),$$

$${\mathbf{a}}_1=$$

$${\mathbf{a}}_2=$$

$${\mathbf{a}}_3=$$

$${\mathbf{\tau }}_1^{(j)}=$$

$${\mathbf{\tau }}_2^{(j)}=$$

$${\mathbf{b}}_1=$$

$${\mathbf{b}}_2=$$

$${\mathbf{b}}_3=$$

$$\mathbf{K}=$$

$${\psi }_{n,\mathbf{k}}(\mathbf{r})=$$

$$⟨{\psi }_m⟩\mathcal{H}{\psi }_n$$

$$\frac{1}{N}\sum _{l,j}\exp \left[i\mathbf{k}\cdot ({\mathbf{R}}_l-{\mathbf{R}}_j)-\frac{ie}{\hslash }{\int }_{{\mathbf{R}}_l}^{{\mathbf{R}}_j}\mathrm{d}\mathbf{r}\cdot \mathbf{A}(\mathbf{r})\right]$$

$$\times {\int }_{{\mathbb{R}}^3}\mathrm{d}{r}^3{\phi }_m^{\ast }(\mathbf{r}-{\mathbf{R}}_j)\mathcal{H}{\phi }_n(\mathbf{r}-{\mathbf{R}}_l),$$

$$\mathcal{H}$$

$$f_1(\mathbf{k})$$

$$f_1^{\pm }(\mathbf{k})$$

$$f_2^{\pm }(\mathbf{k})$$

$$\mathcal{H}(\pm \mathbf{K}+\mathbf{k})$$

$${\kappa }_0=\pm k_{\pm },{\kappa }_1=\pm (k_{\pm }\pm i{b}_{\pm }),{\kappa }_2=\pm (k_{\pm }\mp i{b}_{\pm }),b_{\pm }=\frac{ed}{2\hslash }(B_x\pm i{B}_y),$$

$${\tilde{\gamma }}_0=\frac{\sqrt{3}a}{2}{\gamma }_0,{\tilde{\gamma }}_0^{\prime }=\frac{3a^2}{4}{\gamma }_0^{\prime },{\tilde{\gamma }}_{31}=\frac{\sqrt{3}a}{2}{\gamma }_3,{\tilde{\gamma }}_{32}=\frac{a^2}{8}{\gamma }_3,{\tilde{\gamma }}_4=\frac{\sqrt{3}a}{2}{\gamma }_4.$$

$${\mathcal{H}}_{\text{eff}}=$$

$$+\left\{\frac{{\tilde{\gamma }}_0^2-{\gamma }_1{\tilde{\gamma }}_{32}+{\tilde{\gamma }}_4^2}{{\gamma }_1}(k_y^2-k_x^2){\tau }_0-\left[{\tilde{\gamma }}_{31}{k}_x+e\frac{{\tilde{\gamma }}_0{\tilde{\gamma }}_4}{{\gamma }_1^2}({\mathcal{E}}_x{b}_x-{\mathcal{E}}_y{b}_y)\right]{\tau }_z\right\}\otimes {\sigma }_x$$

$$+\left\{\left[-{\tilde{\gamma }}_{31}{k}_y+e\frac{{\tilde{\gamma }}_0{\tilde{\gamma }}_4}{{\gamma }_1^2}({\mathcal{E}}_x{b}_y+{\mathcal{E}}_y{b}_x)\right]{\tau }_0+2\frac{{\tilde{\gamma }}_0^2-{\gamma }_1{\tilde{\gamma }}_{32}+{\tilde{\gamma }}_4^2}{{\gamma }_1}{k}_x{k}_y{\tau }_z\right\}\otimes {\sigma }_y$$

$$+\left\{2\frac{{\tilde{\gamma }}_0{\tilde{\gamma }}_4+{\tilde{\gamma }}_0^{\prime }{\gamma }_1}{{\gamma }_1}(b_y{k}_x-b_x{k}_y){\tau }_0+e\frac{{\tilde{\gamma }}_0^2+{\tilde{\gamma }}_4^2}{{\gamma }_1^2}({\mathcal{E}}_y{k}_x-{\mathcal{E}}_x{k}_y){\tau }_z\right\}\otimes {\sigma }_z.$$

$${\mathcal{H}}_{\text{ax}}=$$

$${\xi }_{\parallel }=$$

$${\xi }_{\parallel }\approx$$

$$\theta_{\mathrm{BLG}}\equiv\frac{2\pi h}{e}\,\frac{2\xi_{\parallel}+\xi_{z}}{3d}\,n_{\mathrm{v}}\approx 6.7\times 10^{-3}\,n_{\mathrm{v}}\big{[}10^{10}\mathrm{cm}^{-2}\big{]}\,\,.$$

$${\tilde{\Delta }}^{\pm }$$

$${\sigma }_{\mu \nu }(\omega )=$$

$$\text{Re}\left[{\sigma }_{\mu \nu }^{\pm }(\omega )\right]=$$

Full text

I Introduction

Although having been extensively studied for more than half a century O’Dell (1970), the magnetoelectric (ME) effect has recently returned to the spotlight. The renewed attention is motivated by the discovery of new material systems. The traditionally predominant focus on band insulators with intrinsically broken space-inversion and time-reversal symmetries, known as multiferroics Fiebig (2005); Ramesh and Spaldin (2007); Spaldin (2017); Fiebig et al. (2016), has expanded to topological insulators Qi et al. (2008); Essin et al. (2009); Rosenberg and Franz (2010); Qi and Zhang (2011), Weyl semimetals Grushin (2012); Zyuzin et al. (2012); Zyuzin and Burkov (2012); Vazifeh and Franz (2013), and very recently to metals Ma and Pesin (2015); Varjas et al. (2016); Zhong et al. (2016); Thöle and Spaldin (2018). The latest finding of bilayer graphene being a ME medium further extends the list by adding an example from the rapidly expanding class of two-dimensional materials that have their own unique ME properties Zhang et al. (2011); Zülicke and Winkler (2014); Winkler and Zülicke (2015); Zülicke and Winkler (2017).

A particularly appealing feature of ME media is the possibility to manipulate magnetic properties in a solid by electric fields and vice versa. This enables us to engineer device architectures with novel functionalities not achievable with other materials Ramesh and Spaldin (2007); Fiebig et al. (2016); Spaldin (2017). It also establishes an inspiring connection between ideas and methods from condensed-matter physics, high-energy physics, and even cosmology Franz (2008); Spaldin (2017). In a nutshell, the ME effect becomes manifest in a mixing between electric and magnetic fields, $\bm{\mathcal{E}}$ and $\mathbf{B}$, in the expansion of the free energy $\mathcal{F}$. In leading-order, this coupling is linear in both fields and described by the ME tensor $\alpha$, whose components are defined as Fiebig (2005); Spaldin et al. (2008)

[TABLE]

It is common practice to decompose the ME tensor as Essin et al. (2010); Spaldin et al. (2008)

[TABLE]

with a pseudo-scalar $\alpha_{\theta}$ and a traceless tensor $\widetilde{\alpha}_{ij}$. The latter can be further split into a traceless symmetric part $\widetilde{\alpha}^{\text{S}}_{ij}$ and an antisymmetric part $\widetilde{\alpha}^{\text{A}}_{ij}$. These three terms are associated with the ME monopole, quadrupole and toroidal moments, respectively Spaldin et al. (2008); Thöle and Spaldin (2018).

Of particular interest is the pseudo-scalar part,

[TABLE]

with $e>0$ denoting the elementary charge, $h$ the Planck constant, and the dimensionless parameter $\theta$ that can be related to the axion field from astroparticle physics Franz (2008). It yields an isotropic coupling of electric and magnetic fields and can thus be interpreted as a condensed-matter realization of the axion electrodynamics Wilczek (1987); Hehl et al. (2008); Visinelli (2013). The latter is characterized by modified Maxwell equations that emerge from adding the term $\mathcal{L}_{\text{ax}}=\alpha_{\theta}\;\bm{\mathcal{E}}\cdot\mathbf{B}$ to Maxwell’s Lagrangian of classical electromagnetism. The resulting modifications only lead to physical effects if the axion field $\theta$ varies in space or time. A spatial variation of $\theta$ arises naturally by the presence of interfaces of ME and non-ME materials where, as a consequence, an anomalous Hall conductivity appears Qi et al. (2008); Essin et al. (2009). Also, the axion field $\theta$ possesses two fundamental properties: (i) it is invariant under the shift $\theta\rightarrow\theta+2\pi$; and (ii) due to the distinct transformation properties of $\bm{\mathcal{E}}$ and $\mathbf{B}$ under time reversal and spatial inversion, the axion field has to be odd with respect to both symmetry operations. In traditional ME media Fiebig (2005); Ramesh and Spaldin (2007); Spaldin (2017); Fiebig et al. (2016), the latter property is realized by the coexistence of an intrinsic ferromagnetic and ferroelectric order. In contrast, this precondition is circumvented in topological insulators without broken time-reversal and inversion symmetries by the equivalence of $\theta=\pm\pi$ due to property (i). Yet, this constrains the axion field to only appear in a quantized form, i.e., $\theta\in\{0,\pi\}$ Qi and Zhang (2011).

It has recently been established that the ME effect is also present in bilayer graphene and shows intriguing properties Zülicke and Winkler (2014); Winkler and Zülicke (2015); Zülicke and Winkler (2017). Here, the corresponding ME Lagrangian is anisotropic and may be approximated as Zülicke and Winkler (2014)

[TABLE]

where the delta distribution $\delta(z)$ locates the bilayer to be in the $x$-$y$ plane. It involves the valley-isospin density $n_{\mathrm{v}}$, i.e., the difference of electron densities in the two ($\mathbf{K}$ and $\mathbf{K^{\prime}}$) valleys. On the one hand, a finite density imbalance is generated by the axion-like ME coupling, as the latter induces a valley-contrasting potential shift Xiao et al. (2007); Nakamura et al. (2009). At the same time, this density dependence makes the strength of the ME response tunable and establishes a link to valleytronics Schaibley et al. (2016); Sui et al. (2015); Shimazaki et al. (2015), which is not the case in other known ME media. On the other hand, the special transformation properties of the valleys enable the presence of the ME effect even though time reversal and spatial inversion are symmetries of the crystal lattice. Since this observation is general and applies to any multi-valley system Novoselov et al. (2016) with similar symmetries, this indicates that there is another class of ME active materials with bilayer graphene being the first of its kind to be discovered. Aside from this, the strength of the axion field is determined by the material-dependent parameters, where the two-dimensional sheet geometry suggests that there should be generally distinct in-plane (proportional to $\xi_{\parallel}$) and out-of-plane (proportional to $\xi_{z}$) contributions. The parameter $\xi_{z}$ has been evaluated in Ref. Winkler and Zülicke, 2015, but the value of $\xi_{\parallel}$ has until now been unknown.

In this paper, we fill the knowledge gap about the ME response of a Bernal-stacked graphene bilayer in the presence of in-plane homogeneous electric and magnetic fields. Using a tight-binding description and applying quasi-degenerate perturbation theory, we analytically derive an effective low-energy Hamiltonian that comprises all relevant in-plane ME couplings and exhibits the ME equivalence Winkler and Zülicke (2015). This allows us to express the prefactor $\xi_{\parallel}$ in terms of tight-binding parameters. Inserting numerical values, $\xi_{\parallel}$ turns out to be approximately twice as large as $\xi_{z}$ and of the same sign. Additionally, small anisotropic ME terms induced by the electron-hole-symmetry breaking hopping $\gamma_{4}$ are found. To establish a connection to experiment, we discuss the impact of the ME couplings on features exhibited in the optical conductivity. For this purpose, the given system configuration is particularly suitable as, in presence of in-plane fields, the system remains metallic and the axionic response can be more prominent than for perpendicular fields. Also, complications arising from Landau quantization Zhang et al. (2011) can be avoided if the magnetic field is in-plane. We explicitly demonstrate that, due to the axionic term, the minimum optical absorption frequencies become valley-dependent for a non-vanishing chemical potential. Apart from this, the corrections arising from small anisotropic ME couplings lead to a broadening of the absorption peak at zero chemical potential. We treat the response of the graphene bilayer to the static in-plane electric field by invoking the drift-induced Fermi-sea displacement when calculating the optical conductivity.

This paper is structured as follows. In the following section, the general definitions concerning the crystal lattice and tight-binding model of bilayer graphene are briefly reviewed. In Sec. III, we first derive a tight-binding Hamiltonian describing our system of interest in the presence of an in-plane magnetic field by taking into account the arising Peierls phases. Including in-plane electric fields, in the next step, we compute an effective two-band Hamiltonian for the low-energy regime that contains all relevant ME couplings. In Sec. IV, we employ these results to study the impact on the optical conductivity, which turns out to exhibit distinctive features arising from the ME response. Electronic-structure parameters used for numerical calculations in this paper are listed in Table 1.

II Basic theory for the electronic structure of bilayer graphene

II.1 Crystal structure

The crystal lattice of a Bernal-stacked graphene bilayer is defined in accordance with Ref. Winkler and Zülicke, 2015 as illustrated in Figs. 1 and 2. The bilayer is composed of two coupled graphene monolayers which are characterized in real space by the primitive lattice vectors

[TABLE]

with the lattice constant $a$. Note that the distance of two adjacent carbon atoms within each layer is $a/\sqrt{3}$. Each of the coupled monolayers consists of two sublattices, which we define as ($A,B$) for the top and ($A^{\prime},B^{\prime}$) in the bottom layer, forming a hexagonal lattice. In the Bernal-stacked form the atoms are arranged such that the sublattices $A$ and $A^{\prime}$ lie on top of each other, i.e., they are connected by the vector

[TABLE]

where $d$ is the inter-layer distance. In contrast, the other sublattices ($B,B^{\prime}$) are displaced such that the corresponding atom is normal to the center of each hexagon of the other layer. In other words, the top layer can be generated by shifting the bottom layer by the vector $\mathbf{a}_{3}$ followed by a reflection at the $x$-$z$ plane. The point group of the bilayer graphene is $D_{3d}$ McCann and Koshino (2013). The atomic sites ($A,A^{\prime}$) are referred to as dimer and the sites ($B,B^{\prime}$) as non-dimer sites.

The intra-layer nearest-neighbor and second-nearest-neighbor vectors $\bm{\tau}_{1}$ and $\bm{\tau}_{2}$, with respect to, e.g., sublattice $A$, can be written as

[TABLE]

Here, $\mathbf{\tau}_{1}^{(3)}=a\mathbf{\hat{y}}/\sqrt{3}$ and $\mathcal{R}(\phi)$ denotes a rotation about the $z$-axis by the angle $\phi$.

In $\mathbf{k}$-space, the lattice retains its hexagonal shape but is rotated by $\pi/2$ about the $z$-axis. The according primitive reciprocal lattice vectors read as

[TABLE]

The two sublattices give rise to two kinds of inequivalent corner points $\mathbf{K}$ and $\mathbf{K^{\prime}}\equiv-\mathbf{K}$, where

[TABLE]

These corner points or valleys are of fundamental interest as the bandgap is minimal there or vanishes.

II.2 Tight-binding description

The tight-binding model has been used to study the electronic bandstructure for bilayer graphene in many different contexts. For a comprehensive review, see e.g. Refs. Castro Neto et al., 2009; McCann and Koshino, 2013; Rozhkov et al., 2016.

Within the tight-binding approach, the eigenfunctions are linear combinations of the Bloch functions

[TABLE]

where $\mathbf{k}$ denotes the (purely in-plane) wave vector of charge carriers in bilayer graphene, $\varphi_{n}(\mathbf{r}-\mathbf{R}_{l})$ is the $n$th atomic orbital at the lattice site $\mathbf{R}_{l}$, and $N$ is the total number of lattice sites McCann and Koshino (2013). Small corrections that arise from the non-orthogonality of the Bloch functions shall be neglected in the following. As it is the anti-bonding $\pi$ bonds derived from $p_{z}$-orbitals that are relevant for the electronic transport, we consider one $p_{z}$-orbital for each of the four sites within the unit cell, i.e., $n\in\{A,B,A^{\prime},B^{\prime}\}$. With this, we represent the tight-binding Hamiltonian $\mathcal{H}$ in the basis $\{\ket{\psi_{A}},\ket{\psi_{B}},\ket{\psi_{A^{\prime}}},\ket{\psi_{B^{\prime}}}\}$, where we use this very ordering. Including a magnetic field $\mathbf{B}=\nabla\times\mathbf{A}(\mathbf{r})$ with the vector potential $\mathbf{A}(\mathbf{r})$, the Bloch function picks up a Peierls phase Kheirabadi et al. (2016, 2018). A general matrix element thus becomes

[TABLE]

where the vector potential yields a phase given by a line integral between the different lattice sites.

The on-site energies are defined as $\epsilon_{n}=\braket{\psi_{n}}{\mathcal{H}}{\psi_{n}}$, where we assume for simplicity $\epsilon_{A}=\epsilon_{A^{\prime}}=\epsilon_{B}=\epsilon_{B^{\prime}}$, i.e., energy differences between the distinct sublattices as well as dimer and non-dimer sites are neglected. In this paper, we account for the couplings between nearest and second-nearest neighbors that are characterized by the hopping integrals $\gamma_{0}$, $\gamma_{0}^{\prime}$, $\gamma_{1}$, $\gamma_{3}$, and $\gamma_{4}$, as schematically illustrated in Fig. 2. A more detailed definition is given in Appendix A, and numerical values of all parameters are listed in Table 1. Typically, the nearest-neighbor hopping parameter $\gamma_{0}$ (intra-layer) and $\gamma_{1}$ (inter-layer) constitute the dominant couplings. The $\gamma_{3}$ hopping leads to a trigonal warping of the bandstructure, and the parameters $\gamma_{4}$ and $\gamma_{0}^{\prime}$ break the electron-hole symmetry.

For vanishing magnetic fields, the Hamiltonian takes the form Winkler and Zülicke (2015)

[TABLE]

where

[TABLE]

and $f_{2}(\mathbf{k})=|f_{1}(\mathbf{k})|^{2}-3$. Expanding $f_{1}$ in the vicinity of the high-symmetry points $\mathbf{K}$ and $\mathbf{K^{\prime}}\equiv-\mathbf{K}$, $f_{1/2}(\mathbf{k})\mapsto f_{1/2}(\pm\mathbf{K}+\mathbf{k})\equiv f_{1/2}^{\pm}(\mathbf{k})$, gives

[TABLE]

where $k_{\pm}=k_{x}\pm ik_{y}$ and $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$.111Notice that, in Ref. Winkler and Zülicke, 2015, the prefactor of $k_{\pm}^{2}$ in the expression for $f_{1}^{\pm}(\mathbf{k})$ differs by a factor of $2$.

III Magnetoelectric coupling of in-plane electromagnetic fields

Hereafter, we employ the definitions and notations of the previous section and derive a model Hamiltonian that includes the ME couplings due to in-plane electric and magnetic fields.

III.1 Incorporating in-plane magnetic fields into the tight-binding

Hamiltonian

The situation of an in-plane magnetic field has been addressed recently by several authors Pershoguba and Yakovenko (2010); Roy and Yang (2013); Van der Donck et al. (2016); Kheirabadi et al. (2016, 2018). In Ref. Pershoguba and Yakovenko, 2010, the magnetic field was semiclassically included by adding Lorentz-force-induced momentum shifts, and only vertical interlayer hopping $\gamma_{1}$ was considered. This method was extended in Refs. Roy and Yang, 2013; Van der Donck et al., 2016 by taking into account the effect of trigonal warping. Here, we follow a more sophisticated approach Kheirabadi et al. (2016, 2018) that correctly accounts for the arising Peierls phases in the tight-binding model.

A homogeneous in-plane magnetic field $\mathbf{B}=B_{x}\mathbf{\hat{x}}+B_{y}\mathbf{\hat{y}}$ can be associated with a vector potential $\mathbf{A}=z(B_{y}\mathbf{\hat{x}}-B_{x}\mathbf{\hat{y}})$. Selecting this gauge, the translation symmetry is preserved within each layer. Further assuming a symmetric arrangement of the top and bottom layers at $z=\pm d/2$, the Peierls phase for the inter-layer couplings vanishes. The strong dimer-dimer coupling $\gamma_{1}$ opens a gap for the $(A,A^{\prime})$-like states. For low energies $|E|<\gamma_{1}$, the bands of interest are described by the $(B,B^{\prime})$-like states. Hence, it is common to project on the $(B,B^{\prime})$ subspace and perturbatively include the couplings to the other bands. Following this approach, we incorporate the effects of a small magnetic field in the Hamiltonian; that is, we retain terms up to first (second) order in the field or the wave vector on the off-diagonal (diagonal) part of the Hamiltonian. Since the energy dispersion of bilayer graphene turns out to be dominated by terms that are quadratic in the wave vector, we allow, in addition, terms quadratic in the wave vector off-diagonal in the $(B,B^{\prime})$ sector as we project on that subspace. Rearranging further the basis functions as $\{\frac{1}{\sqrt{2}}(\ket{\psi_{A}}+\ket{\psi_{A^{\prime}}}),\frac{1}{\sqrt{2}}(\ket{\psi_{A}}-\ket{\psi_{A^{\prime}}}),\ket{\psi_{B}},\ket{\psi_{B^{\prime}}}\}$ and neglecting the constant energy shift due to $f_{2}$, we obtain

[TABLE]

with

[TABLE]

and

[TABLE]

Disregarding the parabolic terms $\propto\tilde{\gamma}_{0}^{\prime},\tilde{\gamma}_{32}$, this result coincides with the Hamiltonian given in Ref. Kheirabadi et al., 2016 apart from a unitary transformation and the sign of the $\gamma_{3}$-terms. Without magnetic field, this expression corresponds to the Slonczewski-Weiss-McClure Hamiltonian McClure (1957); Slonczewski and Weiss (1958).

The energy dispersion obtained for a finite in-plane magnetic field is displayed in Fig. 3. For better visualization, we used an extraordinary high magnetic field of $1000\text{\,}\mathrm{T}$. In Refs. Pershoguba and Yakovenko, 2010; Roy and Yang, 2013; Van der Donck et al., 2016, it was found that a large in-plane magnetic field produces a change in topology of the band structure similar to the one appearing due to lateral strain Mucha-Kruczyński et al. (2011); Son et al. (2011); Mariani et al. (2012); He et al. (2014); Daboussi et al. (2014). More precisely, the parabolic low-energy dispersion splits into two Dirac cones, cf. Figs. 3(a) and (b), where the new Dirac points appear at wave vectors $\mathbf{k}=\pm(\bm{\hat{z}}\times\mathbf{b})$. However, our more detailed model includes additional hopping terms that further reduce the symmetry and result in a gapped spectrum [cf. Fig. 3(c)]. Due to the presence of trigonal warping, the precise dispersion near the charge-neutrality point is quite complex and depends sensitively on the system configuration.

III.2 Effective low-energy Hamiltonian describing in-plane

magnetoelectric couplings

An effective Hamiltonian that describes the low-energy excitations of bilayer graphene in absence of fields was first derived in Ref. McCann and Fal’ko, 2006. The result was extended in Ref. Winkler and Zülicke, 2015 to the situation where perpendicular, static, and homogenous electric and magnetic fields are present.

In this paper, we consider both static homogeneous electric and magnetic fields applied within the plane of the bilayer. The magnetic field is already included in our tight-binding model. To also account for an electric field $\bm{\mathcal{E}}=\mathcal{E}_{x}\mathbf{\hat{x}}+\mathcal{E}_{y}\mathbf{\hat{y}}$, we add the scalar potential ${V_{\mathcal{E}}(\mathbf{r})=e\,\bm{\mathcal{E}}\cdot\mathbf{r}}$ to the Hamiltonian in Eq. (20). Using quasi-degenerate perturbation theory (QPT) Bir and Pikus (1974); Winkler (2003), we project on the $(B,B^{\prime})$ subspace. For the partitioning of the Hamiltonian, we select the diagonal terms proportional to the magnetic field and all off-diagonal terms as a perturbation. To third order in perturbation theory, this procedure yields the effective two-band Hamiltonian $\mathcal{H}_{\text{eff}}+V_{\mathcal{E}}(\mathbf{r})$, where

[TABLE]

In this notation, the Pauli matrices $\sigma_{0,x,y,z}$ are associated with the sublattice-related pseudospin degree of freedom. On the other hand, $\tau_{0}$ and $\tau_{z}$ are Pauli matrices whose basis states represent the different valleys in the order $(\mathbf{K},\mathbf{K^{\prime}})$. Notably, the effective Hamiltonian at the $\mathbf{K^{\prime}}$ valley can be obtained from the Hamiltonian at the $\mathbf{K}$ valley (and vice versa) by a mirror reflection at the $yz$-plane, i.e., the polar vectors map as $k_{x}\rightarrow-k_{x}$ and $\mathcal{E}_{x}\rightarrow-\mathcal{E}_{x}$ and the axial (pseudo-)vectors as $B_{y}\rightarrow-B_{y}$. In above expression for $\mathcal{H}_{\text{eff}}$, we excluded terms that are of third order in the wave vector or of second order in the magnetic field. (For completeness, terms quadratic in the magnetic field are listed in Appendix B. Additional ME couplings that appear in fourth order QPT are given in Appendix C. To make the connection with results from Ref. Winkler and Zülicke, 2015, we discuss in Appendix D possible couplings between in-plane magnetic and out-of-plane electric fields.) The effective Hamiltonian Eq. (23) exhibits the ME equivalence; i.e., it is form-invariant with respect to interchanging corresponding Cartesian components of the electric and magnetic fields Winkler and Zülicke (2015).

In absence of electric fields, the magnetic-field-dependent terms in Eq. (23) appear only due to the small tunneling amplitudes $\gamma_{0}^{\prime}$ and $\gamma_{4}$ that break the electron-hole symmetry. Hence, to observe the change of the band structure topology Pershoguba and Yakovenko (2010); Roy and Yang (2013); Van der Donck et al. (2016) as shown in Fig. 3, we need to include corrections that are quadratic in the magnetic field, Eq. (41). This implies that these terms should be taken into account for large magnetic fields. On the other hand, realistic parameters require extraordinary large magnetic fields to observe this effect experimentally in bilayer graphene. For instance, the characteristic energy splitting $\Delta_{\text{B}}$ at the $\pm\mathbf{K}$ points can be estimated to be $\Delta_{\text{B}}{}=\frac{3}{8\gamma_{1}}\left(ade\gamma_{0}B/\hbar\right)^{2}$, which yields a splitting of $\Delta_{\text{B}}=($1.64\text{\times}{10}^{-7}\text{,}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{T}^{2}$)\times B^{2}$ for bilayer graphene. A comparison of the low-energy dispersions for a large magnetic field obtained within different models is provided in Fig. 4. Ignoring the effect of trigonal warping, the electron and hole branches have generally two minima with a gap of approximately $4(\tilde{\gamma}_{0}\tilde{\gamma}_{4}/\gamma_{1}+\tilde{\gamma}_{0}^{\prime})\,b^{2}$.222As the impact of trigonal warping is not small near the charge-neutrality point, this value largely underestimates the real gap. If both electric and magnetic fields are present, the ME and purely magnetic terms are competing with each other. One consequence is that, above a critical electric field strength, the electron and hole branches can have one single extremum each. Assuming $\bm{\mathcal{E}}\parallel\mathbf{B}$ and setting $\gamma_{3}=\gamma_{4}=\gamma_{0}^{\prime}=0$, we can estimate the transition to occur at $\|\bm{\mathcal{E}}\|/\|\mathbf{B}\|\approx\gamma_{1}d/(\sqrt{2}\hbar)$. Yet, we stress that the low-energy band structure has much richer features due to the trigonal warping and the gap is highly anisotropic. A more detailed discussion thereof goes beyond the scope of this work and shall be presented elsewhere.

Our main interest concerns the ME terms, which couple the electric and magnetic fields and induce a breaking of the valley degeneracy Xiao et al. (2007); Nakamura et al. (2009); Zülicke and Winkler (2014); Winkler and Zülicke (2015). As the leading contribution, we identify the axionic term

[TABLE]

with $\Delta_{\text{ax}}=e\,\xi_{\parallel}\,\bm{\mathcal{E}}\cdot\mathbf{B}$ involving the material parameter

[TABLE]

Using the values for tight-binding parameters given in Table 1, we estimate

[TABLE]

which has the same sign and is about twice as large as the prefactor for the axionic term involving perpendicular fields; $\xi_{z}\approx$6.0\text{\times}{10}^{-4}\text{,}\mathrm{n}\mathrm{m}\mathrm{/}\mathrm{T}$$ Zülicke and Winkler (2014). Axionic terms are of particular interest as they constitute a condensed-matter realization of axion electrodynamics. With the values for $\xi_{\parallel}$ and $\xi_{z}$ given above, in conjunction with expressions from Eqs. (3a) and (8) of Ref. Zülicke and Winkler, 2014, we find the magnitude of the axionic pseudoscalar in bilayer graphene as

[TABLE]

For comparison, the maximum value of the axionic pseudoscalar measured in the archetypal magnetoelectric Hehl et al. (2008) Cr2O3 is $\theta_{\mathrm{Cr}_{2}\mathrm{O}_{3}}^{(\mathrm{max)}}\approx 0.13$, while in topological insulators, a large pseudoscalar magnitude $\theta_{\text{TI}}=\pi$ is fixed by time-reversal symmetry Qi and Zhang (2011). Furthermore, the axionic terms in bilayer graphene generate an energy shift $\pm\Delta_{\text{ax}}$ of equal magnitude but opposite sign for pseudospin eigenstates in the two valleys $\pm\mathbf{K}$. This leads to a finite valley-isospin density $n_{\mathrm{v}}=n_{\textbf{K}}-n_{-\textbf{K}}$, that is, the difference of charge densities $n_{\pm\textbf{K}}$ in the distinct valleys Xiao et al. (2007); Zülicke and Winkler (2014); Winkler and Zülicke (2015). In the next section, we will show how this feature manifests itself in a valley-dependent minimal absorption frequency in the optical conductivity spectrum when the chemical potential is not at the charge-neutrality point.

Besides the axionic contributions, Eq. (23) contains other ME-coupling terms corresponding to anisotropic contributions of the traceless tensor $\widetilde{\alpha}_{ij}$ in Eq. (2). They are smaller than the uniaxial terms by a factor $\tilde{\gamma}_{4}/\tilde{\gamma}_{0}=\gamma_{4}/\gamma_{0}\approx$4.7\text{\times}{10}^{-2}\text{,}$$. In conjunction with the quadratic-in-magnetic-field corrections, these terms lead to an energy gap

[TABLE]

between the electron and hole branches in the two valleys $\pm\mathbf{K}$. The opposite sign of the purely ME contribution leads to a gap difference between the valleys, i.e., $|\tilde{\Delta}^{+}-\tilde{\Delta}^{-}|\approx 4\gamma_{4}|\Delta_{\text{ax}}|/\gamma_{0}$. As shown in the next section, this property turns out to generate a step-like structure in the optical conductivity when the chemical potential is at the charge-neutrality point.

IV Visibility of magnetoelectric coupling in the optical conductivity

In the remainder of this paper, we explore the possibility to detect the above-discussed ME couplings in bilayer graphene through a transport measurement. Concomitantly with causing the ME effects that are our primary interest, the presence of the static uniform in-plane electric field $\bm{\mathcal{E}}$ will also generate a stationary DC current that is associated with a shifted Fermi sea of charge carriers in bilayer graphene. We envision applying an additional small AC electric field $\delta\bm{\mathcal{E}}(t)$, which results in an AC contribution $\delta\mathbf{j}(t)$ to the current density. The tensor $\sigma_{\mu\nu}(\omega)$ of the frequency-dependent (optical) conductivity relates the Fourier-transformed AC quantities [current density $\delta\mathbf{j}(\omega)$ and electric field $\delta\bm{\mathcal{E}}(\omega)]$ via the linear relation ${\delta j_{\mu}(\omega)=\sigma_{\mu\nu}(\omega)\,\delta\mathcal{E}_{\nu}(\omega)}$. Fundamental properties of the electronic bandstructure can give rise to distinctive features in the optical conductivity, making the latter a prime tool for the study of unconventional materials including graphene. It has been used to study various kinds of systems such as single or multilayer graphene,Stauber et al. (2008); Mak et al. (2008); Koshino and Ando (2008, 2009) bilayer graphene with and without asymmetry gap Nicol and Carbotte (2008), and recently even taking into account a finite twist angle between the two layers. Stauber et al. (2018a, b) In our present case, the frequency dependence of $\sigma_{\mu\nu}(\omega)$ will not only be affected by the ME-effect-related manipulation of the bandstructure, but also by the non-equilibrium distribution of charge carriers within this bandstructure. In the following, we elucidate both effects in turn.

IV.1 Kubo formalism to calculate the frequency-dependent electric

conductivity

We employ the Kubo formula to calculate the conductivity tensor $\sigma_{\mu\nu}(\omega)$. The electric-field perturbation is considered to be spatially homogeneous, time-dependent, and parallel to the bilayer plane, i.e., $\delta\bm{\mathcal{E}}(t)=\delta\bm{\mathcal{E}}(\omega)\exp[-i(\omega+i\eta)t]$. Here we have included an infinitesimally small quantity $\eta\in\mathbbm{R}^{+}$. To avoid generating time-dependent contributions to $\Delta_{\mathrm{ax}}$, the oscillating electric field $\delta\bm{\mathcal{E}}(t)$ should be applied perpendicular to the static in-plane magnetic field $\mathbf{B}$. The conductivity tensor can be expressed in terms of the set of single-particle eigenstates $\{\ket{n}\}$ in the frequency domain as Kammermeier (2018); Wenk (2011)

[TABLE]

with the volume $\mathcal{V}$, the current operator $\textbf{J}=-e(\nabla_{\mathbf{k}}\mathcal{H})/\hbar$, and the single-particle eigenenergies $\epsilon_{i}$. Moreover, the function $f(\epsilon_{i})=\{1+\exp[\beta(\epsilon_{i}-\tilde{\mu})]\}^{-1}$, where $\beta=1/(k_{B}T)$, represents the Fermi-Dirac distribution with the Boltzmann constant $k_{B}$, the temperature $T$, and the chemical potential $\tilde{\mu}$. We are interested in the dissipative part which is given by the real part of the conductivity tensor $\text{Re}\left[\sigma_{\mu\nu}(\omega)\right]$. Here, we can further distinguish two terms: (i) The intra-band ($n=l$) contribution, which determines the DC Drude conductivity, and (ii) the inter-band ($n\neq l$) contribution, which determines the optical absorption at finite frequencies.

To illustrate the emergence of features in the optical conductivity associated with magnetoelectricity, we now focus on the inter-band contribution to the optical conductivity for bilayer graphene in the presence of in-plane ME couplings in the clean limit. The low-energy bandstructure of this system is determined by the single-particle Hamiltonian $\mathcal{H}_{\mathrm{eff}}$ displayed in Eq. (23). Its eigenstates are of the form $\ket{n}=\ket{\mathbf{k}}\otimes\ket{\pm}\otimes\ket{\lambda}_{\mathbf{k}}$ where $\mathbf{k}$ denotes the wave vector, $\pm$ the valley-index, and $\lambda$ distinguishes the different electron and hole branches. Since the Hamiltonian $\mathcal{H}_{\mathrm{eff}}$ is diagonal with respect to $\ket{\mathbf{k}}$ and $\ket{\pm}$, the inter-band optical conductivity of charge carriers from the individual valleys $\pm\mathbf{K}$ becomes Bernád et al. (2010)

[TABLE]

where $\mathcal{H}_{\mathrm{eff}}^{\pm}=\braket{\pm}{\mathcal{H}_{\mathrm{eff}}}{\pm}$, $\sigma_{0}=2e^{2}/h$, and we introduced an additional factor $2$ to account for spin degeneracy. Here and in the following, the $\mathbf{k}$-integration is restricted to the domain around the valleys as otherwise spurious solutions may occur due to the small wave vector expansion of the Hamiltonian $\mathcal{H}_{\text{eff}}$. The total optical conductivity of the system is the sum of contributions from the individual valleys, i.e., $\text{Re}\left[\sigma_{\mu\nu}(\omega)\right]=\sum_{\pm}\text{Re}\left[\sigma_{\mu\nu}^{\pm}(\omega)\right]$.

In contrast to the usual situation, the distribution of charge carriers in the unperturbed state for our case of interest is a uniformly shifted Fermi sea. See Appendix E for a detailed discussion. Accounting for the stationary current-carrying state generated by the static in-plane electric field finally amounts to using a $\mathbf{k}$-dependent chemical potential [cf. Eq. (48)]

[TABLE]

in the expression Eq. (30) for the optical conductivity, where $\tau_{tr}$ is the intra-valley transport relaxation time.

To disentangle the non-equilibrium kinetic effect of the in-plane electric field from features associated with ME couplings, we present below results obtained for $\sigma_{\mu\nu}(\omega)$ both with and without the $\mathbf{k}$-dependent correction to $\tilde{\mu}$ included in the formula Eq. (30). For full consistency, life-time broadening on the scale of $\tau_{tr}$ should also be included in the calculation of $\sigma_{\mu\nu}(\omega)$, and the latter’s salient features need to be sufficiently separated from the intra-band (broadened-Drude-peak) contribution to the optical conductivity to enable experimental observation.

IV.2 Ramification of in-plane magnetoelectric couplings for the optical

conductivity: Discussion

In Section III.2, we identified two major physical implications arising from ME effects involving in-plane electric and magnetic fields: the valley-asymmetric axionic energy shift Eq. (24) and the valley-dependent field-tunable gap Eq. (28). We now discuss the features in the optical conductivity associated with each of these effects.

To start with, we consider the case where ${\gamma_{4}=\gamma_{0}^{\prime}=0}$ and, hence, the in-plane electric field causes only the axionic shift. [The same holds for the in-plane magnetic field if the parabolic terms Eq. (41) are disregarded.] Neglecting the $\mathbf{k}$-dependent correction to the chemical potential for now, the expression Eq. (30) for the conductivity tensor can be factorized into two parts,333This factorization is generally valid Bernád et al. (2010) for any electron-hole-symmetric bandstructure.

[TABLE]

with valley-dependent chemical potentials $\tilde{\mu}_{\pm}=\tilde{\mu}\mp\Delta_{\text{ax}}$ appearing only in the function

[TABLE]

and the part

[TABLE]

being fully determined by the bandstructure and independent of the valley-index which is therefore suppressed. In Fig. 5 it is illustrated that, depending on the value of the chemical potential $\tilde{\mu}$, the minimum transition frequencies $\omega_{\pm}=2|\tilde{\mu}\mp\Delta_{\text{ax}}|/\hbar$ for inter-band transitions at a fixed wave vector are in general valley-dependent. For $\tilde{\mu}=0$, i.e., with the chemical potential at the charge-neutrality point, the minimal frequencies are equally large and determined by the axionic shift. For $\tilde{\mu}\neq 0$, both peaks separate, and the difference $|\omega_{+}-\omega_{-}|$ becomes maximal as soon as the chemical potential exceeds the axionic energy shift, i.e., $|\tilde{\mu}|\geq|\Delta_{\text{ax}}|$. As this property becomes manifest in the optical absorption spectrum, the magnitude of the axionic term is, in principle, accessible in experiment.

We now discuss the situation when the $\gamma_{4}$-terms are included. These result in a breaking of the electron-hole symmetry and a valley-contrasting energy gap due to additional ME couplings (cf. Sec. III.2). This leads to valley-dependent corrections to the minimum transition frequencies, which can be approximated (by neglecting $\gamma_{3}$, $\gamma_{0}^{\prime}$, as well as quadratic terms in the electric or magnetic field) as $\omega_{\pm}\approx{2(1+2\zeta\gamma_{4}/\gamma_{0})|\tilde{\mu}\mp\Delta_{\text{ax}}|/\hbar}$ where $\zeta=1$ if $\pm\Delta_{\text{ax}}>\tilde{\mu}$ and $\zeta=-1$ otherwise. The effect becomes particularly pronounced for a vanishing chemical potential where a broadening of the minimal-frequency absorption peak of $|\omega_{+}-\omega_{-}|=8\gamma_{4}|\Delta_{\text{ax}}|/(\gamma_{0}\hbar)$ occurs.

The case where we neglect electron-hole asymmetry and trigonal warping by setting ${\gamma_{3}=\gamma_{4}=\gamma_{0}^{\prime}=0}$ facilitates further analytical treatment. Including the parabolic terms in the magnetic field, Eq. (41), the energy eigenvalues read as

[TABLE]

where the upper and lower sign corresponds to the electron ($e$) and hole ($h$) branches, respectively. Here, we introduced $\phi_{\mathbf{a}}$ to indicate the polar angle that the in-plane vector $\mathbf{a}$ encloses with the $x$ axis, and we omitted the constant axionic shift as this is absorbed into the valley-dependent chemical potentials $\tilde{\mu}_{\pm}$. Focusing on the longitudinal conductivity correction (selecting $\sigma_{xx}$ without loss of generality), we can simplify Eq. (34) to

[TABLE]

The distinctive features appearing in the optical conductivity due to the axionic shift are illustrated in Fig. 6, both for the idealized case ${\gamma_{3}=\gamma_{4}=\gamma_{0}^{\prime}=0}$ and for the real system at zero temperature. For comparison, the optical conductivity spectrum without static electric and magnetic fields is shown in Fig. 6(a). In Figs. 6(b) and (c) we assume $\mathcal{E}_{x}=0$, $B_{x}=0$, $\mathcal{E}_{y}=$0.05\text{,}\mathrm{V}\mathrm{/}\mathrm{n}\mathrm{m} and $B_{y}=$100\text{\,}\mathrm{T}. For this situation, the axion shift is $\Delta_{\text{ax}}=$5.0\text{,}\mathrm{m}\mathrm{e}\mathrm{V}$$, and the electron and hole branches have a single extremum (cf. Sec. III.2). Differently colored curves represent results obtained for different chemical potentials $\tilde{\mu}$. The vertical red solid grid lines mark the minimum transition frequency $2\Delta_{\text{ax}}/\hbar$ in the case of $\tilde{\mu}=0$ and electron-hole symmetry. In Fig. 6(b) the vertical red dashed grid lines depict the small broadening of the minimum frequency transition peak of $8\gamma_{4}|\Delta_{\text{ax}}|/(\gamma_{0}\hbar)\approx$1.8\text{,}\mathrm{m}\mathrm{e}\mathrm{V}$/\hbar$ due to the electron-hole asymmetry and the ME couplings $\propto\gamma_{4}$. For increasing chemical potential the minimum absorption peaks shift to higher or lower frequencies depending on the valley index. As soon as the chemical potential exceeds the axion shift, both peaks move simultaneously to higher frequencies while retaining the constant difference of approximately $2\Delta_{\text{ax}}/\hbar$ (cf. the orange curve for $\tilde{\mu}=20\times 10^{-3}\gamma_{1}$). Notice also that for both $\tilde{\mu}=13\times 10^{-3}\gamma_{1}$ (light green curve) and $\tilde{\mu}=17\times 10^{-3}\gamma_{1}$ (yellow curve) the chemical potential lies within the gap of one valley and therefore the respective absorption peaks are identical. In Fig. 6(c) the conductivity spectrum is enlarged for the lowest three chemical potentials. Dashed curves show results obtained from the above-discussed simplified theory where ${\gamma_{3}=\gamma_{4}=\gamma_{0}^{\prime}=0}$, yielding Eq. (36). In this approximation, the minimum transition peaks are sharp and occur precisely at $\omega_{\pm}=2|\tilde{\mu}\mp\Delta_{\text{ax}}|/\hbar$ as emphasized by the vertical gray grid lines. Lastly, we illustrate with the horizontal gray grid lines in both plots, that asymptotically the conductivity recovers the universal value of $\text{Re}[\sigma_{xx}]=\frac{\pi}{2}\sigma_{0}$ (or $\text{Re}[\sigma_{xx}^{\pm}]=\frac{\pi}{4}\sigma_{0}$ for either valley), which is in agreement with the established high-frequency limit without fields.Nicol and Carbotte (2008)

As discussed briefly above and in greater detail in Appendix E, we can account for the corrections arising from the carrier drift due to the static electric field $\bm{\mathcal{E}}$ by considering an effectively $\mathbf{k}$-dependent chemical potential $\tilde{\mu}(\mathbf{k})$, as given in Eq. (31). The resulting modifications of the optical-conductivity spectrum are displayed in Fig. 7(a) for various relaxation times and $\tilde{\mu}(0)=0.1\gamma_{1}$. The vertical grid lines correspond to the the valley-dependent minimal absorption frequencies for neglected drift current. In Fig.7(b) we plot the effective chemical potential $\tilde{\mu}(\mathbf{k})$ with respect to the energy dispersion. It is chosen large enough that the Fermi energy lies above the gap in both valleys and the difference of the valley-dependent minimal transition frequencies is maximal. We see that for an increasing drift, the associated energy window obscures the ME features, albeit sharp features associated with these remain. Further inclusion of life-time broadening into the expression Eq. (30) for the optical conductivity at the assumed scale of $\tau_{tr}$ will likely be deleterious to any such features. Yet, the issue of the drift current can be circumvented if the magnetic field plays the dominant role in the valley-contrasting energy shift $\Delta_{\text{ax}}$ (cf. Appendix E).

While conclusive experimental detection of ME effects in pristine bilayer graphene may be a challenge because of their smallness, the presented illustration for their fingerprints in the optical conductivity may serve as a useful guide to inform experimental studies focused on other two-dimensional materials having larger ME couplings. Future work could also explore the potential use of generally large proximity-induced exchange fields Wei et al. (2016) and strain-generated pseudo-magnetic fields Jiang et al. (2017) for boosting magnetoelectric effects in two-dimensional materials.

V Conclusion and Outlook

We have investigated the in-plane ME couplings in bilayer graphene and how their properties give rise to distinctive features in the optical conductivity. We developed a tight-binding description that correctly accounts for the Peierls phases induced by a magnetic field parallel to the bilayer plane. Taking into account second-nearest neighbor hoppings, in particular, the electron-hole symmetry breaking contributions, it is shown that the spectrum is generally gapped, which was not observed within the simplified models that were employed previously.Pershoguba and Yakovenko (2010); Roy and Yang (2013); Van der Donck et al. (2016) In the next step, we included the effect of an in-plane electric field and derived an effective two-band Hamiltonian that comprises all relevant ME couplings and expresses them in terms of tight-binding parameters. We identify an axion-like pseudoscalar contribution to the ME tensor which has the same sign and about twice the magnitude of the previously determined out-of-plane contribution. In addition, small corrections due to the small skew hopping amplitude $\gamma_{4}$ are found, which correspond to the ME quadrupole and toroidal moment. The Hamiltonian also exhibits the equivalence of the electric and magnetic fields that was previously predicted by means of group-theoretical methods.Winkler and Zülicke (2015) A more realistic calculation should also involve Coulomb-interaction effects. For instance, Refs. Kusminskiy et al., 2009; Cheianov et al., 2012 showed how dielectric screening causes band-structure renormalizations that can also be expected to influence ME couplings.

Lastly, we use the effective two-band Hamiltonian to study the impact of the ME terms on the low-energy optical conductivity. Although in an ideal situation each of the ME contribution yields clear features in the frequency dependence of the conductivity, the results are obscured in bilayer graphene. This is a consequence of the displaced Fermi contour due to a finite drift current and the small magnitude of the axionic ME coupling in this material. Nonetheless, these observations remain pertinent for systems with similar symmetries yet more pronounced ME couplings or in the presence of large magnetic fields that allow to reduce the drift-current-inducing electric fields while maintaining a significant axionic energy shift.

VI Acknowledgement

This work was supported by the Marsden Fund Council from Government funding (contract no. VUW1713), managed by the Royal Society Te Apārangi. P.W. acknowledges funding from the German Research Foundation (DFG) via grant SFB 1277 and project 336985961.

Appendix A Tight-binding parameters

In this paper, we take into account the (intra-layer as well as inter-layer) tunneling from each lattice site to its nearest and second-nearest neighbor as schematically depicted in Fig. 2. These couplings are characterized by the hopping integrals

[TABLE]

and the second-nearest-neighbor intra-layer coupling $\gamma_{0}^{\prime}={}\braket{\varphi_{A}}{\mathcal{H}}{\varphi_{A}}=\braket{\varphi_{B}}{\mathcal{H}}{\varphi_{B}}=\braket{\varphi_{A^{\prime}}}{\mathcal{H}}{\varphi_{A^{\prime}}}=\braket{\varphi_{B^{\prime}}}{\mathcal{H}}{\varphi_{B^{\prime}}}$. In literature there appear various definitions for the hopping parameters with respect to their signs. We selected our employed definitions consistent with Ref. Winkler and Zülicke, 2015. Numerical values for the hopping integrals as well as other system parameters that are used in this work are listed in Table 1.

Appendix B Parabolic correction in the magnetic field

For large magnetic fields, the terms that are of second order in the magnetic field may become important as they can lead to a change of the bandstructure topology (cf. Sec. III). Neglecting wave-vector dependent terms, the leading-order contribution reads as

[TABLE]

where we dropped a global shift $\propto\tau_{0}\otimes\sigma_{0}$. Notably, these terms are already obtained in second-order QPT.

Appendix C Higher-order magnetoelectric couplings

In fourth-order QPT we obtain higher order ME couplings due to trigonal warping ($\gamma_{3}$) which are linear in the electric but parabolic in the magnetic field

[TABLE]

or linear in both fields and the wave vector

[TABLE]

where we neglected in the latter equation subordinate terms $\propto\gamma_{3}\gamma_{4}$.

Appendix D Leading magnetoelectric couplings in presence of perpendicular electric fields

To additionally account for perpendicular electric fields which open up a pseudospin gap,McCann (2006); Ohta et al. (2006); Zhang et al. (2009) we may include in Eq. (20) the potentialWinkler and Zülicke (2015); Konschuh et al. (2012)

[TABLE]

where the ab initio calculations in Ref. Konschuh et al., 2012 yield $\varepsilon_{12}=\varepsilon_{33}=0.048e\,\text{nm}$. From symmetry considerations, this potential should give rise to ME couplings $\propto(\mathcal{E}_{z}B_{y}\sigma_{x}-\mathcal{E}_{z}B_{x}\sigma_{y})$.Winkler and Zülicke (2015) Within our model, however, these terms do not appear, meaning that their prefactor vanishes exactly. Instead, the first non-vanishing terms are obtained in third-order QPT and correspond to higher ME couplings analogously to the in-plane couplings in Appendix C. These terms read as

[TABLE]

Hence, the mixing of electric and magnetic fields perpendicular to each other causes only small corrections to the axionic terms.

Appendix E Preconditions for the linear response model

In this section, we specify the requirements that ensure the validity of the linear response model for the optical conductivity in Sec. IV.

Firstly, the static electric field $\bm{\mathcal{E}}$ induces a finite parallel drift of the carriers and the system is out of equilibrium. This yields a modification of the Fermi-Dirac distribution asRammer (2004)

[TABLE]

where $\tau_{tr}$ denotes the intra-valley transport relaxation time of the carriers and the upper and lower sign corresponds to the electron ($e$) and hole ($h$) branches, respectively. Using a parabolic approximation for the spectrum $\epsilon_{e/h}(\mathbf{k})\approx\pm(\tilde{\gamma}_{0}^{2}/\gamma_{1})k^{2}+const.$, we find

[TABLE]

Hence, we may interpret this alteration as a correction to the chemical potential and replace for convenience in Sec. IV $f[\epsilon_{e/h}(\mathbf{k}),\tilde{\mu}]\rightarrow f[\epsilon_{e/h}(\mathbf{k}),\tilde{\mu}(\mathbf{k})]$ where we substituted

[TABLE]

Apparently, the displacement of the Fermi contours affects the sharpness of the absorption peaks. To resolve the splitting of the valley-dependent absorption peaks, it is required that

[TABLE]

where the relevant wave vectors are to be evaluated at the Fermi energy, i.e., $|\bm{\mathcal{E}}\cdot\mathbf{k}|\leq\mathcal{E}k_{F}$. It is particularly illuminating to relate the above condition to the transport mean free path $l_{tr}=v_{F}\tau_{tr}$ with the Fermi velocity $v_{F}=2\tilde{\gamma}_{0}^{2}k_{F}/(\gamma_{1}\hbar)$ as it simplifies to

[TABLE]

As at the same time $\Delta_{\text{ax}}$ is proportional $\mathcal{E}$, the condition is independent of the magnitude of the static electric field. Selecting $\mathbf{B}\parallel\bm{\mathcal{E}}$ and inserting the numerical value of $\xi_{\parallel}$ for bilayer graphene, Eq. (26), we may further estimate

[TABLE]

Obviously, due to the smallness of $\xi_{\parallel}$ in bilayer graphene this condition is hardly fulfilled for realistic values of the magnetic field and the mean free path.

Secondly, it should be noted that for low chemical potentials, the Fermi energy may lie very close to the charge-neutrality point where the role of disorder is generally more delicate.Li et al. (2016) In the context of optical conductivity, one should account for the broadening $\sim\hbar/\tau_{tr}$ of the Drude peak at $\omega=0$. Hence, the minimum absorption peaks should occur at $|\omega_{\pm}|\gg 1/\tau_{tr}$ which may be controlled by choosing an appropriate chemical potential and the strength of the axion shift.

Thirdly, to apply the linear response theory, the oscillatory external electric field $\delta\bm{\mathcal{E}}(t)$ with amplitude $\delta\mathcal{E}$ should yield a small perturbation to the kinetic energy. In other words, the changes of the Fermi wave vector $\delta k_{F}$ due to the minimal coupling to the electric field $\delta k_{F}=e\delta A/\hbar$, where $\delta\bm{\mathcal{E}}=-\partial_{t}(\delta\mathbf{A})$ yields $\delta A=-\delta{\mathcal{E}}/\omega$, should be small, i.e., $\delta k_{F}/k_{F}\ll 1$. Therefore, we obtain the condition

[TABLE]

Using again a parabolic approximation and disregarding the drift corrections due to the static field $\bm{\mathcal{E}}$, the Fermi wave vector $k_{F}$ obeys the relation $k_{F}\approx\sqrt{\gamma_{1}|\tilde{\mu}\mp\Delta_{\text{ax}}|}/\tilde{\gamma}_{0}$ in the valley $\pm\mathbf{K}$. Thus, the magnitude of the oscillatory external electric field $\delta\bm{\mathcal{E}}(t)$ should be chosen according to this condition. For larger fields, higher-order terms in the Kubo formula should be taken into account.