Theory of electric polarization induced by inhomogeneity in crystals
aa r X i v : . [ c ond - m a t . m e s - h a ll ] N ov Theory of electric polarization induced by inhomogeneity in crystals
Di Xiao, Junren Shi, Dennis P. Clougherty, and Qian Niu Department of Physics, The University of Texas at Austin, Austin, Texas 78712, USA Institute of Physics and ICQS, Chinese Academy of Sciences, Beijing 100080, China Department of Physics, University of Vermont, Burlington, Vermont 05405, USA (Dated: October 31, 2018)We develop a general theory of electric polarization induced by inhomogeneity in crystals. Weshow that contributions to polarization can be classified in powers of the gradient of the orderparameter. The zeroth order contribution reduces to the well-known result obtained by King-Smith and Vanderbilt for uniform systems. The first order contribution, when expressed in a two-point formula, takes the Chern-Simons 3-form of the vector potentials derived from the Bloch wavefunctions. Using the relation between polarization and charge density, we demonstrate our formulaby studying charge fractionalization in a two-dimensional dimer model recently proposed.
PACS numbers: 77.22.-d, 75.30.-m, 05.30.Pr, 71.10.Fd, 71.23.An
Electric polarization is a fundamental quantity in con-densed matter physics, essential to any proper descrip-tion of dielectric phenomena of matter. Theoretically, itis well established that only the change in polarizationhas physical meaning and it can be quantified by usingthe Berry phase of the electronic wave functions [1, 2, 3].In practice, the Berry-phase formula is usually expressedin terms of the Bloch orbitals. It has been very successfulin first-principles studies of dielectric properties of oxidesand other insulating materials.While the existing formulation is adequate in periodicinsulators, a theory of polarization for inhomogeneouscrystals would find numerous important applications; forexample, in a class of recently discovered multiferroics,the appearance of electric polarization is always accom-panied by long-wavelength magnetic structures [4, 5, 6].A number of phenomenological and microscopic theorieshave been proposed to understand this magnetically in-duced polarization [6, 7, 8, 9, 10]; however, quantitativestudies of this type of problem still remain in a primitivestate. The fundamental difficulty lies in the fact that theinhomogeneous ordering breaks the translational symme-try of the crystal so that Bloch’s theorem does not apply.In this Letter we present a general framework to calcu-late electric polarization in crystals with inhomogeneousordering. Our theory is based on the elementary relationbetween the change in polarization and integrated bulkcurrent [2, 3]. The latter can be evaluated using the semi-classical formalism of Bloch electron dynamics [11]. Wefind that, in addition to the contribution previously ob-tained for uniform systems [1], the polarization containsan extra contribution proportional to the gradient of theorder parameter. This extra contribution is expressedusing the second Chern form of the Berry curvatures de-rived from the local
Bloch functions. It can also be recastinto a two-point formula, which depends only on the ini-tial and final states, up to an uncertainty quantum afterspatial averaging. We identify this quantum as the sec-ond Chern number in appropriate units. In addition, several general conditions for the inhomogeneity-inducedpolarization to be nonzero are also derived.To demonstrate our theory, we apply our formula tostudy the problem of charge fractionalization in a two-dimensional dimer model recently proposed [12, 13]. Weshow that in this model fractional charge appears as aresult of the ferroelectric domain walls. By using the re-lation between polarization and charge density, we calcu-late the total charge carried by a vortex in the dimeriza-tion pattern and compare it to previous results [12, 13].Our approach has the advantage that it can be easilyincorporated in a band calculation, while previously onerelied on spectral analysis of the Dirac Hamiltonian per-formed in the continuum limit [12, 13].
General formulation .—Suppose we have an insulatingcrystal with an order parameter m ( r ) that varies slowlyin space. We assume that, at least on the mean-field level, m ( r ) can be treated as an external field that couples toan operator in the Hamiltonian H . Thus, we can formallywrite H [ m ( r )]. As was emphasized in previous work [1,2, 3], only the change in polarization P between twodifferent states has meaning, and it is given by [14] P = Z T dt j ( r , t ) , (1)where j ( r , t ) is the bulk current density as the systemadiabatically evolves from the initial state ( t = 0) to thefinal state ( t = T ). In other words, we assume that thetwo states are connected through a continuous transfor-mation of the Hamiltonian H [ m ( r ); λ ] parameterized bya scalar λ with λ (0) = 0 and λ ( T ) = 1.In order to find the current density j ( r , t ), we adoptthe formalism of semiclassical dynamics of Bloch elec-trons [11], which is a powerful tool to investigate theinfluence of slowly varying perturbations on electron dy-namics. Within this approach, each electron is describedby a narrow wave packet localized around r c and k c inthe phase space. If m ( r ) varies smoothly compared tothe width of the wave packet, it is sufficient to study afamily of local Hamiltonians H c [ m ( r c ); λ ] which assumesa fixed value of the order parameter m ( r c ) in the vicinityof r c . Since H c [ m ( r c ); λ ] maintains the periodicity of theunperturbed crystal, its eigenstates have the Bloch form: | ψ n ( k , r c ; λ ) i = e i k · r | u n ( k , r c ; λ ) i , where | u n ( k , r c ; λ ) i is the cell-periodic part of the Bloch functions. Notethat the r c -dependence of | u n ( k , r c ; λ ) i enters through m ( r c ). We can then expand the wave packet using theselocal Bloch functions. For simplicity, in the followingderivation we shall confine ourselves to the case of non-degenerate bands and hence omit the band index n .It has been previously shown that the wave packet cen-ter satisfies the following equations of motion (hereafterthe subscript c on k c and r c is dropped) [11]˙ r α = ∇ kα ε − Ω krαβ ˙ r β − Ω kkαβ ˙ k β − ˙ λ Ω kλα , (2a)˙ k α = −∇ rα ε + Ω rrαβ ˙ r β + Ω rkαβ ˙ k β + ˙ λ Ω rλα , (2b)where ε is the electron energy and we have introducedthe notation ∇ kα = ∂/∂k α and ∇ rα = ∂/∂r α . Summationover repeated indices is implied throughout our deriva-tion. Here, Ω is the Berry curvature obtained from thevector potential A derived from | u ( k , r , λ ) i . For exam-ple, A kα = h u | i ∇ kα | u i , A rα = h u | i ∇ rα | u i , (3)Ω krαβ = ∇ kα A rβ − ∇ rβ A kα . (4)Other Berry curvatures are similarly defined. It is note-worthy that although the vector potential A depends onthe phase choice of the wave function | u ( k , r , λ ) i , theBerry curvature Ω is a well-defined gauge-invariant quan-tity in the parameter space ( k , r , λ ).We now turn to the derivation of P using Eq. (1). The electronic contribution to polarization is given by P = − e Z BZ d k Z T dt D ( k , r ) ˙ r , (5)where − e is the electron charge, and D ( r , k ) is the elec-tron density of states, which is modified from its usualvalue of 1 / (2 π ) d in the presence of the Berry curvature, D ( k , r ) = (1 + Ω krαα ) / (2 π ) d [15].We can solve ˙ r α from Eq. (2) then insert it into Eq. (5).Collecting terms proportional to ˙ λ and keeping those upto first order in the gradient, we obtain [16] P = P (0) + P (1) , (6)where P (0) is the zeroth order contribution P (0) α = e Z BZ d k (2 π ) d Z dλ Ω kλα , (7)and P (1) is the first order contribution P (1) α = e Z BZ d k (2 π ) d Z dλ (cid:16) Ω krββ Ω kλα − Ω krαβ Ω kλβ + Ω kkαβ Ω rλβ (cid:17) . (8) TABLE I: Comparison between P (0) and P (1) Two-point formula Uncertain quantum P (0) Chern-Simons 1-form First Chern number P (1) Chern-Simons 3-form Second Chern number
These are the central results of this work. We note that P (0) has been obtained by King-Smith and Vanderbiltfor uniform systems [1], whereas P (1) , being proportionalto the gradient of m ( r ), only exists in inhomogeneouscrystals.Two remarks are in order: firstly, although in theabove derivation we have assumed an inhomogeneous or-der parameter, it is obvious that our theory is also ap-plicable when the system is subject to a perturbation ofa spatially-varying external field; secondly, we have onlyconsidered the electronic contribution to P here. Whencomparing with experiment, one should also include theionic contribution, which is relatively easy to calculatebecause of its classical nature. Two-point formula .—We first show that P (1) has thedesired property that it depends only on the initial andfinal states. The gauge-invariance of Eq. (8) allows us toevaluate it with any gauge choice. In order to carry outthe integration over λ , we choose the path-independentgauge by requiring that the phase difference between | u ( k , r , λ ) i and | u ( k + G , r , λ ) i does not depend on λ ,where G is a reciprocal lattice vector [3]. Under thisgauge, Eq. (8) can be recast as [17] P (1) α = e Z BZ d k (2 π ) d (cid:16) A kα ∇ rβ A kβ + A kβ ∇ kα A rβ + A rβ ∇ kβ A kα (cid:17)(cid:12)(cid:12)(cid:12) . (9)We recognize that the integrand in the above equation isnothing but the Chern-Simons 3-form.However, we have paid a price for performing the λ -integration; namely, the spatially averaged polarization h P (1) α i = (1 /V ) R d r P (1) α resulting from this two-pointformula (9) can only be determined modulo a quantum.To find the size of the quantum, we consider a cyclicchange in λ . Let us now assume that the order parame-ter m ( r ) is periodic in r . The integral in Eq. (8) (aftera spatial integration) over a closed manifold spanned by( k α , k β , r β , λ ) is an integer called the second Chern num-ber [18]. Since Eq. (9) does not track the evolution of λ ,there is no information of how many cycles λ has gonethrough. This is the reason why h P (1) α i using Eq. (9)can only be determined modulo a quantum. Assuming m ( r ) depends on y , we obtain the quantum for P (1) x ina three-dimensional system:∆ h P (1) x i = el y a z , (10)where l y is the period of m ( y ) and a z is the lattice con-stant along ˆ z .Similarly, the zeroth order contribution P (0) can alsobe cast into a two-point formula and the uncertain quan-tum is given by e/ ( a y a z ) [1]. First-principles calculationsshow that in real materials P (0) is usually smaller thanthis quantum. Hence the ratio between P (0) and P (1) isroughly on the order of l y /a y . The similarities between P (0) and P (1) are summarized in Table I. Minimal conditions for a finite P (1) .—We now eval-uate Eq. (8) using a particular path of λ . We write H [ m ( r ); λ ] = H [ λ m ( r )] so that λ acts like a “switch”of the order m ( r ), i.e., when λ = 0 the system is order-less and when λ = 1 the order is fully developed. Usingthe relation ∇ rα = ∇ rα m µ ∇ mµ and ∇ λ = ( m µ /λ ) ∇ mµ , wecan recast Eq. (8) as P (1) α = em µ ∇ rβ m ν Z BZ d k (2 π ) d Z dλλ (cid:16) Ω kmαµ Ω kmβν − Ω kmαν Ω kmβµ + Ω kkαβ Ω mmµν (cid:17) . (11)As we shall see below, this equation is very useful inassessing the general properties of P (1) .Beside having the crystal be inhomogeneous, there arethree general conditions for P (1) to be nonzero accord-ing to Eq. (11): (i) the system must be two-dimensionalor higher; (ii) the order parameter m ( r ) must have twoor more components; and (iii) the wave function mustdepend on four or more independent parameters. Theseconditions can be obtained by realizing that the inte-grand in Eq. (11) is actually the second Chern 4-formΩ ∧ Ω given in its local expression with respect to the co-ordinates ( k α , k β , m µ , m ν ). It is antisymmetric in k α and k β , and in m µ and m ν , hence condition (i) and (ii). Con-dition (iii) follows from the fact that all 4-forms vanishidentically in three or less dimensions. Based on con-dition (iii) we can further deduce that dim( H ) >
2. Ifdim( H ) = 2, H has four components. However, sinceshifting and scaling energy has no effect on wave func-tions, the wave function can depend on only two indepen-dent parameters (for example, the spherical coordinateson a 2-sphere S ) and P (1) vanishes in this case. Thisset of conditions puts powerful constraints on possiblemicroscopic models that display finite P (1) . Conditions(i) and (iii) can also be obtained directly from Eq. (8).Let us consider a two-dimensional “minimal” modeland assume that both the space of m ( r ) and coordinatespace are two-dimensional. Because of its antisymmet-ric properties, we can write the integrand of Eq. (11) as ǫ αβ ǫ µν χ . Then Eq. (11) takes the following form P (1) = eχ [( ∇ · m ) m − ( m · ∇ ) m ] , (12)Here χ , as a function of m ( r ), can be spatial dependent.Interestingly, if we identify m ( r ) with the magnetiza-tion order parameter M ( r ), the above result is consistentwith the Landau-Ginzburg theory of polarization inducedby spiral magnetic ordering [7]. However, our result (12) ! /m Q πππ πππ πππ FIG. 1: (color online). The total charge carried by the ferro-electric vortex domain wall m ( r ) e iθ = m x + im y as a functionof ∆ /m , where ∆ is the staggered sublattice potential, and m is the dimerization order parameter. As ∆ increases, thedifference between the result from the continuum limit (reddashed line) and that based on the band calculation (bluesolid line) becomes significant. The insert shows the two-dimensional dimerized square lattice with π -flux per plaque-tte in the absence of the vortex. The hopping amplitude is t (1 ± m x,y ) along the x, y -direction. is a direct consequence of the minimal dimensionality andwe did not invoke any symmetry analysis. For higher di-mensions, one will have to carry out a careful symmetryanalysis of the magnetic groups of the crystal [19]. Degenerate bands .—So far, our derivation is for non-degenerate bands. The generalization to degeneratebands is straightforward [20, 21]. The vector potentialand Berry curvature become matrix-valued and are de-fined by ( A a ) mn = h u m | i ∇ a | u n i , (13)Ω ab = ∇ a A b − ∇ b A a − i [ A a , A b ] , (14)where a, b ∈ ( k , r , λ ) and | u m i and | u n i are degeneratebands. We then need to take the trace of Eqs. (7) and(8) for the zeroth and first order contributions to P . Thetwo-point formula in Eq. (9) also takes the non-AbelianChern-Simons form. Fractional charge .—To demonstrate our theory, weconsider the problem of charge fractionalization in a re-cently proposed two-dimensional dimer model [12, 13],shown schematically in the inset of Fig. 1. Introducing γ i = σ i ⊗ σ z , γ = ⊗ σ x and γ = ⊗ σ y , we can writethe Hamiltonian as H = h α γ a , where h = t (cos k x , cos k y , ∆ , m x sin k x , m y sin k y ) , (15) t ∆ is the staggered sublattice potential, t (1 ± m x ) and t (1 ± m y ) are the dimerized hopping amplitudes along x and y direction. We choose the Landau gauge so thatthe effect of the π flux is represented by alternating signsof the hopping amplitudes along adjacent rows. It turnsout that this model is a minimal one satisfying all ourthree conditions: (i) it is two-dimensional; (ii) the orderparameter m = ( m x , m y ) has two components; and (iii) h (after scaling) can be mapped onto a unit sphere S with four independent spherical angles.It can be verified that the energy spectrum of thisHamiltonian consists of two doubly degenerate levels;therefore, the non-Abelian formalism is necessary. TheBerry curvature has SU(2) symmetry [18, 22, 23]; hence, P (0) always vanishes since the non-Abelian version ofEq. (7) has vanishing trace. Thus, we will only consider P (1) in what follows.Suppose there is a vortex in the dimerization pattern:namely, m x + im y = m ( r ) e inθ . According to Eq. (12)together with the fact that ρ ( r ) = − ∇ · P , this ferroelec-tric vortex domain wall will carry a polarization chargeof Q = R d r ρ ( r ) = nm R π dθ χ [7], which is in generalfractional.To compare with previous results, we shall first eval-uate χ in the continuum limit. Expanding the Hamilto-nian around the Dirac point ( π/ , π/ χ = 32 n Z d k (2 π ) Z dλ ∆ λ ( k + m λ + ∆ ) / . (16)Since at large k the integrand decays as k − , we canextend the integration range of k to infinity and obtain χ = n πm (1 − ∆ √ ∆ + m ) , (17)and the total charge carried by the vortex is given by Q = n e − ∆ √ ∆ + m ) , (18)where n is the winding number. This result agrees withthe spectral analysis of the Dirac Hamiltonian [12, 13].The above derivation provides a simple picture ofcharge fractionalization in this type of system: it is adirect consequence of the ferroelectric domain wall, andthe breaking of the sublattice symmetry (∆) allows itto be irrational. A detailed report including both 1Dand 2D cases will be reported elsewhere. We also calcu-late the total charge based on a band calculation usingEq. (15), shown in Fig. 1. As ∆ increases, the devia-tion between the band calculation and continuum limitbecomes significant.In summary, we have developed a general theory ofpolarization induced by inhomogeneity in crystals. Ourresult lays the foundation for quantitative studies of thistype of problem. In connection to multiferroics, the min-imal conditions for a finite P (1) point to general direc-tions to aid in the search for microscopic models. Inaddition, we have illustrated our theory by showing that the fractional charge in certain models can be understoodas the polarization charge accompanying ferroelectric do-main walls.DX thanks D. Culcer for useful discussions. DX wassupported by the NSF (DMR-0404252/0606485), JRSby the NSF of China (No. 10604063), DPC by DARPA(No. MDA0620110041), and QN by the Welch Foun-dation, DOE (DE-FG03-02ER45958), and the NSF ofChina (No. 10740420252). [1] R. D. King-Smith and D. Vanderbilt, Phys. Rev. B ,1651 (1993).[2] R. Resta, Rev. Mod. Phys. , 899 (1994).[3] G. Ortiz and R. M. Martin, Phys. Rev. B , 14202(1994).[4] T. Kimura, T. Goto, H. Shintani, K. Ishizaka, T. Arima,and Y. Tokura, Nature , 55 (2003).[5] N. Hur, S. Park, P. A. Sharma, J. S. Ahn, S. Guha, andS.-W. Cheong, Nature , 392 (2004).[6] G. Lawes, A. B. Harris, T. Kimura, N. Rogado, R. J.Cava, A. Aharony, O. Entin-Wohlman, T. Yildirim,M. Kenzelmann, C. Broholm, et al., Phys. Rev. Lett. , 087205 (2005).[7] M. Mostovoy, Phys. Rev. Lett. , 067601 (2006).[8] H. Katsura, N. Nagaosa, and A. V. Balatsky, Phys. Rev.Lett. , 057205 (2005).[9] I. A. Sergienko, C. Sen, and E. Dagotto, Phys. Rev. Lett. , 227204 (2006).[10] J. Hu, arXiv:0705.0955.[11] G. Sundaram and Q. Niu, Phys. Rev. B , 14915 (1999).[12] B. Seradjeh, C. Weeks, and M. Franz, arXiv:0706.1559.[13] C. Chamon, C.-Y. Hou, R. Jackiw, C. Mudry, S.-Y. Pi,and A. P. Schnyder, arXiv:0707.0293.[14] Strictly speaking, there is still an ambiguity in the defi-nition of P given here because j = ∂ P /∂t + ∇ × M alsocontains a contribution from the magnetization current.As a result, P is only defined up to a divergence-freefield. One can of course fix the gauge by imposing that P only have a longitudinal component. However, this isnot a critical issue because the ambiguity can be removedby spatial averaging.[15] D. Xiao, J. Shi, and Q. Niu, Phys. Rev. Lett. , 137204(2005).[16] The integral of terms that do not contain ˙ λ is given by − e R BZ d k R T dt [ ∇ kα ε + (Ω krββ ∇ kα ε − Ω krαβ ∇ kβ ε + Ω kkαβ ∇ rβ ε )].The integral of ∇ kα ε over the entire Brillouin zone obvi-ously vanishes. After integration by parts and making useof the Bianchi identity ∇ kα Ω rkββ + ∇ kβ Ω krαβ + ∇ rβ Ω kkβα = 0,one can show that the integral of the last three terms con-tributes a divergence-free part, which can be discarded.[17] The formal derivation of Eq. (9) is lengthy and will bereported elsewhere. In the simplest case where A is well-defined everywhere in ( k , r , λ ), one can prove Eq. (9)by using Stokes’ theorem R M Ω ∧ Ω = R M d ( A ∧ d A ) = R ∂ M A ∧ d A .[18] J. E. Avron, L. Sadun, J. Segert, and B. Simon, Phys.Rev. Lett. , 1329 (1988).[19] A. B. Harris, Phys. Rev. B , 054447 (2007).[20] D. Culcer, Y. Yao, and Q. Niu, Phys. Rev. B , 085110 (2005).[21] R. Shindou and K.-I. Imura, Nucl. Phys. B , 399(2005).[22] R. Shankar and H. Mathur, Phys. Rev. Lett. , 1565 (1994).[23] E. Demler and S.-C. Zhang, Ann. Phys.271