Magnetoelectric generation of a Majorana-Fermi surface in Kitaev’s honeycomb model

We study the effects of static magnetic and electric fields on Kitaev’s honeycomb model. Using the electric polarization operator appropriate for Kitaev materials, we derive the effective Hamiltonian for the emergent Majorana fermions to second-order in both the electric and magnetic fields. We find that while individually each perturbation does not qualitatively alter Kitaev spin liquid, the cross-term induces a finite chemical potential at each Dirac node, giving rise to a Majorana-Fermi surface. We argue this gapless phase is stable and exhibits typical metallic phenomenology, such as linear in temperature heat capacity and finite, but non-quantized, thermal Hall response. Finally, we speculate on the potential for realization of this physics in Kitaev materials.


I. INTRODUCTION
Topological states of matter have attracted broad interest due to their fundamental importance in our understanding of many-body systems [1], as well as their potential practical importance in storing and manipulating quantum information [2]. A key role in our understanding has been played by exactly solvable models, such as the toric code [3], where the nature of the ground state and fractionalized excitations is indisputable. However, finding and exploring the physics of topological phases of matter, such as spin liquids, in more realistic models is difficult, with fewer tractable systems to study [4].
The effect of an electric field on the Kitaev spin liquid is much less well-understood. Its response is encoded in the (effective) electric polarization operator [57] appropriate for Ki- taev materials, as was worked out in detail in Refs. [58][59][60]. Using these polarization operators, the dielectric response can be computed [59], providing a natural route to subgap response in the optical conductivity. However, only the dynamic (linear) response has been considered, leaving open the question: how does the Kitaev spin liquid respond to a static electric field?
In this article, we address this question, considering the effect of combined magnetic and electric fields on Kitaev's honeycomb model. Using a generic, symmetry constrained polarization operator [58,60], we apply degenerate perturbation theory to compute the effective Hamiltonian to secondorder in both the magnetic and electric fields. We find that while the pure electric and magnetic contributions do not fundamentally alter the Kitaev spin liquid at this order, the leading Magnetoelectric effect induces a chemical potential at the Dirac touching points and give rise to a Majorana-Fermi surface [61]. This gapless spin liquid phase has no instabilities with respect to arbitrary perturbations and manifests in sig-natures in thermodynamic quantities. We further explore the interplay of this second-order cross-term with third-order contribution from the magnetic field, which stabilizes the gapped chiral spin liquid phase. We find that at low temperature it gives rise to a finite, but non-quantized [62], thermal Hall response, lim T →0 κ xy /T , in the gapless phase.
Magnetoelectric effects [63,64] in ferromagnets have long attracted intense interest due to the potential for electrical control of magnetism (and vice versa) and for a variety of applications in spintronics. Applications in frustrated magnets (and anti-ferromagnets more broadly [65]) are not as well explored [1]. However, several results have established the potential utility of electrical probes, from subgap optical response in gapless spin liquids [59,66,67] to allowing electric control of fractionalized excitations in spin ice materials [68,69]. Here we offer Magnetoelectric route to generating a Majorana-Fermi surface. Since this is due to application of external fields, this has several advantages over more drastic perturbations; avoiding, for example, some of the complications of the quantum chemistry involved with doping or application of pressure.
The appearance of such Majorana-Fermi surfaces in Kitaev-like models has been discussed in several contexts, however each case carries some fundamental difficulty. These include instability towards nodal phases [29,31] (with nodal lines or points), sensitivity to symmetry-allowed perturbations [70] (fine-tuning) or realization through models with very unconventional lattices or interaction terms [71][72][73][74][75] or with site-dependent (staggered) applied magnetic fields [76]. In contrast, our result offers a natural avenue towards a Majorana-Fermi surface in a model that is directly related to realistic models of Kitaev materials. With this in mind, we discuss what kind of electric field strengths would be necessary to observe this physics in an idealized realization of Kitaev's honeycomb model; we find that, while large, the required electric fields are not far outside experimental reach.
The article is structured as follows: in Sec. II we define the Kitaev model and outline the derivation of the symmetry allowed polarization operator. In Sec. III we review the exact solution of the pure Kitaev model and its symmetries to establish our notation and conventions. Sec. IV describes our main results, covering the derivation of the effective Hamiltonian to second-order in the magnetic field alone (Sec. IV A), in the electric and magnetic fields (Sec. IV B) and in the electric field alone (Sec. IV C). This effective Hamiltonian is then solved in Sec. V using the Majorana representation, and its spectrum, including the appearance of the Majorana-Fermi surface, is discussed in Sec. VI. In Sec.VII we explore the dependence on the direction of the electric and magnetic fields, focussing on the competition between the gap-opening term at third-order in the magnetic field and the Majorana-Fermi surface inducing second-order contribution. Finally, in Sec. VIII we discuss relation to some recent works on perturbed Kitaev models, the effects of Majorana interactions, and provide some rough estimates of electric field strengths required to realize this physics, before concluding in Sec. IX with an outlook and some perspective on future directions.

II. ELECTRIC POLARIZATION IN KITAEV MATERIALS
To begin, we define an effective j = 1/2 pseudo-spin Kitaev model on a honeycomb lattice, as might appear in an ideal Kitaev material [5] where i j γ are the (labelled) nearest-neighbour bonds (see Fig. 2). Generically, in the presence of both an electric and magnetic field we expect this Hamiltonian to acquire two new terms, taking the form where S i ≡ σ i /2 are the pseudo-spins, B is the magnetic field, M ≡ gµ B i S i is the magnetization operator, E is the electric field, and P is the electric polarization operator.
To proceed, we need to express the electric polarization operator, P, in terms of the pseudo-spins appropriate for the Mott insulating regime. This can be done perturbatively from the atomic limit, as discussed in Ref. [57]. For Kitaev materials specifically, the structure of these polarization operators was worked out from similar microscopic considerations in Refs. [58][59][60], taking into account the details of the physics of 4d or 5d transition metal oxides. Rather than embark on such a microscopic approach, we will follow the discussion in Ref. [60], and parametrize the polarization operator using a generic form only constrained by lattice symmetries.
First, let us quickly review the derivation of the form of the polarization operators allowed by symmetry. Since the electric polarization operator is time-reversal even, inversion odd and translational invariant, we can immediately see that it must take the form where (by convention) we order each bond i j so that i belongs to sublattice A and j to sublattice B. We have truncated the expansion of this operator to nearest neighbours and only two-spin terms due to the structure of the perturbative expansion [57], with other contributions appearing at higher order in t/U. To go further, we need to know the details of the lattice symmetries. We use the conventional basis for the pseudospins defined such that, in the idealized limit, the octahedral cage of ligands are located along the ±x, ±ŷ, ±ẑ directions. For simplicity, both the electric and magnetic fields are also defined with respect to this basis. Trigonal distortion of the ligand cage generally lowers the local site symmetry of the transition metal ion to D 3d (ignoring any small monoclinic distortions). The remaining symmetries of the crystal are the C 3 symmetry along the direction perpendicular to the honeycomb plane, and the C 2 symmetries along the nearest neighbour bonds. The three-fold symmetry links the three components of P, so we can simply focus on P z , recovering P x and P y by applying C 3 rotations. This leaves nine parameters in p z γ . Under the action of the bond aligned C 2 symmetries one has that P z → −P z for the z bond. To proceed, define the bond directionsû γ aŝ We further define an orthonormal frame for each bond (û γ ,û γ ,ŵ γ ) whereŵ γ =γ andû γ ≡ŵ γ ×û γ . Under its respective bond C 2 symmetry one has thatû γ is invariant whilê u γ andŵ γ change sign. The bond symmetry then implies that u z · p z z = 0 and u z · p z x = −û z · p z y ,û z · p z x = +û z · p z y ,ŵ z · p z x = +ŵ z · p z y .
These four relations then leave us with five parameters, which we denote as m 1 , . . . , m 5 . After some rearrangement we can write the terms in the polarization operator [Eq. (3)] as [60] p µ γ ≡ m 1û µ γûγ +v µ γ m 2ûγ + m 4ŵγ +ŵ µ γ m 3ŵγ + m 5ûγ . (5) The large number of free parameters allowed in the polarization operator echoes the same freedom in the generic exchange Hamiltonian (four parameters) due to the relatively low bond symmetry [11]. Before moving on to the effects of the electric field on the physics of the Kitaev model, let us quickly note that one might consider including contributions to the magnetization operator at next to leading order (three-spin terms), or so-called orbital contributions [77]. However, since these both appear at higher order than the two-spin terms that appear in the polarization operator, we will neglect them here.

III. REVIEW OF THE EXACT SOLUTION
We present here the structure of the exact solution which was originally proposed by Kitaev, as discussed in Ref. [5], to review the general ideas needed to discuss the effects of polarization and to establish our notation and conventions.
Consider the pure (isotropic) Kitaev model defined by the following Hamiltonian where we have transitioned to using the Pauli operators σ i , rather the spin-1/2 spins S i , with J ≡ −K/4 to simplify some of the later algebra. For simplicity we will assume K < 0 so that J > 0, without any loss of generality [5]. This model can be exactly solved, (partly) due to the large number of conserved "flux" operators that commute with both the Hamiltonian and with each other. We can define these flux operators for each honeycomb plaquette as a product of the spin operators going around the plaquette, with the spin component being that of the outward pointing bond type. Explicitly, for a plaquette p where the indexes start from the leftmost site and run clockwise (see Fig. 2). It is easy to verify that when defined this way each flux operator commutes with every other flux operator and with the Hamiltonian, so that for any plaquettes p and p . Since W 2 p = 1, each flux operator has eigenvalues given by w p = ±1 and is a Z 2 degree of freedom. Consequently, along with the Hamiltonian, they can be simultaneously diagonalized and the full Hilbert space can be decomposed into different flux sectors, each corresponding to the choice of the W p eigenvalues {w 1 , w 2 , ..., w n }. If there are N sites, then there are N/2 plaquettes, and so fixing the flux sector halves the dimension of the Hilbert space, leaving 2 N/2 degrees of freedom. The remaining "half" degree of freedom per site immediately suggests that these could be Majorana fermions.
To make this observation manifest, we follow Kitaev [5] and write the Hamiltonian directly in terms of a Majorana representation and c i are each Majorana fermions. To project back into the physical Hilbert space of the spins we must impose the constraint that Practically, this allows for multiple representations of the spins, e.g σ i ≡ D i σ i = −ib i × b i . Since this constraint commutes with the Hamiltonian, the physical eigenstates of the spin model can be obtained from the eigenstates in the extended Hilbert space by projection.
Using these Majorana fermions the Kitaev model in the extended Hilbert space can be written as where we have defined the Z 2 gauge field operators U i j ≡ ib γ i b γ j on each link. Similar to the Z 2 flux operators, these gauge field operators all commute with each other and with the Hamiltonian and thus can be simultaneously diagonalized in the extended Hilbert space. Since U 2 i j = 1 their eigenvalues are simply u i j = ±1. In terms of the u i j gauge-field the eigenvalues of the flux operators can be interpreted as the corresponding Z 2 gauge flux w p = u p 1 p 2 u p 2 p 3 u p 3 p 4 u p 4 p 5 u p 5 p 6 u p 6 p 1 .
The extended Hilbert space then decomposes into spaces where the u i j are fixed and the effective Hamiltonian is a free Majorana problem given by One can show [5,78] that the ground state sector is such that the Z 2 gauge fluxes are equal to one and thus the u i j can be chosen to be uniform with u i j = +1 when going from the A sublattice to the B sublattice (up to gauge redundancy). This describes free Majoranas on a honeycomb lattice with only nearest-neighbour hopping; the spectrum is thus identical to that of graphene [79], with the dispersion having linear touching points at the corners of the Brillouin zone. This can be made precise by defining the Fourier transformed operators on each sublattice where we have defined each site i by a unit cell r and a sublattice α = A or B and we note that c † k,α = c −k,α . The free Majorana Hamiltonian for the ground state sector is then where the sum runs over half the Brillouin zone and we have defined where a 1 ≡ (3x + √ 3ŷ)/2, a 2 ≡ (3x − √ 3ŷ)/2 are the basis vectors of the honeycomb lattice. By diagonalizing this matrix we obtain the spectrum (k) ≡ ±| f (k)|. As in graphene, we is a corner of the Brillouin zone, to obtain One thus has a linear spectrum (K + q) ≈ ±v|q|, with Dirac velocity v ≡ 3J.

A. Projective Symmetries
Before moving on to the effects of perturbations on the Kitaev liquid, we first review how symmetries of the spin model act in this Majorana basis. We focus our attention on the constraints imposed by inversion symmetry (broken by an electric field) and time-reversal symmetry (broken by a magnetic field) within the ground state flux sector. The key property of both the symmetries is that they are implemented projectively [5,80,81], that is the application of the symmetry operation must be followed by a Z 2 gauge transformation.
Consider first time-reversal: note that time-reversal, T , is anti-unitary and maps σ i → −σ i . In the Majorana representation a perfectly valid time-reversal operation is simply b i → b i and c i → c i with the imaginary prefactor giving the change in sign. However, this changes the link variables, as u i j → −u i j due to their imaginary prefactor. Since we would like to work within the fixed gauge sector with uniform u i j = +1 we can undo this via a gauge transformation with a staggered sublattice sign. Our final (effective) time-reversal operator is then where (−1) i is +1 on the A sublattice and −1 on the B sublattice. The treatment of inversion symmetry, I, is essentially the same: instead of being anti-unitary, it interchanges the two sublattices, changing the sign of u i j in the same way. The effective action of inversion is then where I(i) is the site to which i is mapped to under inversion. These symmetry operations constrain the terms that can be generated by the electric and magnetic fields. For example, since B is odd under time-reversal any free Majorana terms like ic i c j that are generated must respect that symmetry. This means, e.g. that any O(B 2 ) terms must connect different sublattices, while any O(B 3 ) terms must only connect the same sublattice. Similarly, all of the pure electric field terms must connect different sublattices, as they are all time-reversal even.
Inversion symmetry is less restrictive; though it does not necessarily preserve the pair of sites in question (unlike timereversal), one can still make some (more limited) statements. For example, since the first and third nearest-neighbour bonds are preserved by inversion they cannot be generated at odd order in E. For second nearest neighbour bonds, inversion only relates the hopping on one bond to the distinct inverted bond. Finally, let us mention the cross term, at O(EB) -the main focus of our work -which is odd under both time-reversal and inversion and thus appears first in the second neighbour bonds, but is distinct from the usual O(B 3 ) contribution.

IV. PERTURBATION THEORY IN ELECTRIC AND MAGNETIC FIELDS
We now review the usual perturbative approach to obtaining effective Hamiltonians for the ground state flux sector when small perturbations are added. To set the stage, we will first recap known results [5] for the effects of a magnetic field at O(B 2 ) in Sec. IV A. We then proceed to derive the O(EB) contributions in Sec. IV B and the O(E 2 ) contributions in Sec. IV C. Relevant aspects of the O(B 3 ) contributions are reviewed [5] in App. A.

A. Magnetic Field
The effect of magnetic field is encapsulated in the piece of the Hamiltonian where we define the (reduced) magnetic field h ≡ gµ B B/2. The form of perturbation theory to be used is motivated by the observation that the action of a single-spin operator always changes the flux sector. To see this consider the effect of the single spin σ z i acting on the system -one small piece of the magnetization M z . Using the commutation relations of the spins, we can see that when acting on state from the ground state flux sector where W µ is the plaquette operator opposite to the µ-bond connected to site i. We thus see that we have added fluxes on a pair of hexagons that are connected to the site where we acted the spin operator. When considering the effect of this perturbation, the virtual states generated will thus not be within the ground state flux sector, but will necessarily mix in the two-or higher-flux sectors. Graphically, the flux configurations generated by the three spin components can be illustrated as where a filled hexagon indicates that w p = −1 on that plaquette and an empty hexagon indicates w p = +1. We thus consider a form of quasi-degenerate perturbation theory, where we derive an effective Hamiltonian within the ground state flux sector. Given the dimension of this Hilbert space, it also naturally admits a description in terms of a single Majorana per site, i.e. the c i fermions. Formulating this perturbation theory strictly requires consideration of the full multi-particle spectra of the virtual flux states -which is a much more challenging task. We will instead follow the approach of Kitaev [5] and assume that most of the weight in the virtual processes comes from the single-particle excitations, an assumption that has been made plausible by more recent numerical studies [43,82]. Practically, this means that we will replace any resolvents in our perturbation theory with a single energy scale -the relevant flux gap.
To see how this is carried out, define the projection operator P 0 that projects into the ground state flux sector. In a magnetic field the effective Hamiltonian would then be (at second order) where H 0 is the Kitaev model and we have used that the resolvent R is approximately given by R = (1−P 0 )/∆, with ∆ being the gap to creating two neighbouring flux pairs. From Ref. [5], one can estimate this to be ∆ ≈ 0.2672|J| = 0.067|K|. Since P 0 MP 0 = 0 due to the change in flux, we see the effective Hamiltonian is then At this order, from the above considerations of flux generation [Eq. (20)], we can further see that The first term describes adding two fluxes by applying a spin σ µ i at one site, and removing them using the same operator yielding an unimportant constant. The second term describes removing the added fluxes by applying the nearest-neighbour σ µ i+µ and yields something non-trivial. We thus obtain The leading effects of the field are thus to renormalize the Kitaev couplings to render them (potentially) anisotropic, depending on the field direction. The presence of these O(h 2 ) contributions also implies a finite magnetic susceptibility [34] at zero temperature. Note that the factor of two arises as the operators adding and removing the flux are different, and thus can be applied in two different orders.

B. Electric and Magnetic Field
We now consider the effects of the electric polarization operator, following the same perturbative scheme that we used for the magnetic field (Sec. IV A). The leading, and most interesting, term will be the combination of the electric and magnetic field at O(EB). The first contribution from the electric field alone appears at O(E 2 ) and is somewhat more complicated, without changing the essential physics of the leading term; we will cover it in detail in Sec. IV C. To make the bookkeeping simpler, we introduce the notation (25) where we have defined ε γ ≡ µ E µ p µ γ /4. Explicitly, in terms of the parameters of the polarization operator [Eq. (5)], we have where α, β, γ are a permutation of x, y, z.
We first need to confirm that key property of the magnetic field perturbation that motivated our quasi-degenerate perturbation theory: that the action of the polarization operator changes the flux sector. To see how this works, consider the contribution of a single z-bond i j z to E · P, which takes the form From the structure of this term we can see that while two spin operators are involved in each (in contrast to the single spin for the magnetic field), they are always different spin components due to the anti-symmetry imposed by inversion symmetry. This means that the fluxes generated by one spin operator are not removed by the other -if one enumerates all the possible combinations, one can see that the two pairs of fluxes must only share at most a single plaquette, and therefore we are left with a pair of fluxes, just like in the magnetic field case. This then immediately implies the first order term is zero, with P 0 PP 0 = 0.
The specific combinations of fluxes can be directly inferred from Eq. (20), but are most clearly illustrated graphically. We delineate two types of flux configurations: those that give rise to a pair of nearest neighbour fluxes (type I) and those that give second-neighbour fluxes (type II). For the z-bond discussed above, four of the operators give type I flux configurations, as illustrated below Note that these type I configurations are only generated by the m 1 , m 2 and m 5 parts of the polarization operator. The remaining two operators, coming fromẑ · (σ i × σ j ), and generated by the m 3 and m 4 polarization operators, give the type II flux configurations Note that we have not explicitly written the signs or prefactors in these expressions, we are simply illustrating the flux content of the generated states. The related patterns for the other types of bonds can be inferred using the three-fold symmetry; for a µ-bond, then the operators corresponding tô ν · (σ i × σ j ) where ν µ give the type I configurations, while the operators fromμ · (σ i × σ j ) give the type II configurations.
To make this more explicit, write the O(EB) correction as where again we have used that the resolvent reduces to R = (1 − P 0 )/∆ for the type I intermediate states of Eq. (27). Since P 0 MP 0 = 0 we can see that this reduces to We can use our knowledge of the intermediate (type I) flux states to simplify this further; we will work out one case explicitly, deriving the rest using symmetry. Focus on the contributions to a set of three sites (i, j, k) that define a second nearest neighbour z-bond type bond, as shown in Fig. 3. There are two processes that involve these three sites and both give a non-trivial contribution to the effective Hamiltonian. Together they give where the overall factor of two accounts for the Hermitian conjugate processes where the magnetic field is applied first. This can be generalized and other bond types can be obtained by cyclically permuting the components of all vectors; one finds the final O(EB) Hamiltonian to be − 2 where 2 i j α(β)γ indicates a second nearest-neighbour bond of type β from i to j, i.e. we get from i to j by traversing an α-bond to an intermediate site, then on to j via a γ-bond (see Fig. 6).

C. Electric Field
We now move to the processes that arise at second order due to the electric field alone. We identify three types of distinct processes: one that generates only two-spin interactions, renormalizing the nearest neighbour couplings and two further processes that generate four-spin interactions that link third and fourth nearest neighbours. These four-spin interactions are closely related to those introduced phenomenologically in Ref. [70], see Sec. VIII for a more detailed discussion.

Two-Spin Process (Type II-Type II)
The first process of interest is generated by the type II flux configurations that did not contribute at O(EB), as shown in Eq. (28). Due to their geometrical arrangement, we see that only the operators that generate these further neighbour flux pairs [Eq. (28)] are the same ones that can remove them. For example, for a i j z bond the two options are σ x i σ y j and σ y i σ x j . Since using the identical operator simply results in a constant, only the cross terms give non-trivial contributions to the effective Hamiltonian; we write where ∆ ≈ 0.2372|J| ≈ 0.0593|K| is the gap for creating two further neighbour pairs [5]. This can be done for each bond type, leading to the total effective Hamiltonian contribution Similar to the case of the O(B 2 ) contributions, this simply renormalizes the bare Kitaev couplings, and (potentially) renders them anisotropic. Note that such processes do not exist for the type I flux configurations: the flux pair generated by each operator (on the same bond) is unique, and so the only way to remove them is by applying the original operator again. As in the type II case, this simply gives an unimportant constant.

Four-Spin Process (Third Neighbour)
We now consider processes that involve the type I flux configurations, but with operators on different bonds. We first consider the process where the two pieces of the polarization operator are separated by a nearest neighbour bond and are non-parallel. Concretely, we can consider the processes illustrated in Fig. 4 that can be associated with a z-type third neighbour bond (see Fig. 6). The first such process contributes (taking into account the reversed, or Hermitian conjugate, process as well) where i, j, k, l are the four sites going clockwise along the top of the hexagon (see Fig. 4). A similar process can be written for the sites running along the bottom of the hexagon, labeled i, r, s, l going counter-clockwise, giving the final contribution Identical contributions exist for each third nearest-neighbour bond. We can then write − 2 where 3 i j αβ(γ) is a third-neighbour bond of type γ (associated with the corresponding nearest-neighbour bond, see Fig. 6).

Four-Spin Process (Fourth Neighbour)
Finally, we consider the process where the two pieces of the polarization operator are separated by a nearest neighbour bond and are parallel. Explicitly, we can consider the processes illustrated in Fig. 5 that can be associated with a zy-type fourth neighbour bond (labeled using the composing nearestneighbour bonds, see Fig. 6). This contribution gives where the path i, j, k, l is shown in Fig. 5. Unlike the previous type of process (Sec. IV C 2), this is the only contribution that involves these two endpoints. We thus can write the full set of contributions as where 4 i j αβ(γ) is an αβ-type fourth neighbour bond.

V. SOLUTION OF EFFECTIVE HAMILTONIAN
With the effective Hamiltonian in the zero-flux sector worked out to second order in both the electric and magnetic fields, we now move on to the solution of this Hamiltonian using the Majorana representation.
The simplest terms are simply those that renormalize the nearest neighbour couplings. It is useful to define the (induced) anisotropic Kitaev exchanges J γ for each bond as (36) with d α being the three (outward) nearest-neighbour bond directions of the honeycomb lattice, starting from the A sublattice.
Next we consider the O(EB) contributions: the three-spin term generated at this order [Eq. (31)] is a close analogue to the three-spin interaction that was obtained by Kitaev [5] when going to third-order in the magnetic field. However, a key difference is in the staggered sublattice sign which is required for the operator to be odd under inversion, and thus appear at O(EB). These additional signs qualitatively change the effects of the interaction on the Majorana spectrum. Instead of opening a gap, as the O(B 3 ) term does, this O(EB) term will instead yield a one-dimensional manifold of zero energy states -a Majorana Fermi surface.
To see this, start by writing the O(EB) terms defined in Sec. IV B [Eq. (31)] using the Majorana fermions as where we have made use of the constraint through the replacement Now since u i j u jk = −1 in the zero-flux sector (due to bond orientations), the O(EB) contribution to the effective Hamiltonian for the Majorana fermions is given by This is simply a second neighbour hopping, anisotropic in space and opposite in sign between the two sublattices. In Fourier space this can be written where we have defined We have used that, in general, the anti-symmetry of the Majorana operators imposes that g(k) = −g(−k), which is indeed satisfied by the above definition. Note that all of the other second order terms connect different sublattices, and so (at this order) this is the only contribution to g(k). This is true generally for the terms generated by the electric-field alone, given they are time-reversal even, they cannot provide any contributions to g(k). Finally, we consider the two types of four-spin interactions that are generated at O(E 2 ). Start first with the third-neighbour type (Sec. IV C 2), looking at the contributions to a z-type bond (as illustrated in Fig. 4 The same manipulations on the second process yield identical results. Generalizing to the full set of these four-spin terms, we therefore have the contribution − 4 where we have reversed the sign, by ordering the bonds so that i ∈ A and j ∈ B. In Fourier space this gives a contribution to f (k) of where the remaining indices are such that α, β γ. A similar procedure can followed for the fourth neighbour bonds (Sec. IV C 3); we simply quote the final result − 2 Again, similarly, in Fourier space this gives a contribution to f (k) that goes as The final result for f (k) ≡ f 1 (k) + f 3 (k) + f 4 (k) can be summarized as where the J γ depend on the fields Eq. (35) and the sums in the final two terms have the same meaning as in Eqs. [43,45]. The free Majorana Hamiltonian for the ground state sector, including the O(B 2 ), O(E 2 ) and O(EB) contributions, is then given by The spectrum is then shifted from the case without the electric field It will be useful to include some aspects of the O(B 3 ) contributions to the effective Hamiltonian into our analysis, so that we can explore the competition with the gap-opening terms. To this end we add the O(B 3 ) second-neighbour hopping [5] − For details of the derivation of this term, see Appendix A or Ref. [5]. Note that we have not included the additional four-Majorana interaction term that is also generated at O(B 3 ) [5]. In Fourier space this yields where we have defined the new dispersion function where α(β)γ is defined as a sum over β with αβγ = +1. With this term included, the spectrum of the Majorana fermions then takes the form With E = 0, and thus g(k) = 0, this is the usual gapped spectrum expected for the Kitaev model in a small magnetic field [5]. Expanding near the Dirac points, one has The energy gap is then given by 2|h(K)| for E = 0; it will be useful to define a mass of the Majorana fermions as m ≡ |h(K)|.

VI. MAJORANA FERMI SURFACE
How does the electric field perturbation affect the spectrum of the Majorana fermions, in particular the Dirac point? First, we should note that the O(B 2 ) and O(E 2 ) corrections preserve the Dirac touching, but shift them from ±K. Generically, this shift is also accompanied by the introduction of anisotropy and renormalization to the Dirac velocity as well. Consider the Dirac point near K, schematically one has (53) where δK is the shift of the Dirac point, R θ is the rotation of the principal axes of the (anisotropic) Dirac cone, and v 1 , v 2 are the two-independent Dirac velocities. The shift δK and the rotation angle θ are both second-order in the fields, while v n = v + O(E 2 ), O(B 2 ). In principle, the leading corrections can be worked out explicitly from an expansion of Eq. (46), though we will leave this implicit for the sake of brevity.
For the g(k) part things are somewhat simpler; since δK ∼ O(B 2 ), O(E 2 ) and and g(k) is already O(EB), we only need to evaluate it at K, with any corrections from δK being O(EB 3 ) or O(E 3 B). Evaluating this in terms of the physical electric and magnetic fields, E and B, one thus has g(K) being wheren = (x +ŷ +ẑ)/ √ 3 is the direction perpendicular to the honeycomb plane. We thus see that the cross-term introduces a finite chemical potential at the Dirac points -we define µ ≡ −g(K). A similar argument applies for h(k), since it is O(B 3 ), any shifts due to the second-order terms only have higher order effects, and thus we can take |h(K)| ∼ 72 √ 3h x h y h z /∆ 2 .
From the structure of E and B dependence we see that this chemical potential vanishes in several high symmetry configurations, including parallel electric and magnetic fields, as well as with either field being perpendicular to the honeycomb plane. This term is maximal when the electric and magnetic field are crossed, i.e. E· B = 0, and both are in the honeycomb plane. We also see that the dependence field is distinct from other field induced terms, i.e. this term can remain finite when the gap-inducing O(B 3 ) term vanishes. We also see that not all of the contributions to the polarization operator are effective in generating this chemical potential -both m 3 and m 4 do not contribute as they only generate type II flux configurations.
The spectrum near the Dirac point is then given by  Fermi wave-vector The Fermi velocity at this Majorana Fermi surface is inherited from the Dirac cone, with v F ∼ v, due to the linearity of the dispersion. Note that we have ignored the O(B 3 ) mass terms here, as they are parametrically smaller than the second-order terms that generate the chemical potential. Just as for a more conventional Fermi surface, the additional low-energy excitations present qualitatively change the thermodynamic properties of this state at sufficiently low temperatures. For example, for temperatures well below the induced O(EB) terms, the specific heat of the perturbed Kitaev model is O(T ), rather than the O(T 2 ) in the original model, leading to a finite linear specific heat coefficient, C/T ∼ γ 0 as T → 0.
In the next section, we will confirm this simple picture using a more complete calculation of the spectrum, as well as explore what happens when some of these terms vanish due to geometrical effects.

VII. RESULTS
For concreteness, we consider a configuration of the electric and magnetic fields that allows us to use the different geometrical dependences of the contributions to the effective Hamiltonian to isolate different aspects of the physics. We consider a configuration that can tune smoothly between regimes where the O(B 3 ) contributions vanish, and those where the O(EB) contributions vanish. For each case, we present the Majorana spectrum, and a few key physical observables, calculated (numerically) using the complete spectrum [Eq. (51)]. For practical purposes we need to make some choices in our free parameters. First, we fix the polarization constants [Eq. (5)] to all be equal m 1 = m 2 = m 3 = m 4 = m 5 ≡ m 0 . This is a completely arbitrary choice, and is made simply to control the complexity of presentation -we expect that the qualitative features of our results will not depend strongly on different, but still generic choices of the m n . Next, we renormalize the energy scales and the electric and magnetic fields to correspond closely to the reduced variables (h and ε) that appear naturally in Sec. IV. To this end, we set J = 1, and take gµ B /2 ≡ 1, absorbing the g-factor and Bohr magneton into the definition of B. For the electric field, we take m 0 ≡ 4, absorbing their units into E, and choosing the pre-factor to compensate for the corresponding factor in the definition of ε.
More explicitly, the electric and magnetic fields are taken to be E ≡ E cos ψX + sin ψn , (57a) The orientation of these fields relative to the honeycomb plane is shown in Fig. 8. These choices yield a chemical potential which vanishes for ψ = π/2 and is maximal for ψ = 0. Numerically, we expect then the radius of the Majorana-Fermi surface to be q F ≈ |µ|/v ≈ 31.76 · EB cos 2 ψ The effective mass that splits the two Majorana bands then takes the form m = 24 which vanishes for ψ = 0 but is maximal for ψ = π/2. Note that this function also vanishes for ψ 1 ≈ 0.282047π, in between these two limits.
We characterize the behaviour of the Majorana spectrum as a function of the angle ψ by looking at two quantities: the spectral gap, which we define as the energy of the lowest lying excitation, ( or Dirac point) and the band gap, which we define as minimum energy between the two Majorana bands (which is induced by the O(B 3 ) from the magnetic field). Explicitly, we define the band gap as min k | + (k) − − (k)| and the spectral gap as min k,± | ± (k)|. These quantities are shown in Fig. 9 as a function of the angle ψ, with the Majorana-Fermi surface and spectrum near the band gap shown in Fig. 10 for a handful of representative angles. Starting from ψ = 0 where the chemical potential is maximal and the O(B 3 ) term vanishes, both the spectral and band gaps are zero. For small ψ, the band gap becomes finite, while the spectral gap remains zero (labelled as region I). The band gap reaches a maximum as a function of ψ before going to zero at the special value ψ 1 where the O(B 3 ) term vanishes (the spectral gap is zero throughout). Increasing ψ 1 past this point, the Majorana-Fermi surface shrinks (region II) and then disappears near ψ 2 ≈ 0.3487π, and the spectral gap then becomes finite (region III). A gapped spectrum, similar to what is found with only the magnetic field, is recovered as ψ approaches π/2. These qualitative features can be diagnosed by looking at the behaviour of the thermal Hall conductivity of the Majorana-Fermions as T → 0 (Fig. 9), computed via [62] lim T →0 where F (n) x is the Berry curvature of the n th band, V is the volume of the system and A (n) µ (k) = −i n(k)|∂ µ |n(k) is the associated Berry potential. This can be naturally expressed in units of κ 0 = πk 2 B /(6 ) and the total (occupied) Berry curvatures for each band [83]. When the spectrum is gapped, the Ω (n) are the Chern numbers of each band and thus the thermal Hall response is quantized in units of κ 0 .
In region I, with ψ 0, the finite O(B 3 ) contribution to the spectrum induces finite Berry curvature at the bottom of the two Majorana bands, near the location of the band gap. This Berry curvature gives a non-zero contribution to the thermal Hall conductivity, following Eq. (60). However, due the presence of the Majorana-Fermi surface the effectively "occupied" states do not include the full Brillouin zone. Thus in region I κ xy /T is not integer or rational valued when expressed in units of κ 0 , varying continuously as a function of ψ. At the special angle ψ 1 the O(B 3 ) term vanishes, changing sign; the thermal Hall conductivity follows suit, remaining finite and unquantized in region II, but with opposite sign relative to region I. In region III, for ψ > ψ 2 , there is no Majorana-Fermi surface, and the band gap is open and one thus recovers the half-quantized thermal Hall conductivity expected from Ref. [5].

VIII. DISCUSSION
We first discuss the relevance of these results to the body of related theoretical works on the Kitaev model and address the effects of interactions on the Majorana-Fermi surface (Sec. VIII A). We then discuss some estimates for the electric fields that may be required to observe this physics, as well as some experimental hurdles in reaching these field strengths (Sec. VIII B).

A. Theoretical
Let us begin with the relationship of the gapless states in this work to some recent studies in the literature. A set of four-spin interactions, similar to those derived in Sec. IV C at O(E 2 ), were studied in Refs. [70,84]. These are not identical to those studied here; generically the four-spin interactions generated at O(E 2 ) are anisotropic and there is a fixed relationship between the coefficients of the two different types of operators; however, given their similarity, some discussion is in order. In Ref. [70] it is shown that for large couplings these operators lead to a change in the ground state flux sector of the model, stabilizing a rich variety of "vison crystals" [70]. In our analysis, we have explicitly restricted ourselves to the zero-flux sector, assuming that perturbations from the electric and magnetic fields are sufficiently small to leave the flux sector unaffected. The results of Ref. [70] can give a rough idea of the range of validity of this approximation; using their notation, the coefficients of the four-spin terms should satisfy K 3 /K 1 0.1 and K 3 /K 1 0.4. In our notation, ignoring the (complicated) direction dependence of E, this implies the loose criterion m 2 0 E 2 /∆ 0.1 (where E is the electric field strength) to preserve the zero-flux sector. Alternatively, one may take the view that applying a strong electric field may be a route to stabilizing the "vison crystal" ground states described in Ref. [70]. However, given this likely would push our perturbation theory to the limit of its regime of validity altogether, one must be cautious.
Some of these concerns could be addressed by more de-tailed numerical studies. For example, a numerical search for the true ground state flux sector [70] through quantum Monte Carlo studies at finite temperature [17,18]. The degenerate perturbation theory itself could be validated using numerical techniques (such as DMRG or exact diagonalization) on the original model [Eq. (2)].
We also note the work of Ref. [76], which discusses the presence of Majorana-Fermi surfaces in a Kitaev model in the presence of site-dependent magnetic fields and off-diagonal exchanges Γ and Γ [11]. In their proposal, these fields are generated by zigzag magnetic ordering in adjacent honeycomb layers, providing a source of both time-reversal and inversion symmetry breaking. While magnetic ordering in some layers, with a Kitaev spin liquid in other layers may not seem particularly natural, the simultaneous breaking of these symmetries plays a similar role to the presence of E and B in this work.

Stability
An important question that must be addressed for the Majorana-Fermi surface we find in this work is stability to interactions. Similar states with gapless surfaces in threedimensional Kitaev models [85][86][87][88][89][90][91][92], as "Bogoliubov-Fermi" surfaces. For these systems, the absence of this kind of (effective) nesting symmetry precluded such an instability.
Generically, in the case of interest for the Kitaev model in electric and magnetic fields, all of the rotational symmetries of the system are broken, along with time-reversal symmetry. Thus, similar to the case of time-reversal breaking noncentrosymmetric superconductors, the system should be stable to the inclusion of weak short-range interactions. This can be seen from the symmetries of the problem: with finite, generic E and B the only remaining symmetries are discrete translations. Without any nesting vectors, interactions can only link a finite number of patches of the Majorana-Fermi surface and thus the enhancement that typically leads to instability is absent [93,94]. Analogues of superconducting instabilities are also avoided [95], given there is no U(1) symmetry of the fermions left to spontaneously break [91].
However, this argument has a few subtleties, given that the perturbations that generate the Majorana-Fermi surfaces also generate the interaction terms, just at higher order in perturbation theory. To see why this could complicate the analysis, note that in Sec. VII we showed that at second order in E and B the Majorana-Fermi surfaces are ellipses, which enjoy a kind of effective inversion symmetry through their centers. At higher-order, contributions would render these surfaces non-elliptical (see Fig. 7), but would also include Majorana-Majorana interactions on equal footing. We leave the (potential) competition between these two effects, and thus the ultimate fate of these Majorana-Fermi surfaces to future studies.
We note that at the special angle ψ 2 where the spectral gap closes, the Majorana spectrum is that of a quadratic band touching, with a spectrum ∝ |k| 2 . While the thermodynamic signature, C ∝ T , this is the same as the Majorana-Surface, (due to its finite density of states) the effects of interactions FIG. 11. Illustration of z-type coverings. A product of W p over either the shaded plaquettes (or the unshaded plaquettes) yields the operator Σ z used in Eq. (61). The Σ x and Σ y operators are defined similarly, with the plaquettes in the product sharing x-type or y-type bonds. may be different. A similar quadratic band touching point appears in AB stacked bilayer graphene [96]. This has attracted attention due to its instability to weak short-range interactions, which leads to the emergence of a variety of new phases. We leave the discussion of the interaction effects on this quadratic Majorana spectrum, and any potential instabilities, to future work.

Correlations
Finally, we note that the gapless nature of the spectrum in our model should lead to algebraic correlations in the spinspin correlation functions. However, these may be "hidden" [97][98][99]; for example, as in the ideal Kitaev model the spin-spin correlations are ultra-short range in spite of the gapless Dirac spectrum of the Majorana fermions. Although the result of Ref. [97] is not directly applicable in our effective Hamiltonian, it can be applied directly to the original model, Eq. (2), including the electric polarization operator. Explicitly, using the criterion from Ref. [97] if we consider the electric and magnetic field perturbation as V ≡ −E · P − B · M from Eq. (2), then where Σ µ is defined to be the product of flux variables over strings of plaquettes [97], as shown in Fig. 11. We thus expect that gaplessness of the Majorana-Fermi surface to manifest in power-law spin-spin correlations. Alternatively, this can be seen by computing the associated canonical transformation of the spin operators when carrying out the degenerate perturbation theory of Sec. IV, as discussed in Ref. [99] in a related context.

B. Experimental
We now address the potential applicability of these results to the growing family of "Kitaev materials" [6,7,13] that are believed to have dominant Kitaev exchange. Ideally, one would like to apply electric and magnetic fields of magnitude that are realistically attainable and generate a Majorana-Fermi surface that could be observable in thermodynamic, transport or spectroscopic probes.
To do so, we must determine a reasonable range of estimates for the coefficients m i that appear in the (effective) electric polarization operator, as well as take a typical value for the Kitaev exchange K. The latter is straight forward: in transition metal Kitaev materials one typically expects K ∼ 5 meV [7,100]. We thus expect J = K/4 ∼ 1.25 meV and thus the relevant flux gaps are ∆ ∼ ∆ ∼ 0.35 meV. Estimating the m i is more complex, with several distinct microscopic mechanisms potentially contributing at the same order [59,60]. In Bolens et al. [60] these are estimated as the sum of two contributions: m n ≡ aA n + B n where a is the nearest-neighbour distance (we will take it to be as in α-RuCl 3 , a ∼ 3.5Å [101]). The first contribution depends on the detailed choice of atomic and hopping parameters, but is estimated to be as large as where e is the electron charge. The contribution B n is not estimated in Ref. [60], so we will simply neglect it in this discussion. We thus have the estimate of m 0 ∼ 3.5 · 10 −2 eÅ = 3.5 · 10 −9 meV/(V/m). This estimate is very rough and errors as large as an order of magnitude would not be surprising, given the uncertainty in the microscopic physics.
Given this value for m 0 , we can now estimate the electric field strengths required to give a Majorana-Fermi surface of a given size. This size is set by the chemical potential which, for crossed fields (see Sec. VII) is µ ∼ 3 √ 2gm 0 µ B EB/(2∆) [restoring some constants in Eq. (54)]. Taking g ∼ 2 and B ∼ 10 T we arrive at µ ∼ 2.4 · 10 −8 meV/(V/m) E.
Since we are estimating J ∼ 1.25 meV the Dirac velocity is v ∼ 3.75 meV and so the radius of the Fermi surface is given by [Eq. (56)] To get a q F that is large, say q F a ∼ 0.1 we would thus need that E ∼ 10 7 V/m ≡ E 0 . This is a very large field [102]; though we note that static electric fields of of magnitude ∼ 10 6 − 10 7 V/m have been used in studies of magnetic systems [103][104][105][106][107]. This naïve analysis must be supplemented with a number of caveats. First, it must be self-consistent: we have that at these large fields m 0 E 0 ∼ 0.1 meV which is not significantly smaller than the flux gap ∆ ∼ 0.35 meV, so our perturbation theory may be approaching the edges of its validity. Further, we have the requirement that we must remain in the zero-flux sector [70]. For this large electric field m 2 0 E 2 0 /(J∆) ∼ 0.02; reasonably far from where one might expect a transition to a new flux sector (see discussion in Sec. VIII). [108] Similar considerations must be applied to the magnetic field as well, concerning the validity of its perturbation theory and its effect on the phase boundary of the zero-flux phase. One route to lowering the required electric field would be via larger magnetic fields along a direction where the O(B 3 ) terms vanish. Alternatively, one could search for Kitaev materials where K or ∆ are smaller; rare-earth Kitaev materials [109,110] may be promising alternatives, if the size of electric field coupling is not dramatically changed. Kitaev materials based on metalorganic frameworks [111], which have longer bond-lengths, may also change the balance of the exchange and electric field interactions.
More practically, at large static fields such as these one must also be aware of effects on the physical system that are not included with our minimal model [Eq. (2)]. This includes potential structural distortions of the lattice along the field, the introduction of charge carriers through electrostatic doping as well as the modification of the microscopic exchange processes due to the field. More dramatically, the material itself could experience dielectric breakdown; for Mott insulators breakdown electric fields of order ∼ 10 8 −10 9 V/m are not atypical, given the expected relationship to the on-site repulsion [112]. One potential route to avoid some of these issues is to consider large electric field pulses, perhaps through the application of laser light, that can reach these field strengths over short time windows. Realizing this physics only within a short time window (long with respect to the spin dynamics) would however limit the experimental probes available to characterize the system.
An alternative route to large static electric fields could be through the engineering of heterostructures of Kitaev materials. Recently, heterostructures of graphene (single and bilayer) and (bulk and monolayer) α-RuCl 3 have been reported in the literature [113][114][115][116]. In each case, the associated potential difference (large effective electric field) results in significant doping of both the graphene and the α-RuCl 3 . While the effect of these fields appears to be too strong for our purposes, the possibility of van der Waals heterostructures of Kitaev materials as a platform to explore the electric field physics discussed in this work remains an intriguing possibility. We note that several recent theoretical works have explored the (related) effect of tunnelling in such heterostructures [117][118][119][120] through the effect of localized electric fields.
While each of these complications to achieving the large electric fields needed here are important, these need not be insurmountable. From the properties of the Majorana Hamiltonian these surfaces can appear whenever we break timereversal and inversion symmetry simultaneously; combined electric and magnetic fields are simply the minimal perturbation that does so. Thus, while we have derived our results within a well-defined framework for the pure Kitaev model with a specific (though generic) electric polarization operator, we expect that generation of a Majorana-Fermi surface not to be strictly dependent on these implementation details. For example, physics similar to that introduced here may be introduced by disorder that breaks inversion symmetry. Explicitly, something akin to charged defects could introduce in-plane electric fields and, in regions where these fields are sufficiently uniform, our results would apply. Introducing a magnetic field would then generate a Majorana chemical potential in concert with these defect fields. We leave the exploration of these possibilities for future work.
Finally, let us make a connection to more conventional Magnetoelectric effects that have been extensively discussed in multiferroic materials. A first step to gauge the importance of these kinds of magnetoelectric interactions, before embarking on a detailed search for a Majorana-Fermi surface, will be to look at the bulk magnetization as a function of electric field strength. As in conventional Magnetoelectric materials, due to the presence of the O(EB) terms in the effective Hamiltonian, one expects that at leading order (fixing directions and coupling coefficients) where χ is the susceptibility and α a Magnetoelectric coupling. A similar term, with E and B switched, would appear in the electric polarization. This effect could be measurable, even if making a sufficiently large enough Majorana-Fermi surface is not feasible, and could provide some guidance on the size of the unknown couplings, m n , that appear with the polarization.

IX. CONCLUSION
In summary, we have analyzed the effects of combined static electric and magnetic fields on Kitaev's honeycomb model. Starting from the symmetry-allowed effective polarization operator [59,60], we used degenerate perturbation theory to derive an effective Hamiltonian to second-order in both the magnetic and electric fields. This effective Hamiltonian is solvable and describes a set of (free) Majorana fermions with a Majorana-Fermi surface over a wide range of parameters, including the neighbourhood of the Kitaev spin liquid. We explored the effects of the O(B 3 ) gap-opening perturbation, showing how it competes with the generation of the Majorana-Fermi surface depending on the relative orientations ofÊ and B. The spectrum of this spin liquid with Majorana-Fermi surface was characterized in the weak field limit, where we derived the low-energy form of the dispersion, and determined the thermodynamic properties and thermal Hall coefficient. Finally, we discussed some related work, speculated on the effects of interactions and provided rough estimates of the magnitude of electric and magnetic fields that are required to render the size of the Majorana-Fermi surface significant.
We hope that our results further motivate the use of electric fields as potentially useful perturbations in Kitaev materials, such as α-RuCl 3 . Not only in linear response (as in Refs. [59,60]), but as a tuning parameter that could shed light on whether a Kitaev spin liquid is present and potentially generate a new spin liquid with a Majorana-Fermi surface. Many questions remain, such as addressing the dynamics of this gapless liquid, the effect of temperature and the role of dilution and bond-disorder on this state. The answers to these questions have proven a rich source of insight in Kitaev's model. We hope that future work on liquids with Majorana-Fermi surfaces, as in the models presented in this work, can be similarly fruitful.
where, as in the main text, we have used that the resolvent reduces to R = (1 − P 0 )/∆. From the flux patterns generated by each piece (see Fig. 12), we see that to obtain a combination of operators that leads us back to the zero state flux sector, one needs µ, ν, λ to be a permutation of x, y, z. All of these permutations yield the same operator, and thus can be accounted for with an overall factor of 3! = 6. Noting that P 0 MP 0 = 0 we can see that the only non-zero terms are − 6h x h y h z ∆ 2 i, j,k P 0 σ x i σ y j σ z k P 0 .
There are two choices for {i, j, k} which will give us non-zero contributions. One is given by which just gives an unimportant constant. The second nonzero contribution is generated by a configuration such as the one shown in Fig. 12.
To make this more explicit, write the O(B 3 ) correction (illustrated in Fig. 12) as − 6h x h y h z ∆ 2 i∈A P 0 σ x i+x σ y i+y σ z i+z P 0 + i∈B P 0 σ x i−x σ y i−y σ z i−z P 0 .