When does Magnetization Depend on the Boundary?
In this paper, we investigate in general how thermodynamic quantities such as the polarization, magnetization and the magneto-electric tensor are affected by the boundaries. We show that when the calculation with periodic boundary conditions does not involve a Berry's phase, the quantity in question is determined unambiguously by the bulk, even in the presence of gapless surface states. When the calculation involves a Berry's phase, the bulk can only determine the quantity up to some quantized value, given that (i) there are no gapless surface states, (ii) the surfaces do not break the symmetries preserved by the bulk, and (iii) the system is kept at charge neutrality. If any of the above conditions is violated, the quantity is then determined entirely by the details at the boundary. For components of response functions such as the magneto-electric tensor that are entirely determined by the bulk, the value of those components should connect smoothly to finite frequency and momentum.
magnetization.tex
%% ****** Start of file template.aps ****** %
%%
%%
%%   This file is part of the APS files in the REVTeX 4 distribution.
%%   Version 4.0 of REVTeX, August 2001
%%
%%
%%   Copyright (c) 2001 The American Physical Society.
%%
%%   See the REVTeX 4 README file for restrictions and more information.
%%
%
% This is a template for producing manuscripts for use with REVTEX 4.0
% Copy this file to another name and then work on that file.
% That way, you always have this original template file to use.
%
% Group addresses by affiliation; use superscriptaddress for long
% author lists, or if there are many overlapping affiliations.
% For Phys. Rev. appearance, change preprint to twocolumn.
% Choose pra, prb, prc, prd, pre, prl, prstab, or rmp for journal
%  Add 'draft' option to mark overfull boxes with black boxes
%  Add 'showpacs' option to make PACS codes appear
%  Add 'showkeys' option to make keywords appear
\newcommand{\ba}{\begin{array}}
\newcommand{\ea}{\end{array}}
\newcommand{\bean}{\begin{eqnarray*}}
\newcommand{\eean}{\end{eqnarray*}}
\newcommand{\rm}[1]{\mathrm{#1}}
% my commands:
\newcommand{\Dslash}{\mathcal{D} \!\!\!\!\slash\,}
\newcommand{\pslash}{\partial \!\!\!\slash}
\newcommand{\slsh}[1]{#1 \!\!\! \slash}
\newcommand{\ts}{\v{c}}
\newcommand{\cc}{\v{c}}
\renewcommand{\c}{\'{c}}
\newcommand{\tj}{\'{c}}
\newcommand{\z}{\v{z}}
\newcommand{\s}{\v{s}}
\newcommand{\ep}{\epsilon}
\newcommand{\p}{\partial}
\newcommand{\tr}{\mathrm{tr}}
\newcommand{\Tr}{\mathrm{Tr}}
\renewcommand{\O}{\mathcal{O}}
\renewcommand{\d}{\partial}
\newcommand{\dd}{\,d}
\newcommand{\tg}{\,\mathrm{tg}\:}
\newcommand{\ctg}{\,\mathrm{ctg}\:}
\newcommand{\bra}[1]{\left\langle {#1} \right|}
\newcommand{\ket}[1]{\left|  #1 \right\rangle}
\newcommand{\bracket}[3]{\left\langle {#1} \left| {#2} \right| {#3} \right\rangle \;}
\newcommand{\braket}[2]{\left\langle #1 | #2 \right\rangle \;}
\newcommand{\abs}[1]{\left| #1  \right|}
\newcommand{\VEV}[1]{\bracket{0}{#1}{0}}
\newcommand{\La}{\mathcal{L}}
\newcommand{\Op}[1]{\mathcal{O}(p^{#1})}
\newcommand{\bmx}[1]{\left(\begin{array}{*{#1}{c}}}
\newcommand{\emx}{\end{array}\right)}
\newcommand{\bmxw}[1]{\renewcommand{\arraystretch}{2}\left(\begin{array}{*{#1}{c}}}
\newcommand{\bmxww}[1]{\renewcommand{\arraystretch}{2.5}\left(\begin{array}{*{#1}{c}}}
\newcommand{\bdet}[1]{\renewcommand{\arraystretch}{1.2}
	\left|\begin{array}{*{#1}{c}}}
\newcommand{\edet}{\end{array}\right|\renewcommand{\arraystretch}{1}}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\EV}[1]{\left\langle #1 \right\rangle}
\newcommand{\iinf}{\int_{-\infty}^{\infty}}
\newcommand{\ifour}[1][x]{\int \! d^4 \! #1}
\newcommand{\ithree}[1][x]{\int \! d^3 \! #1 \,}
\newcommand{\ifourn}[2][x]{\int {#2} \, d^4 \! #1}
\newcommand{\ithreen}[2][x]{\int {#2} \, d^3 \! #1}
\newcommand{\gu}[1][\mu]{\gamma^{#1}}
\newcommand{\gd}[1][\mu]{\gamma_{#1}}
\newcommand{\bilin}[2][\psi]{\bar{#1}(x){#2}{#1}(x)}
\newcommand{\D}{\mathcal{D}}
\newcommand{\ds}{\displaystyle}
\newcommand{\chancery}{\fontfamily{pzc}\fontseries{m}\fontencoding{T1}
      \fontshape{it}\selectfont}
\newcommand{\dditem}[1]{\item[$\diamond$] \textbf{#1}\hfill}
\newcommand{\ditem}[1]{\item[$\diamond$]}
\newcommand{\bitem}[1]{\item {\textbf{#1}}}
\newcommand{\bit}{\begin{itemize}}
\newcommand{\eit}{\end{itemize}}
\newcommand{\eV}{\ensuremath{\:\mbox{eV}}}
\newcommand{\keV}{\ensuremath{\:\mbox{keV}}}
\newcommand{\MeV}{\ensuremath{\:\mbox{MeV}}}
\newcommand{\GeV}{\ensuremath{\:\mbox{GeV}}}
\newcommand{\fm}{\ensuremath{\:\mbox{fm}}}
\newcommand{\sech}{\ensuremath{\,\mbox{sech}\,}}
\newcommand{\reff}[1]{(\ref{#1})}


\newcommand{\half}{\frac{1}{2}}
\newcommand{\x}{\hat{\imath}}
\newcommand{\y}{\hat{j}}
\newcommand{\zz}{\hat{k}}
\newcommand{\rarrow}{\qquad\Longrightarrow\qquad}
\newcommand{\eab}{\begin{eqnarray}}
\newcommand{\eae}{\end{eqnarray}}
\newcommand{\nnb}{\nonumber}
\def\l{\label}
\def\C{{\cal C}}
\def\D{{\cal D}}
\def\F{{\cal F}}
\def\G{{\cal G}}
\def\H{{\cal H}}
\def\L{{\cal L}}
\def\S{{\cal S}}
\def\K{{\cal K}}
\def\T{{\cal T}}
\def\V{{\cal V}}
\def\Z{{}}
\def\M{{\cal M}}
\def\H{{\cal H}}
\def\O{{\cal O}}
\def\P{{\cal P}}
\def\dF{{\partial\F\over \partial p}}
\def\dS{{\partial\S\over \partial q}}
\def\dT{{\partial\T\over \partial p}}
\def\U{{\cal U}}
\def\W{{\cal W}}
\def\Q{{\cal Q}}
\def\V{{\cal V}}
\def\ldl{\Lambda\partial_\Lambda}
\def\s{{\tt s}}
\def\t{{\tt t}}
%\documentclass[aps,prb,preprint,groupedaddress]{revtex4-1}
\documentclass[aps,prb,twocolumn,superscriptaddress]{revtex4-1}
%\documentclass[aps,prb,twocolumn,groupedaddress]{revtex4}
%\usepackage{subfigure}
\usepackage{graphics}
\usepackage{epsfig}
\usepackage{amssymb}
\usepackage{color} 
\usepackage{subfigure}
\usepackage{ulem}
%\usepackage{feynmp}
%\usepackage[pdftex]{graphicx}
%\DeclareGraphicsRule{*}{mps}{*}{} 

% You should use BibTeX and apsrev.bst for references
% Choosing a journal automatically selects the correct APS
% BibTeX style file (bst file), so only uncomment the line
% below if necessary.
%\bibliographystyle{apsrev}

\begin{document}
%\begin{fmffile}{fgraphs}

% Use the \preprint command to place your local institutional report
% number in the upper righthand corner of the title page in preprint mode.
% Multiple \preprint commands are allowed.
% Use the 'preprintnumbers' class option to override journal defaults
% to display numbers if necessary
%\preprint{}

%Title of paper
\title{When does Magnetization Depend on the Boundary?}

% repeat the \author .. \affiliation  etc. as needed
% \email, \thanks, \homepage, \altaffiliation all apply to the current
% author. Explanatory text should go in the []'s, actual e-mail
% address or url should go in the {}'s for \email and \homepage.
% Please use the appropriate macro foreach each type of information

% \affiliation command applies to all authors since the last
% \affiliation command. The \affiliation command should follow the
% other information
% \affiliation can be followed by \email, \homepage, \thanks as well.
\author{Kuang-Ting Chen}
%\homepage[]{Your web page}
%\thanks{}
%\altaffiliation{}
\affiliation{Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139}
\author{Patrick A. Lee}
%\homepage[]{Your web page}
%\thanks{}
%\altaffiliation{}
\affiliation{Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139}

%Collaboration name if desired (requires use of superscriptaddress
%option in \documentclass). \noaffiliation is required (may also be
%used with the \author command).
%\collaboration can be followed by \email, \homepage, \thanks as well.
%\collaboration{}
%\noaffiliation

\date{\today}

\begin{abstract}
In this paper, we investigate in general how thermodynamic quantities such as the polarization, magnetization and the magneto-electric tensor are affected by the boundaries. We show that when the calculation with periodic boundary conditions does not involve a Berry's phase, the quantity in question is determined unambiguously by the bulk, even in the presence of gapless surface states. When the calculation involves a Berry's phase, the bulk can only determine the quantity up to some quantized value, given that (i) there are no gapless surface states, (ii) the surfaces do not break the symmetries preserved by the bulk, and (iii) the system is kept at charge neutrality. If any of the above conditions is violated, the quantity is then determined entirely by the details at the boundary. For components of response functions such as the magneto-electric tensor that are entirely determined by the bulk, the value of those components should connect smoothly to finite frequency and momentum.

\end{abstract}

% insert suggested PACS numbers in braces on next line
%\pacs{
% insert suggested keywords - APS authors don't need to do this
%\keywords{}

%\maketitle must follow title, authors, abstract, \pacs, and \keywords
\maketitle

% body of paper here - Use proper section commands
% References should be done using the \cite, \ref, and \label commands

\section{Introduction}
%Bulk-edge correspondence is a very important characteristic for topological order. In integer and fractional quantum Hall systems, the different properties of the chiral edge states at the Fermi level signify the topological nature of the bulk states. For $Z_2$ topological band insulators (TBI) defined under time reversal (TR) symmetry, the gapless edge states are even the major evidence for the existence of the TBI state.
Recently there has been a resurging interest in understanding the general orbital magneto-electric (ME) response\citep{Essin2009,Essin2010,Malashevich2010,Coh2011}. This is due to the fact that the isotropic magneto-electric effect or the so-called $\theta$-term, $\mathcal{L}_\theta=(\theta e^2/2\pi h)E\cdot B$ with $\theta=\pi$, is suggested to describe the the time reversal invarint (TRI) topological band insulator (TBI) in three spatial dimensions\cite{Qi08} (3d). Usually the signature of the 3d TBI is the existence of an odd number of Dirac cones on the boundary surfaces, when the time reversal symmetry (TRS) is preserved.\citep{Hasan2010} When the surface states are gapped out by breaking the TRS locally on the boundaries, a half integer Quantum Hall effect will take place, and give rise to a quantized bulk magneto-electric response.\cite{Qi08} It is later shown that this isotropic response is only a part of the more general anisotropic orbital ME tensor defined in the bulk\cite{Essin2010}:
\beq
\alpha_{ij}=\alpha_{\rm \theta}\delta_{ij}+\alpha_{{\rm 3d}ij},
\eeq
we define $\alpha_\theta=\frac13\Tr\alpha_{ij}$ and $\alpha_{\rm 3d}$ is therefore traceless. $\alpha_{ij}$ describes either the orbital magneto-polarizability (OMP) or the orbital electric susceptibility (OES):
\beq
 \alpha_{ij}={\rm d}P_i/{\rm d}B_j={\rm d}M_j/{\rm d}E_i;
\eeq
OMP and OES are equal via a Maxwell relation. 
 
One peculiarity of the ME tensor is that $\alpha_\theta$ is only determined up to integer multiples of $e^2/h$ by the bulk band structure.\cite{Malashevich2010,Essin2010} The specific value of $\alpha_\theta$ depends on the details at the boundary. From the polarization perspective, if we attach an integer quantum Hall (IQH) layer with filling $\nu=\pm 1$ respectively on the top and bottom surface of a cylinder, the orbital magneto-polarizability (OMP), i.e., ${\rm d}P/{\rm d}B$, along the axis of the cylinder will change by $e^2/h$, due to the density locking to the magnetic field of the top and the bottom IQH layer. From the magnetization perspective, if we attach an IQH layer with filling $\nu=1$ on the side surfaces of the cylinder, the orbital electric susceptibility (OES), ${\rm d}M/{\rm d}E$, will also change by $e^2/h$, due to the extra Hall current flowing on the surface in response to the electric field. Either way, the ME effect is only determined up to an integer multiple of $e^2/h$. In the bulk, this ambiguity corresponds well to the fact that $\theta$ as a coefficient in front of $(E\cdot B)$ is an angle only defined up to integer multiple of $2\pi$, because $\int E\cdot B {\rm d}^3 x{\rm d}t$ is quantized. $\alpha$ is odd under TRS, but this ambiguity makes it possible for $\alpha_\theta$ not to vanish with TRS preserved, as $\alpha_\theta=\pm e^2/2h$ is differed from each other by $e^2/h$. $\alpha_{\theta}=0$ and $\alpha_{\theta}=e^2/2h$ then describe two different insulating states of matter under TRS.

 However, the above properties raise some questions. From the polarization perspective this ambiguity from the bulk is not so surprising, since the zero field polarization is already ambiguous with periodic boundary conditions.\cite{King93,Zak1989} From the magnetization perspective, however, this ambiguity is a bit more puzzling, because one commonly regards magnetization as a bulk property.\cite{Thonhauser05,Ceresoli06} In particular, with periodic boundary conditions there seems to be no reason to expect any ambiguity in the magnetization, whereas the ambiguity in the polarization is easily understood. 

Before answering this rather specific question, we note that there is a even more general one: to what extent are thermodynamic quantities such as polarization, magnetization, and ME response determined by the bulk? Unlike the conventional thermodynamic quantities which are entirely independent of the boundary, we have already seen that boundary sometimes plays a role. How can we tell when will a thermodynamic variable depend on the boundary and when will it not?

In the following, we will discuss case by case from the ground state polarization, orbital magnetization, to the magneto-electric tensor. We will argue through Gedanken experiments that some of them depends on the boundary while others don't. We will verify our argument with numerical simulations. By matching the observations with our previous calculation done with periodic boundary conditions, we can then directly tell from the calculation with periodic boundary conditions how different thermodynamic quantities depend on boundaries.
 
 
\section{Ground state polarization}
The ground state polarization is given by the following formula with periodic boundary conditions:\cite{King93}
\beq
P=-ie\int_{BZ}\frac{{\rm d}^d k}{2\pi^d}\sum_{\alpha(k)\in occ}\bra{\alpha(k)}\frac{\p}{\p k}\ket{\alpha(k)}.
\eeq
In one spatial dimension (1d), the polarization is defined modulo $e$ with periodic boundary conditions: $P=P_0+ne$, with $n$ an integer. This corresponds to the observation that with periodic boundary conditions, we can move every electron to the next unit cell and return to the original state, while the two states should by definition have polarization differed by $e$. With two ends, the polarization will take one specific value, depending on the number of charges we put at the two ends.

However, if there are zero modes at the two ends, the polarization is then ambiguous, as theoretically we can consider superposition of states of different occupancy of the zero modes. The bulk value of the polarization thus depends entirely on the boundary. 

 

In 3d it is a bit more interesting. For simplicity let us assume the system sits on a cubic lattice of size $a$. Now the bulk formula has an ambiguity of $e/a^2$, which also corresponds well to the fact that we can move every electron to the next unit cell and return to the same state. However, with boundary surfaces the situation becomes quite different. Consider a capacitor setup. We are allowed to put any number of charges on each of the opposing surfaces, resulting in a change of the polarization in units of $e/A$ (A is the total surface area). In the thermodynamic limit, we can put any finite density of charges on the surface, and the polarization in the bulk can take any value. Our bulk formula is thus no longer valid. To accommodate the charge on the surface, however, the system needs to either be in a metallic state near the boundary, or to break the lattice translation symmetry in the two in-plane directions. If neither condition is satisfied, then we can only add an integer number of electrons per unit cell, and the bulk formula is recovered, with the remaining ambiguity determined by the surface. 

How can the bulk formula become invalid? We note that the ground state polarization can be understood as a Berry's phase when one adiabatically turns on the electric field. Firstly, in order for the Berry's phase to make any sense, the system has to be gapped. This is the reason why a metallic surface can render the bulk formula invalid. Secondly, if we break the lattice translation symmetry in the two directions perpendicular to the electric field, we can no longer integrate over the momentum in those directions but should instead sum over a large number of sub-bands labelled by the remaining momentum along the direction of the electric field. The polarization will have an ambiguity of $e/A$ in this case. This is different from the conventional thermodynamic quantity, which will require a symmetry breaking in the bulk to change its value. The Berry's phase is thus a rather fragile thermodynamic quantity.


\section{Ground state orbital magnetization}
It is not immediately obvious that the orbital magnetization is independent of the boundary. In the bulk the operator $\hat M\propto({\bf r\times v})$ is ill-defined with periodic boundary conditions, and seems to be growing as one goes near the boundary. Indeed, when one numerically compute $\langle\hat M\rangle$ summing over the local orbitals, there is a finite contribution from the boundary orbitals, which renders the total orbital magnetization different from the naive bulk value.~\cite{Ceresoli06} Nevertheless, it has been shown\cite{Ceresoli06} that the boundary contribution is in fact independent of the details at the boundary via the use of local Wannier functions, in an insulator with zero Chern number. 

However, in a Chern insulator, a local Wannier function can not be found\cite{Thonhauser06,Thouless1984}, because the Bloch functions cannot be periodic and smoothly defined over the Brillouin Zone. To see that even in this case the orbital magnetization is still independent of the boundaries, we can consider the following setup:

Suppose we have an insulator with a non-vanishing Chern number in two dimensions. Let us imagine putting an auxiliary layer of insulator on top, with an opposite Chern number, without any interaction with the original one. The new insulator as a whole is then of total Chern number zero. We can therefore make a local Wannier orbital, by a linear combination of orbitals from the two layers.\cite{Soluyanov2011} The argument then goes through for the insulator as a whole, and the total orbital magnetization should be independent of the boundary. Now since there is no interaction between the two layers, the total magnetization is just the sum of the magnetization of the original insulator and the auxiliary insulator. We now consider a particular boundary condition, where the two insulators couples to independent boundary terms that do not interact with each other as well. Let us only vary the boundary terms that couple to the original insulator. The total magnetization cannot change, and neither the contribution from the auxiliary insulator. We thus have to conclude that even for a Chern insulator, the orbital magnetization is independent of the boundaries. 

From this abstract point of view, the generalization to Chern insulators seems rather trivial. However, the presence of gapless chiral edge states may cause one to worry. Suppose we can gate the material to supply a constant chemical potential, what will happen if we turn up the electric potential on the edge? Will the edge current decrease because fewer edge states are occupied, or will it stays the same as required for the bulk magnetization not to change?

\begin{figure}[htb]
	\centering
	\subfigure[]{\includegraphics[width=3.5cm]{chern_1.eps}}
	\subfigure[]{\includegraphics[width=3.5cm]{chern_2.eps}}
   \subfigure{\includegraphics[width=1cm]{legend1.eps}}
	\caption{We take our Hamiltonian to be \\$H=\sum_n c_n^\dag(\tau_z-i\tau_x)c_{n+\hat x}+c_n^\dag(\tau_z-i\tau_y)c_{n+\hat y}+mc_n^\dag\tau_zc_n+h.c.$,\\ where $\tau's$ are the Pauli matrices. At half filling with $m=1.5$, the band carries a Chern number $C_1=1$. If we set the chemical potential $\mu=0$, the ground state has no magnetization. We put the Hamiltonian on a $10\times 10$ lattice, and take open boundary conditions in both directions. The current on the vertical links is plotted. We relate the current to the magnetization by $I^b=\ep^{ab}\p_aM$, and take the magnetization at the middle to represent the bulk magnetization. (a) $\mu=0.5$. As expected, some edge states are occupied and give rise to a bulk magnetization. (b) If we set $\mu=0$ but locally apply an electric potential $V=-0.5$ to the first two rows at the boundary, the edge states are again occupied. However, in the region next to those layers, a counterpropagating current takes place. The bulk magnetization remains zero (barring some finite size effect).
}
	\label{chern}
\end{figure}  

We do a straightforward numerical simulation to resolve this paradox. The result is shown in Fig. \ref{chern}. We can see that while shifting the overall chemical potential creates circulating currents, altering the electric potential locally at the edge does not change the bulk magnetization. If we look closer, while the current right at the edge is changed, there is a counter-propagating current near the edge, which keeps the total current localized near one edge constant. The counter-propagating current is just the integer quantum Hall response to the electric potential gradient. This bulk quantum Hall current exactly compensates for the current carried by the now-unoccupied edge states, and leaves the bulk magnetization insensitive to the change of the potential local near the edge.
 

%We mention in passing that one related problem is to consider the orbital susceptibility to a Zeeman field.\cite{yichou} It seems tempting to consider such effects in TBIs in 3D, as the surface states are the only states at low energy that can contribute. However, similar to what we have seen above, the contribution need not only come from states near the Fermi energy. The underlying bands will also contribute, so the total orbital Zeeman susceptibility can not be used as a signature to the existence of such surface states.


A very similar puzzle arises in the $S_z$ conserved spin Hall insulator. On the edge there are counter-propagating TR-paired edge states. When we apply a uniform Zeeman field $H_z$, there will be a net circulating current from the edge states. We can therefore deduce a bulk orbital magnetization response to the Zeeman field. We call this the orbital-Zeeman susceptibility. However, one can locally break the $S_z$ conservation together with the TR symmetry near the edge, to gap out the edge states. In this case, will there still be a bulk magnetization response to the Zeeman field? 


The numerical result is shown in Fig. \ref{spinhall}. Here we can see that even though the edge states are gapped out by the local perturbations, the total current flowing near the edge remains the same. The local perturbation transfers the current from the states at the Fermi level, to the occupied bands. In the end, while local properties can affect the gapless states, the total current near the edge in the is unaffected.

\begin{figure}[htb]
	\centering
	\subfigure[]{\includegraphics[width=4cm]{sh_1.eps}}
	\subfigure[]{\includegraphics[width=4cm]{sh_2.eps}}
	\subfigure[]{\includegraphics[width=3.5cm]{sh_3.eps}}
	\subfigure[]{\includegraphics[width=3.5cm]{sh_4.eps}}
\subfigure{\includegraphics[width=0.9cm]{legend2.eps}}	
	\caption{We now think of the previous model as from spin up electrons and pair it with its time reversal. We applied a unifrom Zeeman field $\delta H_z=0.2\sum_n c_n^\dag S_zc_n$. (a) We plot the eigenstate energies in ascending order. The edge states live inside the gap. (b) By applying a time-reversal as well as $S_z$ symmetry breaking term near the boundary $\delta H=\sum_{\rm n\in edge}c^\dag_n S_x c_n$, we can gap out the edge states. (c)-(d) We look at the current on the vertical links. While the current distributes slightly differently with or without the symmetry breaking term at the edges, the contributions to the bulk magnetization are identical.}
	\label{spinhall}
\end{figure}  

We therefore conclude that the orbital magnetization, as well as the orbital-Zeeman susceptibility is independent of the boundary for an insulator. While the circulating current may be carried by the edge states, the total amount is entirely insensitive to the local boundary conditions. One can understand this from a calculation with periodic boundary conditions: the magnetization is calculated as an energy density in a magnetic field. The total energy, unlike the Berry's phase, is a truly extensive property, so that the boundary contribution is irrelevant in the thermodynamic limit. The energy density in the bulk is thus entirely independent of the boundaries far away enough, whether there are gapless states or not.

%In addition we can also comment on the problem raised in Refs.~\cite{dunno}. The situation in 3D is similar. While the surface state will carry a circulating current when a Zeeman field is applied, it is not the only contribution. The total current is constrained by the bulk magnetization. Unfortunately, without specific knowledge to the current carried by the occupied bands, we cannot characterize the topological insulator by its bulk orbital-Zeeman susceptibility.

\section{magneto-electric effect} %Recently, the general theory of orbital magneto-electric effect is developed, to relate to the $\theta$-term, which characterizes the TI in 3D when $\theta=\pi$. The key new finding is that the diagonal response of either the magnetization to the electric field, or the polarization to the magnetic field, is ambiguous in the bulk, up to any units of $e^2/h$. In the case of polarization, the ambiguity is the same as the intrisic ambiguity of the ground state polarization in the bulk, if we consider the enlarged magnetic unit cell. For the magnetization, one can easily see that on a cylinder, the orbital electric susceptibility (OES) will change by $ne^2/h$ when one attaches an interger quantum Hall state with filling factor $n$ on the side surface; nevertheless, it is harder to understand this ambiguity from the bulk perspective. Here we would like to address this ambiguity, by understanding which category does the OES belong. 
%Now we consider the orbital magnetization in an electric field. The OES is thought to be related to the OMP via the Maxwell relation. Let us first briefly discuss the ME effect from the OMP perspective. In our previous work\cite{mywork,mybigwork} we explained that with periodic boundary conditions, the OMP can be understood as a Berry's phase proportional to the magnetic flux. Due to the fact that the Berry's phase is defined only modulo $2\pi$, and the fact that the magnetic flux is quantized in units of $h/e$, the OMP is only defined modulo $e^2/h$. When there are boundary surfaces, one can show that the OMP takes one definite value among the ambiguities, according to the details at the boundary. In a slab geometry, the polarization is the surface charge density. Two paired surface IQH layers on the top and the bottom surfaces change the OMP by $\nu e^2/h$, as discussed in the introduction. In a cylinder with top and bottom surface identified, the magnetic flux is no longer quantized. While the Berry's phase is ambiguous, the derivative of the phase with respect to the magnetic field is not. For a given setting we therefore have a definite value. Similarly, an IQH layer on the side of the cylinder will also change the OMP by $\nu e^2/h$: when we turn on a magnetic field, the electric field from Faraday's law will generate a surface Hall current, integrated to give a polarization of $(\nu e^2/h) B$. The OMP thus belongs to class (iia) for a finite size sample.
After the discussion of the polarization and the magnetization and seeing that they are thermodynamic quantities with very different behaviors, it is thus a natural question to ask the same question about the ME tensor; in addition, about how the Maxwell relation can be maintained. Before going into details of the boundary dependence, however, let us first show that the anisotropic part $\alpha_{3d}$ is independent of the boundaries.

In terms of electronic Green's functions and with periodic boundary conditions, we have derived the ME tensor from the OMP perspective, as a Berry's phase in a magnetic field:\cite{Tim2011c}
\bea\label{ompg}
\alpha_{ij}&=&(\alpha_{\rm wzw}+\alpha_{\rm 3d})_{ij},\nonumber\\
{\alpha_{\rm wzw}}_{ij}&=&-\frac{\pi i}6 \ep_{abcd}\Tr^S({\bf g}\p_a{\bf g}^{-1}{\bf g}\p_b{\bf g}^{-1}{\bf g}\p_c{\bf g}^{-1}{\bf g}\p_d{\bf g}^{-1}{\bf g})\delta_{ij};\nonumber\\
{\alpha_{\rm 3d}}_{ij}&=&-\frac{i}6\ep_{abj}\Tr\big(g\p_ig^{-1}g\p_ag^{-1}g\p_bg^{-1}g-h.c.\big).
\eea
The traces include the frequency and momentum integral divided by factors of $(2\pi)$; the symbol $\Tr^S$ denotes the integral and trace in one extra dimension in momentum space, with the original Brillouin zone and a trivial test system as the boundary. While the entire ME tensor is derived as a Berry's phase, $\alpha_{\rm 3d}$ does not depend on the Green's function extended to the extra dimension. Without considering boundaries directly, we can show that $\alpha_{\rm 3d}$ is independent of the boundaries, by showing it extends smoothly to finite frequency and momentum. 

At finite frequency and momentum, the ME response is understood as a term in the effective action which is proportional to $E^i(q,\omega)B^j(-q,-\omega)$. Unlike the uniform ME response however, this term can no longer be understood as OMP or OES, due to the fact that unlike uniform electromagnetic fields, the electric and magnetic fields at finite frequency and momentum are related by Faraday's law. The term nevertheless affects properties of the propagating electromagnetic waves. For our purposes, it suffices to show that the effective Lagrangian is continuous from $q=0$ to $q\rightarrow 0$. At any $q\neq 0$, we can calculate the effective Lagrangian by the conventional diagrammatic method. Calculated in the Appendix, the bubble diagram gives 
\beq
S_{ME}
=-\int\frac{{\rm d^4}q}{(2\pi)^4} B^\ell(q) E^k(-q)\alpha_{{\rm 3d}k\ell}+\mathcal{O}(q).
\eeq
Comparing with Eq.~(\ref{ompg}), we see that $\alpha_{\rm 3d}$ is continuous, whereas $\alpha_{\rm wzw}$ is entirely absent at finite momentum. One might worry that we have missed $\alpha_{\rm wzw}$ in momentum space due to the fact that it is a total derivative in real space, which Fourier transforms to zero and cannot be seen in momentum space. However, one can evaluate the diagram in real space, and it is still absent. Fundamentally this is due to the fact that the conventional perturbation theory is perturbative in orders of the gauge field, which breaks down with uniform field strength. Nevertheless, combining the two calculation, we can still say that $\alpha_{\rm 3d}$ is a bulk property and is independent of the boundaries. $\alpha_{\rm wzw}$, on the other hand, is similar to the polarization: it does depend on the boundary, but when there is no boundary, it presents itself as a Berry's phase. Note that one benefit of using the Green's function is that the separation of the local terms and boundary terms matches exactly how the expression depends on the extra dimension or not. This is not the case if we use the density matrices, either to calculate the same Berry's phase\cite{Tim2011a}, or to calculate a current response to a pumping procedure\cite{Essin2010}. In both calculations the ME tensor naturally separates into two terms, with the first term independent of the energy gap: 
\beq
\alpha=\alpha_{\rm cs}+\alpha_{\rm G};
\eeq
$\alpha_{\rm cs}$ is isotropic, but $\alpha_{\rm G}$ is not traceless. While $\alpha_{\rm G}$ can be uniquely determined by the bulk band structure and is independent of the boundaries, its trace is actually not measurable in the bulk.

Let us now focus at the isotropic part $\alpha_{\rm wzw}$. In terms of polarization in a magnetic field, the ambiguity is no surprise. However, how does the ambiguity of the orbital magnetization in a electric field come about?

One origin of the ambiguity is from the fact that the perturbation of a uniform electric field grows with distance. It therefore naturally depends on the boundary, when there is one. When we consider periodic boundary conditions, however, it becomes less clear.

%Let us now turn to the OES. Which class does the OES belong? In Ref.~\cite{Mala}, the OES tensor is calculated with open boundary conditions in all the spatial directions. {\color{red} They claim the isotropic part of the OES is ambiguous up to $e^2/h$, which implies that the OES should belong to class (ii). However, there is one conceptual difficulty: they derive the formula of the OES in momentum space without any reference to the boundary, which is only possible when the OES belongs to class (iii). (They follow a similar argument to Ref.~\cite{ceresoli} to claim it is the case.) In their calculation, the ambiguity of the OES seems to arise from the fact that the perturbation expansion in the external electric field is intrinsically ambiguous with periodic boundary conditions\cite{sth}; it is unclear how the ambiguity is reltated to the physical ambiguity from the boundary. Furthermore, following their calculation, it seems that the OES with periodic boundary conditions is ambiguous as well. On physical grounds it is very puzzling, because the magnetization in zero field belongs to class (iii) as we have shown. It is very hard to imagine how the ambiguity in the bulk may arise in the electric field, without any boundaries. 

%Here we avoid the question of how to derive an analytic formula for the OES, augmented with possible boundary contributions. We focus instead on the physical question: we know the OES takes one of the possible values given by Ref.~\cite{Mala} with open boundary conditions. Does the OES take the same value with periodic boundary conditions, and does it have any intrinsic ambiguity? In the following, we shall focus on the isotropic part of the OES, which is the problematic part.}

In order to study the OES with periodic boundary conditions, we first have to properly define the magnetization with periodic boundary conditions. Without the current at the boundary, one sensible definition of the magnetization is from the relation $B=H+M$. That is, in the absence of applied current (which generates $H$), the magnetization simply equals the measured magnetic field. Note that with periodic boundary conditions and a finite volume, the magnetic field is quantized, because the total magnetic flux through the sample is quantized in units of $h/e$. In this case we take the perspective that the magnetic field will take the closest quantized value to the magnetization while the magnetization itself is still continuous. %The usual definition $M={\rm d}F/{\rm d}B$ encounters difficulty when we apply an electric field with periodic boundary conditions.\cite{explain}

%In thermodynamic limit this definition agrees with the definition such that $M={\rm d}F/{\rm d}B$, where $F$ is the free energy of the system. The latter definition encounters difficulty when we apply an electric field with periodic boundary conditions, since the system is not in a stationary state. 

%Without an electric field, the magnetization can be thought of as the derivative of the free energy with respective to the magnetic field. With an electric field, this definition becomes tricky. Since the electric field cannot be thought of as from a time-independent potential with periodic boundary conditions, the system is not in a stationary state, and the free energy is therefore undefined. A better way to define the magnetization is to require the magnetic field to be the same as the magnetization in the thermodynamic limit\cite{explain}, when we do not put any current by hand in the system. We therefore define the OES in a periodic system, as the response of the magnetic field to the electric field. 

In our previous work\cite{Tim2011a}, we have shown that in a magnetic field, the $\theta$ term, which characterizes the isotropic part of the OMP, changes the quantization condition of the global electric flux. The ground state of the system thus carries an electric flux of -$(\theta e^2/2\pi h)\Phi_B+ne$, where $n$ is some integer that minimizes the flux. Using $0=D=E+P$, the $\theta$ term thus gives an isotropic orbital magneto-polarization response $\frac{\p P}{\p B}=\frac{\theta e^2}{2\pi h}$. However, this result is valid only when $(\Phi_B\theta e/2\pi h)<1$. In the thermodynamic limit this condition is always violated, and instead $\frac{\p P}{\p B}=0$. %Combining this with the well-known fact that the $\theta$ term is measurable with boundary surfaces, we conclude that the isotropic part of the OMP belongs to class (iib) in the thermodynamic limit. At finite size however it belongs to class (iia).

 Similarly, to see whether the same term contributes to the OES of the system, we would like to investigate whether there is a uniform magnetic field, when we constrain the path integral to have a given average electric field in the same direction. However, the electric field and the magnetic field behave in intrinsically different ways, when we formulated our theory assuming the existence of electric charges and the absence of magnetic monopoles: the quantization of the electric flux can change in the presence of the magnetic field, while the quantization of the magnetic flux is fixed at $(h/e)$. When we apply an electric flux, we can always imagine that the system is a coherent state composed of states with integer electric fluxes. The background magnetic field therefore does not have to be different from zero. Therefore, even at finite size, the $\theta$ term does not give rise to the OES. \textit{The Maxwell relation between the the isotropic OMP and the OES are thus violated.} They are only equal in the thermodynamic limit, where the $\theta$ term gives no contribution for both quantities. In other words, the isotropic OES is better thought of as a bulk-induced surface response, which vanishes when there is no boundary surfaces.

Now let us consider geometry with boundaries in some detail. From the result of Ref.~\cite{Malashevich2010}, we know that with open boundary conditions in all directions, the OES has an ambiguity only determined by specific surface boundary conditions. We have also seen in the introduction that in a cylinder geometry, the ambiguity of the OES can come from the quantized Hall current on the side surfaces. 

\begin{figure}[htb]
	\centering
	\includegraphics[width=8cm]{dmde.eps}
	\caption{Here we plot the calculated OES versus the number of layers in the $z$-direction, with the model described by Eq.~(73) in Ref.~\cite{Qi08} with $\theta=0$, $m=c=1$. We take $\theta=0.5\pi$ at the top and the bottom layer to gap out the edge states. (If we take $\theta=\pm 0.5\pi$ on the two surfaces respectively, the whole system will be a Chern insulator and can no long be kept at charge neutrality without closing the gap in a magnetic field.) We put on an electric field such that the potential difference between the top and the bottom layer is $0.2$. The boxes show the calculated values. The solid curve is a fit by assuming a fixed width $w$ of the surface charges when there is a magnetic field, such that $\frac{{\rm d}M}{{\rm d}E}\propto(1-\frac wn)$. The fit gives gives $\frac{{\rm d}M}{{\rm d}E}=0.50 \frac{e^2}{h}$ in the thermodynamic limit and $w=2.54$. The OES changes sign as expected, when we change to $\theta=-0.5\pi$ instead on the boundary.}
	\label{OES}
\end{figure} 

What if there are no side surfaces? Suppose we take periodic boundary conditions only in two directions to get rid of the side surfaces. Does the OES still have the same ambiguity? One naively would expect the situation to be similar to the case with periodic boundary conditions, due to the absence of the possible circulating Hall currents. However, a more careful argument shows it is not the case. In fact, the system will spontaneously generate a magnetic field, which will then generate surface charge density $\sigma=\pm(\nu+\theta/2\pi)e^2B/h$ via the OMP response, to lower the electric energy. Minimizing the total energy as a function of $B$, we then get $B=M=(\nu+\theta/2\pi)e^2E/h$. While at finite size the total magnetic flux is quantized in units of $h/e$ in this setup, in the thermodynamic limit, the magnetic field will converge to the expected value, in contrast to the situations with periodic boundary conditions where it stays at zero. We have numerically confirmed this result by calculating the magnetization in the electric field, using the momentum space formula for the magnetization, derived in Ref.~\cite{Ceresoli06}, as shown in Fig.~\ref{OES}.

 

Before summing up, let us consider how gapless surface states can alter the ME response. Evidently, if we attach a fractional quantum Hall state on the side of the cylinder, the OES is going to change by a fraction of $e^2/h$.\cite{Brian2011} In general the fraction is quite arbitrary, so in this case the bulk value of the isotropic OES is not valid. This corresponds to the fact that the fractional quantum Hall state has ground state degeneracy. In general, we will therefore expect any gapless surface state will destroy the bulk description of the isotropic ME response.

To sum up, The anisotropic part of the ME tensor $\alpha_{\rm 3d}$ is independent of the boundaries. The isotropic part $\alpha_{wzw}$ depends partially on the boundary. While $\alpha_{\rm 3d}$ is a truly local quantity, $\alpha_{\rm wzw}$ only lives at $q=0$. Corroborating with the fact that both isotropic OES and OMP responses vanish with periodic boundary conditions in the thermodynamic limit, it is better to  think of $\alpha_{\rm wzw}$ as a quantized surface effect induced by the bulk.   

\section{conclusion}
 We have thus gone through polarization, magnetization, and magneto-electric responses and see their dependence on the boundary. %Using our formalism, the quantity expressed in terms of Green's functions in the extra dimensions are the ones that are ambiguous and the ambiguity is further resolved by the boundary. Those quantity are more fragile as a bulk quantity: either symmetry breaking just near the boundary or gapless surface states can destroy the bulk constraint. 
A bulk calculation done with periodic boundary conditions contains enough information to predict how the quantity in question can depend on the boundary, however. In particular, using our formalism described in Ref.~\cite{Tim2011c}, any quantity that does not involve an extension of the Green's function to one extra dimension is independent of the boundary. On the other hand, quantities that requires an extension to extra dimension will depend on the boundary. The bulk can determine its value up to some quantized amount, only when (i) there are no gapless surface states, (ii) surfaces break no symmetry that is required to determine the bulk value with periodic boundary conditions, and (iii) the system is kept at charge neutrality. If any of the conditions are violated, the surface contribution will dominate and render the results obtained with periodic boundary conditions invalid.

\section*{Acknowledgement}
We thank X.G. Wen and N. Nagaosa for insightful discussions. This work is supported by NSF grant DMR 1104498.
 %In this geometry, the electric field can be described by a static electric potential. As long as the potential difference of the two ends is smaller than the band gap, our system remains in the insulating state. We calculate the magnetization using a discretized version of the formula
%\beq\label{magnetization}
%M=\frac{i\ep^{ab}}2\int\frac{{\rm d}^2k}{(2\pi)^2}\sum_{n\in {\rm occ}}\bracket{\p_au_n}{H+E_n}{\p_bu_n};
%\eeq 
%the results are summarized in Fig.~\ref{oes}. In a topological insulator, the OES is not zero, with the ambiguity depending on the details at the boundary. To understand how the magnetization arises in this geometry, we note that in deriving Eq.~\ref{magnetization}, we calculate the energy of the system in the presence of a weak magnetic field. In the magnetic field, the top layer and the bottom layer are charged due to the OMP response, which can also be thought of as from the $\pm(\nu+\frac12)$ quantum Hall effect from the top and the bottom surface. The polarization then naturally gives an energy in the electric field. In the thermodynamic limit, the magnetization will generate a magnetic field with equal magnitude. If we contrast this situation to the case with periodic boundary conditions in every direction, we should understand the bulk manetization response to the electric field as generated from the top and bottom surfaces. The situation is very different from the ground state magnetization of a Chern insulator: in that case, the magnetization is present with or without a boundary, despite that one can understand the magnetization as from the edge states when there is a boundary. In our case, the isotropic OES is only present with boundaries. It is thus a surface effect, but still required by the bulk band structure.


%\section{local measurable responses}
%\textit{Local measurable responses}: Now that we know the isotropic OES belongs to class (iib), it becomes sensible to ask whether one can measure the bulk and the surface contribution independently. Specifically, a ``local" measurement in the bulk of the OES should only measure the bulk contribution, since it does not know about the boundary. One can think of a local measurement as a response at a small momentum. The question then becomes that, how does the OES response extend to finite momentum?

% Physically, we can think of the first term describing the current or charge response from the surface, as we have shown in the previous section. The surfaces do not contribute at finite momentum if it is sufficiently far away. For example, on a cylinder along the $z$ direction with radius $r$, the surface quantum Hall effect gives the isotropic OES response described by $\alpha_{\rm wzw}$. If we consider an electric field at finite $q_z$, the magnetization at the same momentum will decay into the bulk following the Biot-Savart law. In the limit where $r>>1/q_z$, the bulk receives no surface contribution at the momentum. 

%The remaining local contribution can only be from $\alpha_{3d}$. While we derived the whole formula using the Berry's phase procedure, the second term $\alpha_{\rm 3d}$ is defined in the physical momentum space. We therefore expect the traceless part of the OES tensor extends smoothly to finite momentum. The locally measurable part of the OES is given by $\alpha_{\rm 3d}$, which is traceless. The formula of the ME tensor in terms of the Green's function therefore naturally decomposes the tensor in a physical way: the isotropic part is ambiguous; physically it comes solely from the surface and cannot be locally measured in the bulk. The remaining part is locally measurable and is insensitive to boundary conditions. 

% Create the reference section using BibTeX:
%\bibliography{basename of .bib file}
%\end{fmffile}
\bibliography{mybib}
%\bibliographystyle{phjcp}

$\;$
\newpage 
\newpage


%\subsection{Numerical Calculation of the OES:}
%The orbital magnetization in two spatial dimensions is given by
%\beq\label{magnetization}
%M=\frac{i\ep^{ab}}2\int\frac{{\rm d}^2k}{(2\pi)^2}\sum_{n\in {\rm occ}}\bracket{\p_au_n}{H+E_n}{\p_bu_n}.
%\eeq
%We numerically evaluate Eq.~\ref{magnetization} using the density matrices expression of the formula, as in Ref.~\cite{gonze}:
%\beq
%M=-\frac{i\ep^{ab}}2\int\frac{{\rm d}^2k}{(2\pi)^2}\Tr\bigg(\p_a\rho\p_b\rho(1-2\rho)H\bigg);
%\eeq
%$\rho(k)=\sum_{n\in{\rm occ}}\ket{u_{nk}}\bra{u_{nk}}$ is the density matrix.The density matrices are gauge invariant, which allows us to approximate the differentiation by finite difference. To check the validity of our numerical method, we apply it on some random two dimensional time-reversal breaking insulator, and compare the result with a straightforward real space calculation by integrating out the current distribution, using $I^b=\ep^{ab}\p_aM$. The parameters of the model is listed in Table \ref{tab}, and the result is shown in Fig~\ref{comp}.

%\begin{table}[htb]
%\centering
%\begin{tabular}{ |l |p{1.0cm} p{2.2cm} p{2.2cm}| }
%    \hline
%    $(x,y)$ & $E$ & $t_x$ & $t_y$ \\ \hline
%    $(-,-)$ & 0.5 & $2.2\exp(-i0.8\pi)$ &$\exp(-i0.5\pi)$\\ \hline
%    $(-,+)$ & -0.5 & $\exp(i\phi_1)$ &$\exp(-i0.25\pi)$\\ \hline%	$(+,+)$ & 0 & $\exp(i0.3\pi)$ &$\exp(i0.3\pi)$\\ \hline
%	$(+,-)$ & 0 & $2.2\exp(-i0.3\pi)$ & $\exp(i0.45\pi)$\\ \hline
% \end{tabular}
%\caption{\label{tab}Parameters of the tight-binding model. We consider a spinless nearest-neighbor tight-binding model on a square lattice. We double the unit cell in both directions, and add random phases to each of the eight independent hopping amplitude, then collapse the four orbitals in the doubled unit cell to the same position. $(x,y)$ denotes the position of the orbital before the collapse. $t_x$ and $3++t_y$ are hopping amplitudes to the next orbital in the $x$ and $y$ direction, respectively. $\phi_1$ is scanned from $0$ to $2\pi$.}
%\end{table}

%\begin{figure}[htb]
%	\centering
%	\includegraphics[width=8cm]{mcomp.eps}
%	\caption{The orbital magnetization of an time-reversal breaking insulator. The tight-binding model is described in Table~\ref{tab}. We plot the resulting orbital magnetization as a function of $\phi_1$. The solid curve is from the calculation in momentum space, and the boxes are values integrated from the current distribution in real space. In real space we have a square lattice with $20\times 20$ sites, with 4 orbitals on each site. we take the magnetization in the middle to be the bulk value.}
%	\label{comp}
%\end{figure}



%We now calculate the orbital magnetization in a slab geometry, with a small applied electric field. We consider periodic boundary conditions in the $x$ and $y$ direction, and open boundary conditions in the $z$ direction. We apply a small electric field in the $z$ direction by a static electric potential, such that the potential difference at the two ends of the sample is smaller than the energy gap.

%We find the energies and the eigenstates of the 1D hopping problem at every $k_x$ and $k_y$. We then construct the density matrix at every $k$, and use the numerical method above to calculate the orbital magnetization in the $z$ direction. Dividing by the electric field we then get the OES.
 
%We use the model described by Eq.~(73) in Ref.~\cite{QHZ} with $\theta=0$ in the bulk, which describes a 3D topological insulator. (In this model, $\theta=\pi$ decribes a trivial insulator.) We also set $\theta=0.5\pi$ at both the top and the bottom layer to gap out the edge states. We plot the OES as a function of the number of layers in the $z$ direction, as shown in Fig~\ref{OES}. While the finite size effect is evident, the OES goes toward $e^2/2h$ when the size gets large.



\begin{widetext}
\section*{Appendix: ME Effect at Finite Momentum}

In the main text we have argued heuristically that the trace of the ME tensor comes entirely from the surface, and therefore does not contribute at finite momentum. The locally measurable ME tensor is therefore traceless in the $q\rightarrow 0$ limit. We can directly calculate the ME tensor at finite $q$: Fourier transforming and expanding the hopping Hamitonian up to second order of $A^\mu$, we have
\beq
\label{H}
\Delta H=\sum_{k,q}c^\dag_{k+q/2}\p_uH(k)c_{k-q/2}A^\mu(-q)+\frac12\sum_{k,q,q'}c^
\dag_{k+(q+q')/2}\p_u\p_\nu H(k)c_{k-(q+q')/2}A^\mu(-q)A^\nu(-q');
\eeq
$H_k\equiv\sum_{d_i}t_{d_i}\exp(ikd_i)$ is a matrix. Integrating out the electrons, the effective action at quadratic order of $A_\mu$ reads

\beq
S_{\rm eff}=\int\frac{{\rm d^4}q}{(2\pi)^4}\frac i2 A^\mu(q) A^\nu(-q)\Tr\bigg(\p_\mu\p_\nu g^{-1}(k)g(k)+\p_\mu g^{-1}(k)g(k+\frac{q}2)\p_{\nu}g^{-1}(k)g(k-\frac{q}2)\bigg);
\eeq

similarly, the trace includes the integral of energy and momentum divided by $2\pi$. The first term in the trace is from the second term in Eq.~(\ref{H}), usually called the paramagnetic current, and does not have $q$ dependence. To compare with Eq.~(\ref{ompg}), we Taylor-expand the second term to second order in $q$ to get the behavior in the $q\rightarrow 0$ limit: (From here on, we drop the dependence on $k$ to avoid cluttering.)
\beq
S_{\rm eff}\sim\int\frac{{\rm d^4}q}{(2\pi)^4}-\frac i4 q^\lambda q^\sigma A^\mu(q) A^\nu(-q)\Tr\bigg(\p_\mu g^{-1}\p_\lambda g\p_{\nu}g^{-1}\p_\sigma g-\frac12\p_\mu g^{-1}\p_\lambda\p_\sigma g\p_{\nu}g^{-1} g-\frac12\p_\mu g^{-1}g\p_{\nu}g^{-1}\p_\lambda\p_\sigma g\bigg)+\mathcal{O}(q^3);
\eeq

To further simplify the expression, let us now take the Coulomb gauge. In the Coulomb gauge, we have to take either $\lambda$ or $\sigma$ to be in the time direction to have the expression contribute to the ME tensor. Since $\p_i\p_\omega g^{-1}=0$, we can integrate-by-part the time derivative. using $\p_\omega g=-g^2$, and rename the indices $i, j, k$, now running through only the spatial directions, we get
\beq\label{sme}
S_{ME}\sim\int\frac{{\rm d^4}q}{(2\pi)^4}\frac i2 \omega q^i A^j(q) A^k(-q)\Tr\bigg(g\p_j g^{-1}\p_i g\p_{k}g^{-1} g\bigg).
\eeq
Now we need to massage the expression a little bit. Let us use $(ijk)$ as a short-hand notation of the expression $\Tr\bigg(g\p_i g^{-1}\p_j g\p_{k}g^{-1} g\bigg)$. Integrating by parts\cite{Tim2011c}, we have the following relation:
\beq
(ijk)+(jki)+(kij)=0.
\eeq
We therefore have
\beq
(jik)=\frac23\big(2(jik)-(ikj)-(kji)\big).
\eeq
In the trace in Eq.~(\ref{sme}), only the part symmetric under the exchange of the index $j$ and $k$ would contribute, as we can change variables from $q$ to $-q$, effectively exchanging $A^j(q)$ and $A^k(-q)$. Therefore, in the expression above, we can exchange $j$ and $k$ freely. We therefore have
\bea
S_{ME}&\sim&\int\frac{{\rm d^4}q}{(2\pi)^4}\frac i{6} \omega q^i A^j(q) A^k(-q)
\big((jik)+(kij)-(ijk)-(kji)\big)\nonumber\\
&=&\int\frac{{\rm d^4}q}{(2\pi)^4}\frac i{6}\omega q^i A^j(q) A^k(-q)\ep_{ij\ell}\ep_{ab\ell}\big((kab)+(bak)\big).\nonumber\\
&=&\int\frac{{\rm d^4}q}{(2\pi)^4}\frac i{6} B^\ell(q) E^k(-q)\ep_{ab\ell}\big((kab)+(bak)\big)\nonumber\\
&\equiv&-\int\frac{{\rm d^4}q}{(2\pi)^4} B^\ell(q) E^k(-q)\alpha_{k\ell}(q\rightarrow 0).
\eea

$\alpha_{k\ell}(q\rightarrow 0)=-\frac i{6}\ep_{ab\ell}\big((kab)+(bak)\big)$ is traceless, as the two terms cancel each other with antisymmetrization. 
\end{widetext}
\end{document}
%
% ****** End of file template.aps ******


magnetization.bbl
%merlin.mbs apsrev4-1.bst 2010-07-25 4.21a (PWD, AO, DPC) hacked
%Control: key (0)
%Control: author (8) initials jnrlst
%Control: editor formatted (1) identically to author
%Control: production of article title (-1) disabled
%Control: page (0) single
%Control: year (1) truncated
%Control: production of eprint (0) enabled
\begin{thebibliography}{16}%
\makeatletter
\providecommand \@ifxundefined [1]{%
 \@ifx{#1\undefined}
}%
\providecommand \@ifnum [1]{%
 \ifnum #1\expandafter \@firstoftwo
 \else \expandafter \@secondoftwo
 \fi
}%
\providecommand \@ifx [1]{%
 \ifx #1\expandafter \@firstoftwo
 \else \expandafter \@secondoftwo
 \fi
}%
\providecommand \natexlab [1]{#1}%
\providecommand \enquote  [1]{``#1''}%
\providecommand \bibnamefont  [1]{#1}%
\providecommand \bibfnamefont [1]{#1}%
\providecommand \citenamefont [1]{#1}%
\providecommand \href@noop [0]{\@secondoftwo}%
\providecommand \href [0]{\begingroup \@sanitize@url \@href}%
\providecommand \@href[1]{\@@startlink{#1}\@@href}%
\providecommand \@@href[1]{\endgroup#1\@@endlink}%
\providecommand \@sanitize@url [0]{\catcode `\\12\catcode `\$12\catcode
  `\&12\catcode `\#12\catcode `\^12\catcode `\_12\catcode `\%12\relax}%
\providecommand \@@startlink[1]{}%
\providecommand \@@endlink[0]{}%
\providecommand \url  [0]{\begingroup\@sanitize@url \@url }%
\providecommand \@url [1]{\endgroup\@href {#1}{\urlprefix }}%
\providecommand \urlprefix  [0]{URL }%
\providecommand \Eprint [0]{\href }%
\providecommand \doibase [0]{http://dx.doi.org/}%
\providecommand \selectlanguage [0]{\@gobble}%
\providecommand \bibinfo  [0]{\@secondoftwo}%
\providecommand \bibfield  [0]{\@secondoftwo}%
\providecommand \translation [1]{[#1]}%
\providecommand \BibitemOpen [0]{}%
\providecommand \bibitemStop [0]{}%
\providecommand \bibitemNoStop [0]{.\EOS\space}%
\providecommand \EOS [0]{\spacefactor3000\relax}%
\providecommand \BibitemShut  [1]{\csname bibitem#1\endcsname}%
\let\auto@bib@innerbib\@empty
%</preamble>
\bibitem [{\citenamefont {Essin}\ \emph {et~al.}(2009)\citenamefont {Essin},
  \citenamefont {Moore},\ and\ \citenamefont {Vanderbilt}}]{Essin2009}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.~M.}\ \bibnamefont
  {Essin}}, \bibinfo {author} {\bibfnamefont {J.~E.}\ \bibnamefont {Moore}}, \
  and\ \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Vanderbilt}},\
  }\href {\doibase 10.1103/PhysRevLett.102.146805} {\bibfield  {journal}
  {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf {\bibinfo {volume} {102}},\
  \bibinfo {pages} {146805} (\bibinfo {year} {2009})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Essin}\ \emph {et~al.}(2010)\citenamefont {Essin},
  \citenamefont {Turner}, \citenamefont {Moore},\ and\ \citenamefont
  {Vanderbilt}}]{Essin2010}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.~M.}\ \bibnamefont
  {Essin}}, \bibinfo {author} {\bibfnamefont {A.~M.}\ \bibnamefont {Turner}},
  \bibinfo {author} {\bibfnamefont {J.~E.}\ \bibnamefont {Moore}}, \ and\
  \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Vanderbilt}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Phys. Rev. B}\
  }\textbf {\bibinfo {volume} {81:205104,2010}} (\bibinfo {year} {2010})},\
  \Eprint {http://arxiv.org/abs/1002.0290} {1002.0290} \BibitemShut {NoStop}%
\bibitem [{\citenamefont {Malashevich}\ \emph {et~al.}(2010)\citenamefont
  {Malashevich}, \citenamefont {Souza}, \citenamefont {Coh},\ and\
  \citenamefont {Vanderbilt}}]{Malashevich2010}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Malashevich}}, \bibinfo {author} {\bibfnamefont {I.}~\bibnamefont {Souza}},
  \bibinfo {author} {\bibfnamefont {S.}~\bibnamefont {Coh}}, \ and\ \bibinfo
  {author} {\bibfnamefont {D.}~\bibnamefont {Vanderbilt}},\ }\href@noop {}
  {\bibfield  {journal} {\bibinfo  {journal} {New J.Phys.}\ }\textbf {\bibinfo
  {volume} {12:053032,2010}} (\bibinfo {year} {2010})},\ \Eprint
  {http://arxiv.org/abs/1002.0300} {1002.0300} \BibitemShut {NoStop}%
\bibitem [{\citenamefont {Coh}\ \emph {et~al.}(2011)\citenamefont {Coh},
  \citenamefont {Vanderbilt}, \citenamefont {Malashevich},\ and\ \citenamefont
  {Souza}}]{Coh2011}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Coh}}, \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Vanderbilt}},
  \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont {Malashevich}}, \ and\
  \bibinfo {author} {\bibfnamefont {I.}~\bibnamefont {Souza}},\ }\href
  {\doibase 10.1103/PhysRevB.83.085108} {\bibfield  {journal} {\bibinfo
  {journal} {Phys. Rev. B}\ }\textbf {\bibinfo {volume} {83}},\ \bibinfo
  {pages} {085108} (\bibinfo {year} {2011})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Qi}\ \emph {et~al.}(2008)\citenamefont {Qi},
  \citenamefont {Hughes},\ and\ \citenamefont {Zhang}}]{Qi08}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {X.-L.}\ \bibnamefont
  {Qi}}, \bibinfo {author} {\bibfnamefont {T.~L.}\ \bibnamefont {Hughes}}, \
  and\ \bibinfo {author} {\bibfnamefont {S.-C.}\ \bibnamefont {Zhang}},\ }\href
  {\doibase 10.1103/PhysRevB.78.195424} {\bibfield  {journal} {\bibinfo
  {journal} {Phys. Rev. B}\ }\textbf {\bibinfo {volume} {78}},\ \bibinfo
  {pages} {195424} (\bibinfo {year} {2008})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Hasan}\ and\ \citenamefont {Kane}(2010)}]{Hasan2010}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.~Z.}\ \bibnamefont
  {Hasan}}\ and\ \bibinfo {author} {\bibfnamefont {C.~L.}\ \bibnamefont
  {Kane}},\ }\href {\doibase 10.1103/RevModPhys.82.3045} {\bibfield  {journal}
  {\bibinfo  {journal} {Rev. Mod. Phys.}\ }\textbf {\bibinfo {volume} {82}},\
  \bibinfo {pages} {3045} (\bibinfo {year} {2010})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {King-Smith}\ and\ \citenamefont
  {Vanderbilt}(1993)}]{King93}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {R.~D.}\ \bibnamefont
  {King-Smith}}\ and\ \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Vanderbilt}},\ }\href {\doibase 10.1103/PhysRevB.47.1651} {\bibfield
  {journal} {\bibinfo  {journal} {Phys. Rev. B}\ }\textbf {\bibinfo {volume}
  {47}},\ \bibinfo {pages} {1651} (\bibinfo {year} {1993})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Zak}(1989)}]{Zak1989}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Zak}},\ }\href {\doibase 10.1103/PhysRevLett.62.2747} {\bibfield  {journal}
  {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf {\bibinfo {volume} {62}},\
  \bibinfo {pages} {2747} (\bibinfo {year} {1989})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Thonhauser}\ \emph {et~al.}(2005)\citenamefont
  {Thonhauser}, \citenamefont {Ceresoli}, \citenamefont {Vanderbilt},\ and\
  \citenamefont {Resta}}]{Thonhauser05}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.}~\bibnamefont
  {Thonhauser}}, \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Ceresoli}}, \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Vanderbilt}}, \ and\ \bibinfo {author} {\bibfnamefont {R.}~\bibnamefont
  {Resta}},\ }\href {\doibase 10.1103/PhysRevLett.95.137205} {\bibfield
  {journal} {\bibinfo  {journal} {Phys. Rev. Lett.}\ }\textbf {\bibinfo
  {volume} {95}},\ \bibinfo {pages} {137205} (\bibinfo {year}
  {2005})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Ceresoli}\ \emph {et~al.}(2006)\citenamefont
  {Ceresoli}, \citenamefont {Thonhauser}, \citenamefont {Vanderbilt},\ and\
  \citenamefont {Resta}}]{Ceresoli06}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Ceresoli}}, \bibinfo {author} {\bibfnamefont {T.}~\bibnamefont
  {Thonhauser}}, \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Vanderbilt}}, \ and\ \bibinfo {author} {\bibfnamefont {R.}~\bibnamefont
  {Resta}},\ }\href {\doibase 10.1103/PhysRevB.74.024408} {\bibfield  {journal}
  {\bibinfo  {journal} {Phys. Rev. B}\ }\textbf {\bibinfo {volume} {74}},\
  \bibinfo {pages} {024408} (\bibinfo {year} {2006})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Thonhauser}\ and\ \citenamefont
  {Vanderbilt}(2006)}]{Thonhauser06}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.}~\bibnamefont
  {Thonhauser}}\ and\ \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Vanderbilt}},\ }\href {\doibase 10.1103/PhysRevB.74.235111} {\bibfield
  {journal} {\bibinfo  {journal} {Phys. Rev. B}\ }\textbf {\bibinfo {volume}
  {74}},\ \bibinfo {pages} {235111} (\bibinfo {year} {2006})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Thouless}(1984)}]{Thouless1984}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~J.}\ \bibnamefont
  {Thouless}},\ }\href {http://stacks.iop.org/0022-3719/17/i=12/a=003}
  {\bibfield  {journal} {\bibinfo  {journal} {Journal of Physics C: Solid State
  Physics}\ }\textbf {\bibinfo {volume} {17}},\ \bibinfo {pages} {L325}
  (\bibinfo {year} {1984})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Soluyanov}\ and\ \citenamefont
  {Vanderbilt}(2011)}]{Soluyanov2011}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.~A.}\ \bibnamefont
  {Soluyanov}}\ and\ \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Vanderbilt}},\ }\href {\doibase 10.1103/PhysRevB.83.035108} {\bibfield
  {journal} {\bibinfo  {journal} {Phys. Rev. B}\ }\textbf {\bibinfo {volume}
  {83}},\ \bibinfo {pages} {035108} (\bibinfo {year} {2011})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Chen}\ and\ \citenamefont
  {Lee}(2011{\natexlab{a}})}]{Tim2011c}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {K.-T.}\ \bibnamefont
  {Chen}}\ and\ \bibinfo {author} {\bibfnamefont {P.~A.}\ \bibnamefont {Lee}},\
  }\href {\doibase 10.1103/PhysRevB.84.205137} {\bibfield  {journal} {\bibinfo
  {journal} {Phys. Rev. B}\ }\textbf {\bibinfo {volume} {84}},\ \bibinfo
  {pages} {205137} (\bibinfo {year} {2011}{\natexlab{a}})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Chen}\ and\ \citenamefont
  {Lee}(2011{\natexlab{b}})}]{Tim2011a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {K.-T.}\ \bibnamefont
  {Chen}}\ and\ \bibinfo {author} {\bibfnamefont {P.~A.}\ \bibnamefont {Lee}},\
  }\href {\doibase 10.1103/PhysRevB.83.125119} {\bibfield  {journal} {\bibinfo
  {journal} {Phys. Rev. B}\ }\textbf {\bibinfo {volume} {83}},\ \bibinfo
  {pages} {125119} (\bibinfo {year} {2011}{\natexlab{b}})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Swingle}\ \emph {et~al.}(2011)\citenamefont
  {Swingle}, \citenamefont {Barkeshli}, \citenamefont {McGreevy},\ and\
  \citenamefont {Senthil}}]{Brian2011}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {B.}~\bibnamefont
  {Swingle}}, \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Barkeshli}},
  \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont {McGreevy}}, \ and\
  \bibinfo {author} {\bibfnamefont {T.}~\bibnamefont {Senthil}},\ }\href
  {\doibase 10.1103/PhysRevB.83.195139} {\bibfield  {journal} {\bibinfo
  {journal} {Phys. Rev. B}\ }\textbf {\bibinfo {volume} {83}},\ \bibinfo
  {pages} {195139} (\bibinfo {year} {2011})}\BibitemShut {NoStop}%
\end{thebibliography}%
