Orbital magnetoelectric coupling at finite electric field

# Orbital magnetoelectric coupling at finite electric field

Andrei Malashevich Department of Physics, University of California, Berkeley, California 94720, USA    David Vanderbilt Department of Physics & Astronomy, Rutgers, the State University of New Jersey, Piscataway, New Jersey 08854, USA    Ivo Souza Centro de Física de Materiales and DIPC, Universidad del País Vasco, 20018 San Sebastián, Spain Ikerbasque, Basque Foundation for Science, E-48011 Bilbao, Spain
July 15, 2019
###### Abstract

We extend the band theory of linear orbital magnetoelectric coupling to treat crystals under finite electric fields. Previous work established that the orbital magnetoelectric response of a generic insulator at zero field comprises three contributions that were denoted as local circulation, itinerant circulation, and Chern-Simons. We find that the expression for each of them is modified by the presence of a dc electric field. Remarkably, the sum of the three correction terms vanishes, so that the total coupling is still given by the same formula as at zero field. This conclusion is confirmed by numerical tests on a tight-binding model, for which we calculate the field-induced change in the linear magnetoelectric coefficient.

###### pacs:
75.85.+t,03.65.Vf,71.20.Ps

Magnetoelectrics are magnetic insulators whose dielectric polarization changes linearly under a small applied magnetic field and, conversely, whose magnetization changes linearly with a small applied electric field .odell70 (); fiebig05 () This linear magnetoelectric (ME) coupling is described by the response tensorexplan-units ()

 αij=∂Mj∂Ei=∂Pi∂Bj, (1)

which is odd under both spatial inversion () and time-reversal () symmetries. Thus ME materials must be acentric and display magnetic order.

In crystals where only one of the two symmetries, or , is present, it may still be possible to induce a linear ME effect by applying an external field which breaks that symmetry. So, for example, a centrosymmetric insulating antiferromagnet placed in a (strong) electric field loses its inversion center. Likewise, a nonmagnetic ferroelectric crystal loses time-reversal symmetry when subject to a magnetic field. In both cases the symmetry is sufficiently lowered that the tensor becomes nonzero.

It is useful to view these field-induced effects as higher-order ME responses of the unbiased crystal. ascher-philmag68 () Two quadratic ME effects can be defined in this way. Going to next order in magnetic field yields the tensor

 βijk=∂αij∂Bk=∂2Pi∂Bj∂Bk, (2)

which is odd under and even under . Going instead to next order in the electric field gives

 γijk=∂αji∂Ek=∂2Mi∂Ej∂Ek, (3)

which is even under and odd under . Reference ascher-philmag68, lists the form of these tensors for all the crystal classes. While most investigations of ME couplings in solids have focused on the linear response for a reference state of the crystal at zero electric and magnetic fields, the quadratic responses and have also been measured in materials where vanishes by symmetry. In particular the electric-field-induced effect, which constitutes the primary focus of this work, was first measured by O’Dell in yttrium iron garnet.odell-philmag67 ()

The ME response can be divided into four contributions, depending on whether the response is frozen-ion (purely electronic) or lattice-mediated, and whether it is spin or orbital in character. We will refer to the frozen-ion part of the orbital response as the orbital magnetoelectric polarizability (OMP).essin09 (); malash2010 () While the OMP is typically a small contribution to the ME response in conventional magnetoelectrics, it was recently realized that, under certain conditions of surface preparation, -odd topological insulatorshasan-rmp10 () should display a large, quantized OMP response.qi08 (); essin09 () This is a remarkable prediction, especially considering that in this class of materials symmetry is preserved in the bulk (it must, however, be broken on the surface). This topological magnetoelectric effect has triggered a great deal of interest in orbital magnetoelectric couplings in solids.

The microscopic theory needed to calculate the OMP at zero electric and magnetic fields from first principles was worked out in Refs. malash2010, ; essin2010, . In addition to the so-called Chern-Simons term responsible for the topological ME effect,qi08 (); essin09 (); coh () it was found that two more (Kubo) terms contribute to the OMP in conventional magnetoelectrics in which and symmetries are broken spontaneously in the bulk.

In this work we generalize the band theory of OMP of periodic insulatorsmalash2010 (); essin2010 () to finite electric fields. That is, we evaluate the coefficient at nonzero ,

 αij(E)=∂Mj∂Ei∣∣∣B=0. (4)

(Henceforth, the condition will be implied throughout. It is also understood that from now on denotes the OMP part of the entire ME response.) A principal result of our work is the conclusion that the zero-field expression for the total OMP remains valid at finite electric field, while the above-mentioned Chern-Simons and Kubo terms separately acquire field-induced contributions. We confirm our formal results by numerical tests on a tight-binding model.

Our derivation of a formula for proceeds along the lines of Ref. malash2010, . We start from the expression given therein for the orbital magnetization of a generic band insulator under a finite electrical bias. It comprises three terms,

 Mj(E)=MLCj(E)+MICj(E)+MCSj(E), (5)

where

 MLCj=−η2ϵjpq∫d3kIm⟨˜∂punk|H0k|˜∂qunk⟩, (6)
 MICj=−η2ϵjpq∫d3kIm{⟨unk|H0k|umk⟩⟨˜∂pumk|˜∂qunk⟩}, (7)

and

 MCSj=eη2Ej∫d3kϵpqrtr[Ap∂qAr−2i3ApAqAr]. (8)

The common prefactor in these formulas is ( is the magnitude of the electron charge), and a sum is implied over repeated Cartesian () and valence-band () indices. The cell-periodic part of the field-polarized Bloch statesouza02 () is denoted by , is the partial derivative with respect to the th component of the wavevector , and the tilde indicates a covariant derivative , where (sum implied over ). The Hamiltonian is defined as

 H0k=e−ik⋅rH0eik⋅r, (9)

where is the zero-field part of the crystal Hamiltonian. In Eq. (8) the symbol denotes the Berry connection matrix

 Amnkp=i⟨umk|∂punk⟩, (10)

and the trace is over the valence bands.

Equations (6) and (7) describe respectively the local and itinerant circulation contributions to the magnetiztion,malash2010 () while Eq. (8) is the Chern-Simons term. At variance with the other two terms, whose dependence on the electric field is only implicit, displays an explicit linear dependence on . It is therefore expedient to introduce a new quantity via the relation

 MCSj(E)≡EjMCS1(E), (11)

where the subscript ‘1’ serves as a reminder that enters the expression for multiplied by to the first power.

All three magnetization terms, , , and , are invariant under gauge transformations within the valence-band manifold, although in the case of this invariance is only modulo a quantum of indeterminacy.qi08 () In the limit that goes to zero, vanishes and Eq. (5) reduces to the expression for the spontaneous orbital magnetization.ceresoli06 ()

As already mentioned, all terms in Eq. (5) can contribute to the linear ME coupling, Eq. (4), so that

 αij(E)=αLCij(E)+αICij(E)+αCSij(E). (12)

The derivation of the expressions for these objects is straightforward though somewhat lengthy. It essentially repeats the steps in Appendix B of Ref. malash2010, , where the derivation was carried out for the LC and IC (“Kubo”) terms under the assumption that (the CS term is trivial at ). At one may show that each of the terms in Eq. (12) consists of a “zero-field” part plus a “field-correction” part having an explicit linear dependence on ,

 αij(E)=α0,ij(E)+Ejα1,i(E). (13)

The field-correction terms for the LC and IC contributions can be traced back to Eqs. (B.7) and (B.8) in Ref. malash2010, , which at acquire extra terms. As for the Chern-Simons contribution, differentiating Eq. (11) with respect to yields and .

Thus, we arrive at the results

 αLC0,ij(E)=ηϵjpqIm∫d3k(⟨˜∂punk|(∂qH0k)|˜Diunk⟩−12⟨˜∂punk|(DiH0k)|˜∂qunk⟩), (14)
 αIC0,ij(E)=ηϵjpqIm∫d3k(⟨˜∂punk|˜Diumk⟩⟨umk|(∂qH0k)|unk⟩−12⟨˜∂punk|˜∂qumk⟩⟨umk|(DiH0k)|unk⟩), (15)
 αCS0,ij(E)=δijηe2∫d3kϵpqrtr[Ap∂qAr−2i3ApAqAr], (16)
 αLC1,i(E)=ηe∫d3kϵpqrRe[⟨˜Diunk|˜∂pumk⟩⟨˜∂qumk|˜∂runk⟩], (17)

and

 αLC1,i(E)=αIC1,i(E)=−12αCS1,i. (18)

In the above expressions, is the partial derivative with respect to the th component of the electric field. The terms containing in Eqs. (14) and (15) are screening corrections which are present in self-consistent calculations.

Equations (14)–(16) for the zero-field terms are essentially rewritten from Ref. malash2010, . It should be emphasized, however, that in the present context these expressions depend on the electric field implicitly via the wave functions. The explicit field dependence is given by the field-correction terms, Eqs. (17) and (18). Remarkably, these terms are not independent and add up to zero when inserted into Eq. (12). We conclude, therefore, that the expression for the total OMP derived in Refs. malash2010, ; essin2010, assuming remains valid for . This constitutes one of our principal results. The explicit expression given in Eq. (17) for the field-correction terms is the other main result of this work. It is useful if one is interested in the field dependence of the separate gauge-invariant contributions to the OMP. Because it contains three derivatives and one field derivative, this quantity is even under and odd under , just like the coefficient defined by Eq. (3). This is reasonable since, as one can see from Eq. (13), gives a contribution to and should therefore have the same symmetry properties.

As a check of our analytic derivation, we have implemented the formula for in a tight-binding model, and used it to calculate the nonlinear ME coefficient at . Since the tensor vanishes in -invariant systems, we need a model where is spontaneously broken, and we chose that of Ref. malash2010, . This is a spinless model with eight sites per primitive cell arranged on a cube, where symmetry is broken by complex nearest-neighbor hoppings, and we have used the same on-site energies and nearest-neighbor hoppings tabulated in that work. (This choice of parameters also breaks , so that the linear ME tensor is nonzero already at , but this is not essential for our present purposes.) As in Ref. malash2010, the two lowest bands were treated as occupied, and the phase of one of the complex hoppings was chosen as a control parameter for plotting purposes.

The technical details of the tight-binding implementation of Eqs. (6)–(8) and (14)–(17) can be found in Ref. malash2010, . The only significant difference with respect to that work is that the field derivative of the cell-periodic Bloch states must be evaluated at finite . Under these circumstances the usual “sum-over-states” formula ceresoli06 () cannot be employed, and one must instead minimize a suitably defined functional. wang-prb07 ()

We shall calculate the component of from the first equality in Eq. (3). Combining with Eq. (12) we find

 γ=γLC+γIC+γCS. (19)

The CS term is the simplest to evaluate, as the derivative of Eq. (13) with respect to can be taken analytically. The zero-field and field-correction terms therein both contribute an amount to . Thus,

 γCSzzz(0)=2αCS1,z(0)=−4αLC1,z(0), (20)

where the second equality follows from Eq. (18). The quantity on the right-hand side can be evaluated directly from Eq. (17). For the LC and IC terms we calculate the derivative of the zero-field terms in Eq. (13) using finite differences and obtain

 γLC/ICzzz(0)≃αLC/IC0,zz(Ez)−αLC/IC0,zz(−Ez)2Ez+αLC1,z(0). (21)

In practice we evaluate the first term from Eqs. (14) and (15), using small positive and negative fields along of magnitude .

The results of the above calculations were compared with a finite-difference determination of the second field derivative of ,

 γzzz(0)=∂2Mz∂E2z∣∣ ∣∣E=0≃Mz(Ez)−2Mz(0)+Mz(−Ez)E2z, (22)

using the -space expressions from Ref. malash2010, for the LC, IC and CS terms in Eq. (5). The results obtained in this manner can be taken as a reference, since the -space expression for has been carefully tested by comparing with real-space calculations on bounded samples cut from the bulk crystal.malash2010 ()

The agreement between the two sets of calculations can be seen in Fig. 1, where the LC, IC, and CS contributions to are plotted separately as functions of . In this calculation is about an order of magnitude smaller than . From Eqs. (20) and (21) it then follows that the field-correction terms contribute little, especially in the case of . Further numerical tests focusing on those small terms are therefore desirable.

In order to isolate the field-correction terms in and , we subtract the zero-field terms from the total:

 ⎡⎢⎣∂2MLC/ICz∂E2z−∂αLC/IC0,zz∂Ez⎤⎥⎦E=0=αLC/IC1,z(0). (23)

In Fig. 2 (a) we plot, as a function of , the two sides of this equation. The field derivatives on the left-hand side are evaluated by finite differences, while the right-hand side is calculated from Eq. (17). It is clear that the field-correction terms in Eq. (13) are nonzero, and the good agreement between the three curves demonstrates that for both LC and IC they are given by Eq. (17).

The CS contribution does not need additional tests since, as noted above, the contributions to from the zero-field and field-correction terms are identical. However, we reproduce in Fig. 2 (b) the CS curve from Fig. 1 multiplied by a factor , so that the correctness of Eq. (18) can be verified by direct visual inspection. This completes the numerical checks of the -space formula for .

To summarize, we have extended the recently developed band theory of orbital magnetoelectric response to treat crystals under a finite electrical bias. The theory presented in this work may be especially useful in calculations of the second-order magnetoelectric effect defined by Eq. (3). While it is possible in principle to calculate the second derivative of by finite differences, the numerical stability is likely to be improved by taking one of the field derivatives analytically, leaving only one derivative to be performed numerically. We have demonstrated that in order to calculate the total OMP at finite electric field, one may use the same equations (14)–(16) that were previously derived for zero field. This is true even though the individual local-circulation, itinerant-circulation, and Chern-Simons contributions do separately acquire field-correction terms. At present, we are not aware of any simple argument that could have anticipated the exact cancellation of these terms in the expression for the total OMP.

We thank Sinisa Coh for stimulating discussions. The work was supported by NSF under Grants No. DMR-0706493 and No. DMR-1005838. Computational resources were provided by NERSC.

## References

• (1) T. O’Dell, The Electrodynamics of Magneto-Electric Media (North-Holland, Amsterdam, 1970)
• (2) M. Fiebig, J. Phys. D 38, R123 (2005)
• (3) In this work we use SI units. For a discussion of various conventions and choices of units used to define a magnetoelectric tensor, see Ref. coh, and references therein.
• (4) E. Ascher, Philos. Mag. 17, 149 (1968)
• (5) T. H. O’Dell, Philos. Mag. 16, 487 (1967)
• (6) A. M. Essin, J. E. Moore, and D. Vanderbilt, Phys. Rev. Lett. 102, 146805 (2009)
• (7) A. Malashevich, I. Souza, S. Coh, and D. Vanderbilt, New J. Phys. 12, 053032 (2010)
• (8) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010)
• (9) X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Phys. Rev. B 78, 195424 (2008)
• (10) A. M. Essin, A. M. Turner, J. E. Moore, and D. Vanderbilt, Phys. Rev. B 81, 205104 (2010)
• (11) S. Coh, D. Vanderbilt, A. Malashevich, and I. Souza, Phys. Rev. B 83, 085108 (2011)
• (12) I. Souza, J. Íñiguez, and D. Vanderbilt, Phys. Rev. Lett. 89, 117602 (2002)
• (13) D. Ceresoli, T. Thonhauser, D. Vanderbilt, and R. Resta, Phys. Rev. B 74, 024408 (2006)
• (14) X. Wang and D. Vanderbilt, Phys. Rev. B 75, 115116 (2007)
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters