\section{Introduction}
\label{sec:intro}

More than a century ago, in 1879 Edwin Hall discovered the Hall effect \cite{hall1879}, which occurs in a sheet of metal in a perpendicular magnetic field. An electric voltage perpendicular to the charge current is induced by the magnetic field, due to the Lorentz force acting on electrons, as is illustrated in Fig. \ref{Hall}(a).
In a two dimensional electron gas (2DEG) under low temperature and strong magnetic field, the Hall conductance becomes quantized in the unit of $\frac{e^2}{h}$, while longitudinal resistance drops to zero. This is known as the integer quantum Hall (IQH) effect \cite{klitzing1980}. The IQH effect occurs when electrons enter an incompressible state, with a robust uni-directional, or ``chiral", edge state propagating along the boundary, and all other states localized, as shown in Fig. \ref{Hall}(d). At a lower temperature, fractional quantum Hall (FQH) effect \cite{tsui1982} can occur, which are realized by strongly correlated electron states with fractionally quantized Hall conductivity (in the same unit of $\frac{e^2}{h}$).

One year after his discovery of the Hall effect, Edwin Hall also noticed that the Hall voltage was ten times larger in ferromagnetic materials than that in non-magnetic conductors, which is later known as the "anomalous Hall effect" \cite{hall1881} (Fig. \ref{Hall}(b)). The physical origin of anomalous Hall effect has been studied extensively in the past 50 years \cite{nagaosa2010}. Although the general mechanism is complicated with multiple contributions, the key ingredients responsible for the anomalous Hall effect that have been studied are magnetism, spin-orbit coupling and disorder effects. With the discovery of the quantum Hall effect, it is natural to ask whether there is a similar quantum version of the anomalous Hall effect, or the quantum anomalous Hall (QAH) effect. More precisely, the question is whether a quantized Hall conductance can be induced by mechanisms different from an external magnetic field, such as magnetism, spin-orbit coupling and/or disorder. An extreme form of this question is whether a state with a quantized Hall conductance can exist without external magnetic field at all. The first step of answering this question came from understanding the general topological reason of the Hall conductance quantization \cite{thouless1982,avron1985,niu1985}. A geometrical gauge field can be defined in the space of momentum or twist-boundary conditions, which has a quantized total flux that determines the Hall conductance quantum. Therefore, the IQH effect can exist in more general systems as long as the geometrical gauge field is nontrivial. In 1988, F. D. M. Haldane proposed the first explicit model realizing quantized Hall conductance without orbital magnetic field \cite{haldane1988}, which is known as the Haldane model. The Haldane model is a tight-binding model defined on a honeycomb lattice with real nearest neighbor hopping and complex next-nearest-neighbor hopping.

Although the Haldane model demonstrated the theoretical possibility of realizing a QAH state, how to realize such a state in realistic materials remained unknown for another two decades. The solution was eventually discovered through a ``detour" into another cousin of Hall effect---the spin Hall effect and the quantum spin Hall (QSH) effect. In a system with spin-orbit coupling, electrons with different spin directions can experience different forces. In particular, the motion of electrons in an external electric field can be biased to different directions depending on their spin direction, as is illustrated in Fig. \ref{Hall}(c). This is named as the spin Hall effect \cite{sinova2004,murakami2003,hirsch1999}, since in the direction perpendicular to the electric field there is a net spin current though no net charge current. Different from the charge Hall effect, the spin Hall effect does not require time-reversal symmetry breaking, and thus can occur in a semiconductor without magnetic field. Since 2005, the quantum version of the spin Hall effect (Fig. \ref{Hall}(f)) was proposed theoretically \cite{kane2005a,bernevig2006a,bernevig2006c} and realized experimentally in HgTe/CdTe quantum wells \cite{koenig2007} soon after the theoretical proposal \cite{bernevig2006c}. The discovery of the QSH effect has lead to the latter discovery of topological insulators and topological superconductors \cite{hasan2010,qi2011}, which greatly expanded our understanding of topological states of matter. The QSH effect is deeply related to the QAH effect, since the simplest models realizing QSH effect describe electrons with opposite spin carrying opposite quantized Hall conductance, as is illustrated in Fig. \ref{Hall}(f). The QH state of each spin component can be either a IQH state with Landau levels \cite{bernevig2006a} or a QAH state \cite{kane2005a,qi2005}. Although the QSH state is only protected by time-reversal symmetry and does not require the two spin components to be decoupled \cite{kane2005b}, the simple QSH states with opposite QAH of each spin component pointed out a new route of realizing the QAH effect. Intuitively one can say that for electrons in a QSH state, each spin component is already in a QAH state, and the only problem is the Hall conductance of two spin components cancels each other due to time-reversal symmetry. Therefore the only remaining task is to remove the QAH effect of one spin component, so that the net Hall conductance becomes nonvanishing. This goal was achieved in Ref. \cite{liu2009C}, in which HgTe/CdTe quantum wells doped with Mn impurities, in proper thickness range, was proposed to realize the QAH state, given that the Mn spins are polarized. However, Mn spins are paramagnetic instead of ferromagnetic in this material, so that a small external magnetic field is still required to polarize them in order to realize the QAH phase. Indeed, the early transport experiment \cite{buhmann2002} has observed a Hall conductance $-\frac{e^2}h$ plateau starting at a quite small magnetic field and persisting for an extremely wide range of magnetic field, agreeing with the theoretical prediction.

More recently, a new QAH state was proposed in magnetically doped topological insulator films in the family of ${\rm Bi_2Se_3}$\cite{yu2010}, which is ferromagnetic and thus finally leads to the realization of the QAH effect completely without external magnetic field. ${\rm Bi_2Se_3}$ and related materials ${\rm Bi_2Te_3}$ and ${\rm Sb_2Te_3}$ are three-dimensional topological insulators \cite{hjzhang2009,xia2009}. Similar to the QSH state, the three-dimensional TI turns out to be also deeply related to the QAH state. A universal feature of three-dimensional TI is that its two-dimensional surface state always obtains a half quantized Hall conductance $\sigma_H=\frac{e^2}{h}\left(n+\frac12\right)$ (with $n$ an integer) once time reversal is broken \cite{qi2008b,fu2007b}. This becomes the key ingredient of realizing QAH state in the thin film of TI doped with magnetic impurities such as Cr or V. The first experimental realization of this proposal was achieved in Cr doped ${\rm Bi(Sb)_2Te_3}$ thin films \cite{chang2013a}. The QAH state occurs when the bulk carrier density is tuned to be closed to zero by applying a gate voltage. The magnetic field dependence of the Hall conductivity shows a hysteresis loop that agrees with the anomalous nature of the QH state. In particular, the Hall conductance is quantized to $\pm \frac{e^2}h$ at zero external magnetic field.  Since then this QAH state has been firmly verified in more and more works \cite{checkelsky2014,kou2014,bestwick2014,kou2015,feng2015,chang2014,kandala2015}, in both Cr doped and V doped ${\rm (Bi,Sb)_2Te_3}$ films. An overview on the QAH effect is provided in Ref. \cite{wang2014c}.

\begin{figure} [h]
\includegraphics[width=8cm]{Hall.eps}
\caption{(Color online) Six members in the family of Hall effect. (a) Hall effect; (b) Anomalous Hall effect; (c) Spin Hall effect; (d) Quantum Hall effect; (e) Quantum anomalous Hall effect; and (f) Quantum spin Hall effect.  } \label{Hall}
\end{figure}

This review article is organized as follows. In Sec. II, we will first present a theoretical toy model, and then discuss the general conditions for the occurence of the QAH effect. Sec. III is devoted to the theoretical prediction of different classes of realistic QAH materials and the recent experimental realization. Sec. IV contains the conclusion and further discussions.

\section{Realistic materials and experiments}
\label{sec:mat}

\subsection{Magnetically doped HgTe and InAs/GaSb quantum wells}
HgTe and CdTe are typical II-VI group compound semiconductors, which possess the zinc-blende lattice structure with both the anion and cation atoms forming two interpenetrating face-centered-cubic lattices. Similar to other II-VI and III-V semiconductors, the bands of HgTe and CdTe near Fermi energy are a s-type band ($\Gamma_6$), and a p-type band split by spin-orbit coupling into a $J=3/2$ band ($\Gamma_8$) and a $J=1/2$ band ($\Gamma_7$), where $J$ denotes the total angular momentum. CdTe has an energy gap ($\sim 1.6$eV) and its band ordering is normal with s-type ($\Gamma_6$) bands above p-type ($\Gamma_7$ and $\Gamma_8$) bands. In constrast, HgTe has an inverted band structure \cite{bernevig2006c}: $\Gamma_6$ bands are below $\Gamma_8$ bands with a negative energy gap of $\sim 0.3$eV. HgTe can form a quantum well (QW) structure with CdTe as a barrier, in which both the conduction and valence bands are splitted into sub-bands, denoted as $|E_n\rangle$ for electron $\Gamma_6$ bands and $|H_n\rangle$ for heavy-hole $\Gamma_8$ bands ($n=1,2,\cdots$), due to the quantum confinement. Because of the opposite band sequence between CdTe and HgTe, a topological phase transition exists in this system at a critical thickness of the quantum well, denoted as $d_c$ ($d_c\approx 6.3$nm), across which the band sequence between $|E_1\rangle$ and $|H_1\rangle$ sub-bands is reversed. For the thickness $d<d_c$, $|E_1\rangle$ state has higher energy than $|H_1\rangle$ state (normal band structure) while for $d>d_c$, $|E_1\rangle$ is below $|H_1\rangle$, leading to a similar inverted band structure as descibed in the last section. However, different from our simple model, we need to consider electron spin states for both $|E_1\rangle$ and $|H_1\rangle$ sub-bands. Thus, under the basis $|E_1,+\rangle$, $|H_1,+\rangle$, $|E_1,-\rangle$ and $|H_1,-\rangle$ (here $\pm$ labels spin states), the Hamiltonian describing the low energy physics in HgTe/CdTe quantum wells is written as
\begin{eqnarray}
      H=\left(
      \begin{array}{cc}
	      h(k)&0\\
	      0&h^*(-k)
      \end{array}
      \right),
	\label{eq:BHZ}
\end{eqnarray}
where $h(k)$ is given by the Hamiltonian (\ref{eq:twoband}). Therefore, this Hamiltonian is just two copies of the two band model with each copy described by the minimal two-band model (\ref{QAH}) and is first described by B. A. Bernevig, T. L. Hughes and S. C. Zhang\cite{bernevig2006c}, known as the BHZ Hamiltonian. These two copies are for spin up and spin down states and related to each other by TR symmetry. The BHZ Hamiltonian forms the basis of the QSH effect, for which edge states also exist and reveal a helical nature.

\begin{figure} [h]
\includegraphics[width=8cm]{HgTe.eps}
\caption{(Color online) (a) Evolution of band structure and edge states upon increasing the spin splitting for Mn doped HgTe quantum wells. Band inversion occurs for one species of spin but not for the other. This figure and more details can be found in Ref. \cite{liu2009C}. The Landau levels are shown as a function of magnetic fields for the case (b) with initial magnetization and (c) without initial magnetization.  } \label{HgTe}
\end{figure}

The band inversion for electrons with both spins does not directly give QAH effect since time-reversal symmetry has not been broken. However, it does make the realization of QAH state easier, since there is actually two QAH effects (for opposite spins) canceling each other, and we simply need to break time reversal and turn off one of them. For this purpose, magnetization is introduced to this system\cite{liu2009C}, so that the two spin blocks are no longer related. Physically, magnetization is introduced by Mn doping. ${\rm Mn}^{2+}$ substitutes ${\rm Hg^{2+}}$, and the localized spin of Mn ion couples with that of itinerant electrons through exchange coupling
\begin{eqnarray}
	H_s=\left(
	\begin{array}{cccc}
		G_E&0&0&0\\
		0&G_H&0&0\\
		0&0&-G_E&0\\
		0&0&0&-G_H
	\end{array}
	\right),
	\label{eq:Hex}
\end{eqnarray}
where spin splitting is $2G_E$ for the $|E_1,\pm\rangle$ band and $2G_H$ for the $|H_1,\pm\rangle$ band. The corresponding energy gap is $E_g^+=2M_0+G_E-G_H$ for spin up block while $E_g^-=2M_0-G_E+G_H$ for spin down block. To realized the QAH effect, we require that (1) the state with one kind of spin is in the inverted regime while the other goes to the normal regime; and (2) the entire system is in the insulating regime. Combining these two requirement leads to a simple condition $G_EG_H<0$, which means that the spin splittings for electron sub-bands and hole sub-bands are opposite. As shown in Fig. \ref{HgTe}(a), in this case, the energy gap for one copy of spin is increased while that for the other spin is decreased. Thus, the low energy physics is just described by the two-band model (\ref{eq:twoband}), which can result in the QAH phase.
In the other regime where $G_EG_H>0$, electrons with opposite spin still acquire different gaps, but the system becomes metallic before the two blocks develop an opposite Hall conductance.

The above discussion suggests that to realize inverted band structure for one spin but normal band dispersion for the other, the key condition is that spin splitting for electron and heavy hole sub-bands has opposite sign ($G_EG_H<0$). Fortunately, in HgTe QW doped with Mn, sp-d exchange coupling indeed gives the opposite signs for spin splitting of $|E_1\rangle$ and $|H_1\rangle$ bands. The underyling reason is different exchange mechanism for the coupling between electron spin and magnetic moments for $\Gamma_6$ and $\Gamma_8$ bands. For $\Gamma_8$ bands consisting of p orbitals, there is a large hybridization between p orbital from anions, such as Te, and d orbital from Mn which replace cation positions. As a result, the so-called p-d exchange mechanism gives rise to antiforromagnetic coupling between $\Gamma_8$ bands and magnetic moments. In contrast, the $\Gamma_6$ bands come from s orbitals of cations, such as Hg, and its hybridization to d orbitals of Mn atoms is forbidden by symmetry. Thus, the direct exchange coupling is dominant here, leading to ferromagnetic coupling. The direct calculation based on the eight band Kane model indeed confirms the opposite exchange coupling for electron and heavy hole sub-bands, suggesting that the required inverted band structure can be satisfied in this system.

However, the difficulty of this system lies in the fact that it is a paramagnetic rather than ferromagnetic material. The magnetization in Mn doped HgTe satisfies the Curie-Weiss rule
\begin{eqnarray}
	\langle\mathcal{S}\rangle=-S_0B_{5/2}\frac{5g_{Mn}\mu_BB}{2k_B(T+T_0)},
	\label{eq:CurieWeiss}
\end{eqnarray}
where $S_0=5/2$, $B_{5/2}$ is the Brillouin function and $T_0$ stands for a weak antiferromagnetic coupling between magnetic moments. Due to $T_0$, Mn doped HgTe reveals spin glass behaviors at low temperatures, which is the main obstacle to observe the QAH effect in this system experimentally. Therefore, a small magnetic field is required to magnetize the system and the orbital effect (Landau levels) is inevitable. However, it is still possible to distinguish the effect of magnetization from the orbital effect by checking the behaviors of Landau levels. We may first check the Landau level with initial magnetization at a zero magnetic field, as shown in Fig. \ref{HgTe}(b). One can see that the gap at a zero magnetic field corresponds to the regimes with $-e^2/h$ Hall conductance in the Landau level fan diagram. This indicates that the QAH effect and the corresponding QH effect are adiabatically connected. For the realistic situation, there is no magnetization at a zero magnetic field, so the zero field gap is connected to the zero Hall conductance regime. However, from Fig. \ref{HgTe}(c), one can see that the energy of the zeroth Landau level of the valence band (black line) first shows a rapid and nonlinear increase in energy at low magnetic field and then follow the linear behavior at higher field. This is because magnetic moments can be easily polarized by magnetic field and saturate at a magnetic field weaker than 1T. Therefore, if Fermi energy is tuned into the bulk gap, we expect the quantized Hall conductance can be observed at a small magnetic field. Indeed, transport experiments performed by Molenkamp's group \cite{buhmann2002} have indeed observed an extremely long quantized Hall plateau with the Hall conductivity $-e^2/h$ with magnetic fields ranging from 0.3T up to 30T. This provides the evidence of the quantized Hall conductance plateau induced by magnetization, rather than the magnetic field in this system. Moreover, different behaviors between exchange coupling of magnetic moments and the orbital effect from magnetic fields can also be unveiled by rotating magnetic fields \cite{hsu2013}.

More recently, it was proposed that the QAH effect can be realized in another diluted magnetic semiconductor heterostructure, Mn doped type II InAs/GaSb quantum wells \cite{wang2014}. Both InAs and GaSb are III-V group semiconductors. The unique feature of these two materials is that the conduction band minimum of InAs has lower energy than the valence band maximum of GaSb, due to the large band offset. As a consequence, for the InAs/GaSb/AlSb quantum wells, in which InAs and GaSb together serve as well layers and AlSb layers are for barriers, the first electron subband $|E_1\rangle$ of InAs layers lies below the first hole subband $|H_1\rangle$ of GaSb layers, when the well thickness is large enough. The effective model for this system is essentially the same as that for HgTe/CdTe QWs and this material has been theoretically predicted to be a QSH system \cite{liu2008}. Therefore, it is natural to expect that magnetic doping can drive this system into the QAH state. It turns out that there are two salient features making this system more attractive than Mn doped HgTe QWs. Firstly, Since the $|E_1\rangle$ band in the InAs layer is lower than the $|H_1\rangle$ band in the GaSb layer, there is a charge transfer between these two layers. Therefore, if one neglects the coupling between these two layers, InAs is electron-doped while GaSb is hole-doped. The coupling between these two layers can open a hybridization gap. However, since two bands are separated into two layers, this hybridization gap is quite small. Correspondingly, in the effective model (\ref{eq:BHZ}), the coefficient $A$ in the linear term is much smaller in this system than that in HgTe QWs. Due to the small $A$ term, energy dispersion reveals a Mexican hat shape for both the conduction and valence bands, resulting in a peak of the density of states around band edges. This band edge singularity is responsible for the enhancement of spin susceptibility discussed below. Secondly, both Mn doped InAs and GaSb are diluted magnetic semiconductors, and the exchange mechanism is consistent with free-hole-mediated ferromagnetism within a mean field Zener model \cite{dietl2000,jungwirth2006,sato2010a,dietl2010}. If one neglects the coupling between InAs and GaSb layers, there are free carriers in both layers due to charge transfer, leading to ferromagnetism in the system. However, when hybridization gap is taken into account, free-hole mediated mechanism does not work. Thanks to the band edge singularity, spin susceptibility is enhanced even for the Fermi energy in the hybridization gap. Thus, the ferromagnetic insulating phase is possible in this system. By combining the first principles calculations and the analytical Kane model, it was found that both conditions can be naturally satisfied for magnetically doped InAs/GaSb quantum wells, which provide good candidates for realizing the QAH insulator at a high temperature.

Mn doped HgTe QWs and InAs/GaSb QWs provide examples of the QAH states in II-VI and III-V compound semiconductor heterostructures. A more general scheme based on p-n junction structure has been proposed by H.J. Zhang {\it et al} \cite{zhang2014} for other II-VI, III-V and IV semiconductors. The authors have shown that for some narrow gap semiconductors, inverted band structure can naturally exist at the interface between p and n-type region. These structures are dubbed ``junction quantum wells'' in the limit where both the p-type and n-type regions are narrow, and comparable in size to the interface region, and can host the QSH state. This scheme can also be applied to conventional diluted magnetic semiconductors, leading to the QAH state.

\subsection{Magnetically doped (BiSb)$_2$(TeSe)$_3$ family of materials}

Tetradymite semiconductors Bi$_2$Te$_3$, Bi$_2$Se$_3$ and Sb$_2$Te$_3$ is a family of three dimensional TR invariant TIs with the helical surface state of a single Dirac cone at each surface\cite{hjzhang2009,xia2009}. It has been shown that a thin film of 3D TI crosses over to a 2D insulator which is either a TI (the QSH state) or a trivial insulator, depending on the thickness of the film in an oscillatory fashion\cite{liu2010a,lu2010}. Therefore, inverted band structures can be realized in this class of materials in a controllable manner. In addition, it is remarkable that the ferromagnetism can be realized in this class of materials by doping with proper transition metal elements, such as Cr, Fe and V, even in the insulating regime, which is predicted by Yu {\it et al} \cite{yu2010}. Based on the first principles calculations and the mean field Zener model, Yu {\it et al} showed that spin susceptibility is significantly enhanced due to the Van Vleck mechanism in the inverted regime. Therefore, inverted band structure can drive the system into a ferromagnetic insulating phase. This  prediction was soon confirmed by the experiment, in which by tuning the chemical composition of Te and Se in Bi$_2$(Te$_{1-x}$Se$_x$)$_3$, a topological phase transition between TIs and normal insulators is accompanied by a magnetic phase transition between ferromagnetism and paramagnetism \cite{zhang2013a}. Therefore, both conditions for the QAH effect are naturally satisfied in this class of materials.

To understand the mechanism of QAH effect in this family of materials, we do not need to start from the full 3D theory. Instead, it is sufficient to start from the two massless surface states at the top and bottom surfaces of a TI film.
For each surface Dirac cone, its effective model is similar to our two band model, but with the mass term vanishing $M({\bf k})=0$ due to time reversal symmetry. When magnetization is introduced in the surface and couples to surface electron spin, the mass $M$ can be non-zero and open a gap for surface Dirac cones. We emphasize that in this case, the quadratic term $B$ in Eq. (\ref{eq:twoband}) is still zero. For this case, we find the $\hat{\bf d}$ vector is perpendicular to the plane for $k=0$ but lies in the plane for $k\rightarrow \infty$. Therefore, it forms half of a skyrmion configuration, known as a ``meron''. For such a "meron" configuration, the $\hat{\bf d}$ vector only covers half of the unit sphere, leading to a winding number $\pm 1/2$ and corresponding to a Hall conductance $\pm \frac{e^2}{2h}$. The sign of the Hall conductance depends on the magnetization direction (or the sign of mass $M$). This half Hall conductance of 2D massive Dirac model has been studied for a long time in high energy physics, known as the ``parity anomaly''\cite{redlich1984,jackiw1984,haldane1988,fu2007b,qi2008b}. We emphasize that this half Hall conductance only exists at one surface of 3D TIs and distinguish topological surface states from any purely 2D states. In a realistic 2D system, this half Hall conductance cannot be observed since there are always even number of 2D Dirac fermions, known as the ``fermion doubling theorem'' \cite{nielsen1981,nielsen1981a}. For our interest, we may consider a TI thin film, for which there are two surfaces, top and bottom surfaces. Thus, the total number of Dirac cones is two in the whole film, as is required by the fermion doubling theorem. We may introduce magnetization to gap both surfaces. When magnetic moments on both surfaces are in parallel, we find both surface states contribute to the Hall conductance with the same sign, resulting in a total Hall conductance of $\pm e^2/h$, as shown in Fig. \ref{TIsurface}(a) \cite{qi2008b}. This exactly corresponds to the QAH effect that we are looking for. When magnetic moments of two surfaces are opposite, the Hall conductance for two surface states cancels each other, so there is no total Hall conductance (Fig. \ref{TIsurface}(b)). This discussion concludes that uniform magnetization in TI thin films can lead to the QAH effect.

\begin{figure} [h]
\includegraphics[width=8cm]{TIsurface.eps}
\caption{(Color online)  Illustration of the Hall conductance in ferromagnet-topological insulator heterostructure. (a) When the magnetizations at the top and bottom surfaces are parallel, there is a quantized Hall conductance accompanied by a chiral edge state in topological insulator films. (b) When the magnetizations are anti-parallel, the Hall conductances from the top and bottom surfaces cancel each other. Here STI is for strong topoogical insulator and FM is for ferromagnetism \cite{qi2008b}.
} \label{TIsurface}
\end{figure}

In a thin film of TI, one cannot neglect the coupling between top and bottom surface states. A low energy effective model, taking into account this type of coupling, was proposed in Ref. \cite{yu2010} as
\begin{eqnarray}
	H_{sf}=\left(
	\begin{array}{cc}
		v_f(-k_x\sigma_y+k_y\sigma_x)&m_k\\
		m_k&v_f(k_x\sigma_y-k_y\sigma_x)
	\end{array}
	\right),
	\label{eq:Hsurf}
\end{eqnarray}
where $\sigma$ is the Pauli matrices for spin and the upper and lower blocks denote two surfaces respectively. Thus, the basis is given by $|t\uparrow\rangle$, $|t\downarrow\rangle$, $|b\uparrow\rangle$ and $|b\downarrow\rangle$ where $t,b$ are for top and bottom surfaces and $\uparrow,\downarrow$ are for spin. $v_f$ is Fermi velocity and $m_k=M+B(k_x^2+k_y^2)$ describes the coupling between top and bottom surfaces. The exchange coupling between magnetic moments and electron spin is given by
\begin{eqnarray}
	H_{ex}=\left(
	\begin{array}{cc}
		G_0\sigma_z&0\\
		0&G_0\sigma_z
	\end{array}
	\right).
	\label{eq:Hex2}
\end{eqnarray}
At the first sight, this Hamiltonian is different from the Hamiltonian (\ref{eq:BHZ}). One can consider the binding and antibinding states of two surface states, $|\pm\sigma\rangle=(|t\sigma\rangle\pm|b\sigma\rangle)/\sqrt{2}$, with $\sigma=\uparrow,\downarrow$. Under these basis, the Hamiltonian is transformed to
\begin{eqnarray}
	H_{sf}+H_{ex}=\left(
	\begin{array}{cc}
		h(k)+G_0\sigma_z&0\\
		0&h^*(k)-G_0\sigma_z
	\end{array}
	\right)
	\label{eq:Hsurf2}
\end{eqnarray}
under the basis $|+\uparrow\rangle$, $|-\downarrow\rangle$, $|+\downarrow\rangle$
and $|-\uparrow\rangle$. Here $h(k)=m_k\sigma_z+v_f(k_y\sigma_x-k_x\sigma_y)$, which is similar to that in the Hamiltonian (\ref{eq:BHZ}). For exchange coupling, compare to the case of Mn doped HgTe QWs, one can notice that the situation here indeed corresponds to the case with $G_EG_H<0$ in the exchange coupling (\ref{eq:Hex}). Thus, the conditions for the inverted band structure and an ferromagnetic insulating state are automatically satisfied and this analysis suggests the existence of the QAH effect in magnetically doped Bi$_2$Te$_3$ family of materials.

\begin{figure} [h]
\includegraphics[width=8cm]{QAHexp.eps}
\caption{(Color online) (a) The Hall resistance and (b) longitudinal resistance are measured as a function of magnetic fields for different gate voltages in Cr doped (Bi,Sb)$_2$Te$_3$ films \cite{chang2013a}. The quantized Hall resistance is found at the gate voltage $-1.5V$. Precise quantization of Hall resistance and negligible longitudinal resistance at a zero magnetic field were observed in both (c) V doped \cite{chang2014} and (d) Cr doped (Bi,Sb)$_2$Te$_3$ films \cite{bestwick2014}.   } \label{QAHexp}
\end{figure}

This theoretical prediction was first confirmed experimentally in the transport measurement of Cr doped (Bi,Sb)$_2$Te$_3$ thin films \cite{chang2013a}. The first experimental attempt was to achieve stable ferromagnetism in Bi or Sb chalcogenides. It turns out that magnetic property of this family of materials dramatically depends on both parent compounds and magnetic dopants. For example, a robust out-of-plane ferromagnetism was observed in Cr, V and Mn doped Bi$_2$Te$_3$ or Sb$_2$Te$_3$ \cite{chang2013b,hor2010b,kou2013b,kou2013a,li2015,dyck2005}, while spin glassy behavior was observed in Fe doped Bi$_2$Te$_2$ \cite{kim2013a}. For Mn doped Bi$_2$Se$_3$, ferromagnetism is only found on the surface and magnetization measurements reveal hysteresis for both the out-of-plane and in-plane magnetizations \cite{bardeleben2013,zhang2012a}. For the purpose of realizing the QAH effect, we require ferromagnetism out of the plane. Therefore, it is natural to consider Cr, V and Mn doped Bi$_2$Te$_3$ or Sb$_2$Te$_3$.

The first experimental observation of the QAH effect was reported in Cr doped (Bi,Sb)$_2$Te$_3$ thin film with five quintuple layers \cite{chang2013a}. As shown in Fig. \ref{QAHexp}(a) and (b), hysterisis loop in Hall and longitudinal resistivity was observed and Hall resistivity was almost $h/e^2$ at a zero magnetic field when the temperature is lowered to $30mK$. At the same time, the longitudinal resistivity drops significantly to $0.098 h/e^2$ ($\sim2.53k\Omega$). Furthermore, the observed plateau at a zero magnetic field is adiabatically connected to the quantum Hall plateau at high magnetic fields. All these observations confirm the existence of the QAH effect in this system. This experimental breakthrough has immediately aroused great interest of researchers in this field and was soon reproduced by several other experimental groups \cite{checkelsky2014,kou2014,bestwick2014,kou2015,feng2015,chang2014,kandala2015}.

\begin{figure} [h]
\includegraphics[width=8cm]{phase.eps}
\caption{(Color online) The global phase diagram of conductance tensor is constructed by (a) varying gate voltage for different magnetic fields \cite{checkelsky2014} and (b) tilting magnetic fields \cite{kou2015}.  } \label{phase}
\end{figure}

Nevertheless, for the first experiment of the QAH effect, a relatively large residual longitudinal resistance (around $2.5k\Omega$) always remains and the corresponding Hall conductance is around $0.987 e^2/h$ at zero magnetic field. Furthermore, this effect only occurs around 30mK and the longitudinal resistance will increase significantly once the temperature is increased to 1K. Therefore, an improvement of sample quality is required to demonstrate the precise quantization of this effect. In Cr doped (Bi,Sb)$_2$Te$_3$ films, the quantization of the Hall conductance was improved to the value within $0.01\%$ of $e^2/h$ and the longitudinal resistance reached as low as $15\Omega$ through adiabatic demagnetization cooling of the syste
\cite{bestwick2014}. More recently, the QAH effect was observed in V doped (Bi,Sb)$_2$Te$_3$ \cite{chang2014} with the zero-field longitudinal resistance down to $0.00013\pm0.00007 h/e^2$ ($\sim3.35\pm1.76\Omega$), Hall conductance reaching $0.9998\pm0.0006e^2/h$ and the Hall angle becoming as high as $89.993\pm0.004^o$ at $T=25mK$. The advantage of this material lies in the fact that it is a hard ferromagnet with a large coercive field ($H_c>1.0T$) and a relative high Curie temperature.
These observations lay a firm foundation for the future experiments and applications of the QAH effect.

More recent experiments start exploring different aspects of the observed QAH effect. Ref. \cite{kou2014} observed the QAH effect in thicker Cr doped (Bi,Sb)$_2$Te$_3$ films up to 10 quintuple layers. Non-local transport measurements\cite{kou2014} have also been performed to reveal the coexistence of chiral edge modes and non-chiral ones, which has been predicted theoretically \cite{wang2013b}. Ref. \cite{checkelsky2014} first explored the global phase diagram of the QAH effect. By mapping the behavior of the conductivity tensor as a function of temperature, magnetic field and chemical potential, they find the delocalization behvior of the QAH effect quantitatively agrees with the theoretical prediction based on renormalization group analysis for the QH state \cite{wang2014a,kivelson1992,khmel1983}. For example, in Fig. \ref{phase}(a), all lines of conductance tensor for different magnetic fields clearly collapse onto the same semi-circle when the contribution of bulk carriers is not significant. The global phase diagram is further clarified by Ref. \cite{kou2015} by tilting the magnetic field. In Fig. \ref{phase}(b), two semi-circles for conductance tensor are revealed with the longitudinal conductance minimum at the $\pm e^2/h$ Hall plateau, as well as the unique zero resistance plateau \cite{kou2015,feng2015}, as predicted by Wang {\it et al} \cite{wang2014a}. These observations establish the equivalence between the QAH effect and QH effect in term of their topological nature.

\subsection{Other classes of materials}

Besides the two classes of materials reviewed above, the QAH effect was also proposed in a large variety of materials, including graphene, silicene, germanene and stanene systems\cite{haldane1988,qiao2010,tse2011,qiao2014,zhang2012b,ezawa2012,pan2014,wu2014,xu2013}, transition metal oxides \cite{xiao2011,wang2014b,cook2014}, and others \cite{dong2014}. Hybrid structures combining magnetic insulators with TIs are also discussed \cite{garrity2013}.

Graphene of a honeycomb lattice structure was the first model system for the study of the QAH effect, as proposed by Haldane \cite{haldane1988}. In Ref. \cite{raghu2008} it was shown that a strong interaction can induce a spontaneous symmetry breaking and drive the system into the QAH state. More sophisticated theoretical calculations suggest that the QAH can also be realized in graphene once exchange coupling and spin-orbit coupling are introduced \cite{qiao2010,tse2011,qiao2014} by, for example, decoration with 5d transition-metal adatoms \cite{zhang2012b}. The next development along this line is about other monolayer materials with honeycomb lattice structure, such as silecene, germanene and stanene, which are formed by a single layer of buckled honeycomb lattice of silicon, germanium and tin, respectively. Since these elements are heavier than carbon, intrinsic spin-orbit coupling in these materials is much stronger than that in graphene. Thus, once exchange coupling is introduced to these systems, the QAH effect is naturally expected \cite{ezawa2012,pan2014,wu2014,xu2013}. It should be emphasized that for germanene and stanene, low energy physics is dominated by the band inversion between s orbitals and $p_x-p_y$ orbitals, rather than $p_z$ orbitals \cite{xu2013}. It has been shown that ferromagnetism can emerge in germanene and stanene by controlling surface functionality instead of magnetic doping \cite{wu2014}.

Transition metal oxides (TMO) also provide another potential family of materials to search for QAH effect since 4d and 5d elements of transition metals naturally carry large magnetic moments and strong spin-orbit coupling. It was first pointed out by D. Xiao \cite{xiao2011} that the (111) bilayer heterostructure of perovskite TMO can be viewed as a buckled honeycomb lattice and has natively inverted band structure. This leads to material proposals of bilayer LaAuO$_3$, LaNiO$_3$ and LaCoO$_3$ heterostructures \cite{wang2014b} and some TMO materials with double perovskite structure \cite{zhang2014b,cook2014}.

Besides these families of materials, hybrid structures \cite{garrity2013,zhang2014a} have also been discussed for the QAH effect and these systems have the advantage of separating topological nature and ferromagnism into two different materials, which can possibly allow the effect to be realized in a stronger ferromagnet and at higher temperature. In particular, Ref. \cite{zhang2014a} proposed the QAH effect in materials with 5d elements.
Potential candidate materials for QAH effect with the Hall conductance higher than $e^2/h$ have also been considered theoretically\cite{wang2013a,zhang2013b,fang2014,dong2014}. The QAH effect has also been proposed in other magnetic structures \cite{liu2013a,liang2013}.

\section{Conclusion and discussions}

In this article, we have reviewed the theoretical prediction and experimental realization of the QAH effect. By making use of spin-orbit coupling and magnetism, a material without external magnetic field can have quantized Hall conductance, as a consequence of the Berry's phase gauge field in momentum space. This effect has been predicted and experimentally verified in magnetically doped topological insulators. The discovery of the QAH effect has significance in both fundamental science and potential applications. Different from conventional QH effect, the QAH effect is a direct demonstration of {\it intrinsic} topological properties of a material, rather than those imposed by an external magnetic field. Viewing from a different angle, the materials with QAH effect are deeply connected to topological properties of both two-dimensional and three-dimensional topological insulators. The observation of QAH effect has provided a direct demonstration of the topological response properties of topological insulators in electrical transport measurements. In addition, realizing QAH state also opens the door for future experiments to verify other topological properties of TI related to time-reversal breaking at the surface, such as the topological magneto-electric (TME) effect\cite{qi2008b}, the image monopole effect\cite{qi2009} and topological Faraday rotation\cite{qi2008b,maciejko2010b}.

There are various generalizations of QAH effects in other systems and physical situations. Similar to the QH effect, the QAH effect also has a fractional QH generalization in strongly interacting systems, which has been recently proposed and studied theoretically, known as the fractional Chern insulators (FCI) \cite{tang2011,neupert2011a,sun2011,sheng2011,qi2011a,regnault2011,parameswaran2012,barkeshli2012,wu2013,bergholtz2013}.
The realization of FCI phases requires strong electron interaction and a narrow topologically nontrivial band that is partially filled.
QAH effect has also been generalized to systems with periodic time-dependent Hamiltonian, known as the Floquet Chern insulators\cite{lindner2011,kitagawa2011,gu2011,wang2013d}.
Effective topological band structure can be induced by the time-dependent coupling, such as optical transition between different bands for an electron in external photon field. Beyond electron systems, the topological band structure leading to QAH effect has also been proposed and even realized in cold atom systems \cite{wu2008,zhai2014,liu2010b,duca2015,goldman2013,jotzu2014}, photonic crystals \cite{haldane2008,wang2008,wang2009,rechtsman2013} and magnonics\cite{zhang2013c,shindou2013,mook2014}.

QAH effect also has potential applications in future electronic devices. A key signature of the QH and QAH states is the chiral edge state, which can carry electric current without dissipation, since all electrons have the same direction of current, and backscattering is impossible. Dissipationless transport is always very important since Joule heating becomes a more and more significant problem when the size of electronic devices is reduced\cite{kim2003}. If QAH edge states are used to transport electric current, the only source of resistance is the contact resistance which is independent from the length of the edge. Consequently, using such a ``chiral interconnect" will be more efficient than an ordinary metal wire as long as the length exceeds a critical value \cite{zhang2012c}. It has been shown that chiral edge states can also carry spin polarization \cite{wu2014a}, thus enabling the potential applications in spintronics. Another mechanism to realize dissipationless transport is, of course, using superconductors. Currently both QAH state and superconductivity requires very low temperature, which limits their practical application. However, one may suspect that it is in general easier to realize a room-temperature QAH state than to realize a room temperature superconductor, because magnetism can easily occur at room temperature in many materials. More experimental and theoretical efforts in seeking QAH materials are certainly required for finding the route towards room temperature QAH effect.

\section*{ACKNOWLEDGMENTS}
We would like to acknowledge the helpful discussions and collaborations with C. Z. Chang, M. Chan, X. Dai, Z. Fang, C. Felser, T. L. Hughes, X. Liu, L. Molenkamp, N. Samarth, J. Wang, K. L. Wang, Y. Y. Wang, Y. S. Wu, G. Xu, Q. K. Xue, H. J. Zhang and all other collaborators of us in related works. SCZ is supported by the US Department of Energy, Office of Basic Energy Sciences under contract DE-AC02-76SF00515 and the NSF under grant numbers DMR-1305677. XLQ is supported by the National Science Foundation through the grant No. DMR-1151786.

\section{Theory of the quantum anomalous Hall effect}
\label{sec:theory}

\subsection{General setup and history}

\begin{figure} [h]
\includegraphics[width=7cm]{Haldane.eps}
\caption{(Color online) (a) The honeycomb lattice with the complex next nearest neighbor hopping. (b) The phase diagram of Haldane model where the Hall conductance is given by $\sigma_{xy}=\nu e^2/h$. From Ref. \cite{haldane1988}.  } \label{Haldane}
\end{figure}

In the QH effect, the Hall conductance is quantized to the integer value of $e^2/h$. Consider a general non-interacting electron Hamiltonian $H=\sum_{\bf k}\sum_{\alpha,\beta=1}^Nc_{{\bf k}\alpha}^\dagger h^{\alpha\beta}({\bf k})c_{{\bf k}\beta}$, with $\alpha,\beta=1,2,...,N$ labeling $N$ orbitals in each unit cell. Diagonalization of single particle Hamiltonian matrix $h({\bf k})$ gives the eigenstates $\left|n{\bf k}\right\rangle$. The components $\langle\alpha|n{\bf k}\rangle=u_{n{\bf k}}^\alpha$ describe the periodic part of Bloch wavefunctions for the $n$-th band. For a band insulator with $M$ of the $N$ bands fully occupied, the Hall conductance can be calculated directly using Kubo formula, which at zero temperature is reduced to the following formula:
\begin{eqnarray}
	\sigma_{H}&=&\frac{e^2}{h}\frac{2\pi}\int d^2{\bf k}\left(\partial_{k_x}a_y-\partial_{k_y}a_x\right)\label{eq:TKNN}
\end{eqnarray}
where
\begin{eqnarray}
a_{x(y)}&=&-i\sum_{n=1}^M\langle n{\bf k}|\frac{\partial}{\partial k_{x(y)}}|n{\bf k}\rangle.
\end{eqnarray}
${\bf a}=(a_x,a_y)$ is a gauge field in the momentum space which describes the Berry's phase obtained by adiabatically changing momentum ${\bf k}$ in the Brillouin zone. The Hall conductance is the net flux of the Berry's phase gauge field ${\bf a}$ times a universal constant $e^2/2\pi h$. This formula was first studied for the IQH state and is known as the Thouless-Kohmoto-Nightingale-Nijs (TKNN) formula \cite{thouless1982}. Mathematically, the flux of a gauge field on a compact manifold, such as the Brillouin zone torus, is always quantized, which is called the first Chern number.

In the QH effect, non-zero Chern number comes from strong magnetic fields and the associated Landau levels. A natural question arises: can non-zero Chern number exist an electronic band without Landau levels? It has been known \cite{nagaosa2010} for a while that Berry curvature in the momentum space can emerge in an electronic band of a metallic system and leads to the so-called intrinsic anomalous Hall effect. Therefore, there is no constraint to forbid the appearance of non-zero Chern number for a band in an insulating system. The first explicit model with quantized Hall conductance in absent of magnetic fields was constructed by F. D. M. Haldane \cite{haldane1988} in the honeycomb lattice. Besides the nearest neighbor hopping between A and B sublattices in the honeycomb lattice, this tight-binding model takes into account the complex second-neighbour hopping induced by a periodic local magnetic flux. The directions of positive phase hopping is indicated by the arrows in Fig. \ref{Haldane}(a). It should be emphasized that the periodic local magnetic flux is chosen to keep zero total flux through one hexagonal unit cell so that the full lattice translation symmetry is respected. The Hamiltonian of this model is given by
\begin{eqnarray}
	&&H({\bf k})=2t_2\cos\phi (\sum_i\cos({\bf k}\cdot{\bf b}_i)){\bf I}\nonumber\\
	&&+t_1\left( \sum_i[\cos({\bf k}\cdot{\bf a}_i)\sigma_1+\sin({\bf k}\cdot{\bf a}_i)\sigma_2 \right)\nonumber\\
	&&+\left[ M-2t_2\sin\phi \sum_i \sin({\bf k}\cdot{\bf b}_i) \right]\sigma_3,
	\label{eq:Ham_Haldane}
\end{eqnarray}
where $\sigma_i$ are the Pauli matrices for A and B sublattices, $t_1$ is for the nearest neighbour hopping, $t_2$ is for the complex second-neighbour hopping, $\phi$ is the phase accumulated by the $t_2$ hopping and $M$ is for on-site energy between A and B sublattices. ${\bf a}_i$ are the displacements from a B site to its three nearest-neighbor A site and ${\bf b}_i$ are displacements for nearest neighbor sites in the same sublattices. When both $M$ and $t_2$ are zero, two Dirac cones appear at the corner of the Brillouin zone (usually denoted as K and K'), which have been well studied ever since the discovery of graphene. As shown in Fig. \ref{Haldane}(b), when the phase $\phi$ of the next nearest neighbor hopping ($t_2$ term) is included, mass terms appear in the Dirac equation and gap the gapless Dirac cones. This gaped system exhibits a non-zero Hall conductivity ($\sigma_{xy}=\nu e^2/h$), and thus realizes the QAH phase.

\subsection{A minimal two-band model}
\label{sec:minimalmodel}

Although Haldane's remarkable model demonstrates the realization of quantized Hall conductance in the absence of magnetic field, it remains theoretical and has not been realized in real materials with honeycomb lattice. In 2005, Ref. \cite{qi2005} discussed the possibility of the QAH effect in semiconductors, and proposed a different two-band tight-binding model with the QAH phase, which is defined on the square lattice and has low energy excitations around momentum ${\bf k}=0$. Compared to the Haldane model, this model turns out to be more closely related to the realistic material proposals of the QAH states\cite{liu2009C,yu2010}. Furthermore, this model is also closely related to the so-called BHZ model of B. A. Bernevig, T. L. Hughes and S. C. Zhang\cite{bernevig2006c}, which describes the first QSH insulator, the HgTe quantum well. This model can be viewed as a lattice regularized Hamiltonian of Dirac fermions, which can be generalized straightfowardly to lattice Dirac models in general dimensions. It turns out that the lattice Dirac models can be used as ``prototype models" which describes all free-fermion topological insulators in different dimensions and symmetry classes\cite{qi2008b}. In the following, we will introduce this simple two-band model and analyze it's physical properties. Based on the explicit understanding gained in this model, we will then discuss the necessary conditions for realizing the QAH phase in realistic systems.

\begin{figure} [h]
\includegraphics[width=7.5cm]{bandinversion.eps}
\caption{(Color online) (a) Normal band structure; (b) Band structure at the critical point; and (c) Inverted band structure. Here we always keep the parameter $B>0$. } \label{bandinversion}
\end{figure}

The minimal two-band model was first introduced by Qi \cite{qi2005}, given by
\begin{eqnarray}
h({\bf k})=\sum_a d_a({\bf k})\sigma^a.\label{QAH}
\end{eqnarray}
where
\begin{eqnarray}
{\bf d} (\mathbf{k})&=&\left(A \sin k_x,-A \sin k_y,\mathbb{M}(\mathbf{k})\right),\nonumber\\
\mathbb{M}(\mathbf{k})&=& M+2B(2-\cos k_x-\cos k_y).\label{eq:twoband}
\end{eqnarray} \noindent
$A,B,M$ are some material related parameters and their physical meaning depends on detailed properties of materials.

The topological property of the model (\ref{QAH}) is determined by the unit-vector $\hat{\bf d}({\bf k})={\bf d}({\bf k})/|{\bf d}({\bf k})|$ in the momentum space.
For generic models with general form (\ref{QAH}), without specifying the form of vector $d_a({\bf k})$ one can directly simplify the TKNN formula (\ref{eq:TKNN}) of the zero temperature Hall conductivity as
\begin{eqnarray}
\sigma_H=\frac{e^2}{h}\frac1{4\pi}\int dk_x\int dk_y\,\hat{\bf d}\cdot
\left(\frac{\partial
\hat{\bf d}}{\partial k_x}\times\frac{\partial \hat{\bf d}}{\partial
k_y}\right).\label{winding2band}
\end{eqnarray}

The expression (\ref{winding2band}) has an explicit geometric meaning. The function $\hat{\bf d}({\bf k})$ defines a mapping from the momentum ${\bf k}$ to the unit vector ($|\hat{\bf d}|=1$) in a sphere, denoted as $\hat{\bf d}({\bf k}): K \rightarrow S^2$. Here $K$ denote the momentum space and in a lattice, it corresponds to the Brillouin zone, belonging to the torus $T^2$. One notices that $\hat{\bf d}\cdot\left(\frac{\partial\hat{\bf d}}{\partial k_x}\times\frac{\partial \hat{\bf d}}{\partial k_y}\right)$ is exactly the Jacobian of this mapping. Thus, the integration over the whole Brillouin zone gives the area which the vector $\hat{\bf d}$ covers, which equals to the integer number times the area of a unit-sphere surface, $4\pi n$. This integer number $n$ is just the winding number for this mapping\cite{semenoff1984,khveshchenko1989,volovik2003,qi2005}.
Consequently, the Hall conductivity is given by $\sigma_H=n\frac{e^2}{h}$.

A more intuitive picture can be obtained by taking the continuum limit of this model around ${\bf k}=0$, given by
\begin{eqnarray}
d_x=A k_x, d_y=-A k_y, d_z=M+B(k^2_x+k^2_y). \label{eq:twoband2}
\end{eqnarray} \noindent
Here we expand the Hamiltonian up to the second order in ${\bf k}$. Without the parabolic term ($B=0$), this Hamiltonian is nothing but massive 2+1D Dirac Hamiltonian with the Dirac mass $M$, which was studied by Jackiw \cite{jackiw1984}. For our purpose, it turns out that the parabolic term is essential, as discussed below. This continuous model corresponds to the case where the gap for large ${\bf k}$ is much larger than that near ${\bf k}=0$, so the corresponding topological property is determined by low energy physics around ${\bf k}=0$. In this case, the Brillouin zone is approximated by a 2D infinite plane. For most general function ${\bf d}({\bf k})$, the quantization of the integral in Eq. (\ref{winding2band}) may fail because ${\bf k}$ is not defined on a compact manifold anymore. However, this is not a problem in the model in Eq. (\ref{eq:twoband}) because ${\bf d}({\bf k})$ has a simple behavior at large ${\bf k}$. For large enough ${\bf k}$, the parabolic term will dominate over other terms, leading to $\hat{\bf d}=sign(B)\hat{\bf z}$. Therefore, there is no contribution to the winding number (\ref{winding2band}) from large momentum, and the quantization of Hall conductivity still holds.
For the Hamiltonian (\ref{eq:twoband}) together with (\ref{eq:twoband2}), one can explicitly show that $\hat{\bf d}=sign(M)\hat{\bf z}$ for $k=0$ and $\hat{\bf d}=sign(B)\hat{\bf z}$ for $k\rightarrow\infty$. Therefore, when $M/B<0$, the $\hat{\bf d}$ vector is perpendicular to the plane and points in opposite directions at $k=0$ and $k\rightarrow\infty$. For other momenta, the $\hat{\bf d}$ vector tilts to the plane and thus forms a skyrmion configuration in the momentum space with a winding number 1, leading to the Hall conductance $e^2/h$. In contrast, for $M/B>0$, the $\hat{\bf d}$ vector points the same direction for $k=0$ and $k\rightarrow\infty$ and gives rise to zero winding number, as well as zero Hall conductance. Thus, the Hall conductivity is quantized to be $e^2/h$ for the regime $M/B<0$ and zero for $M/B>0$. The energy band dispersions for $M/B>0$ and $M/B<0$ are shown in Fig. \ref{bandinversion}. We called the case with $M/B>0$ as normal band structure and that with $M/B<0$ as the ``inverted'' band structure. The ``inverted'' band structure usually possesses topologically non-trivial property. Between these two distinct phases, a topological phase transition occurs at the gapless critical point with $M=0$ and the corresponding system is described by a massless Dirac fermion. The ``inverted'' band structure provides us an intuitive picture to understand how non-trivial topological property occurs in realistic materials and also a helpful guidance to search for realistic topoloigical materials.

Just like the QH effect, one expects chiral edge modes to appear at the boundary of the above two-band model in the inverted regime as a result of non-zero Hall conductance. To see edge states explicitly, let us consider the eigen equation $h(-i\partial_x,k_y)\psi(x)=E\psi(x)$ of the two-band model in a semi-infinite sample ($x>0$) with the boundary condition $\psi(x=0)=0$ and $\psi(x\rightarrow\infty)=0$. We set $k_y=0$ and search for zero energy state ($E=0$). We take the wave function ansatz $\psi\sim e^{\lambda x}\xi$ with two component spinor $\xi$. With this ansatz, the eigen equation is simplified as $(M-B\lambda^2)\xi+A\lambda \sigma_y\xi=0$. Thus, the spinor $\xi$ should be the eigen state of $\sigma_y$, $\sigma_y\xi_s=s\xi_s$ ($s=\pm$), and $\lambda$ satisfies the equation $B\lambda_s^2-sA\lambda_s-M=0$. This suggests that the general solution of the wave function should be given by
\begin{eqnarray}
	\psi(x)&=&(c_{+1}e^{\lambda_{1}x}+c_{+2}e^{\lambda_{2}x})\xi_{+}\nonumber\\
	&+&(c_{-1}e^{-\lambda_{1}x}+c_{-2}e^{-\lambda_{2}x})\xi_{-},
	\label{eq:sol_edge}
\end{eqnarray}
where $\lambda_{i}$ ($i=1,2$) is given by $\lambda_{1,2}=\frac{1}{2B}(A\pm\sqrt{A^2+2BM})$. The boundary condition $\psi(x=0)=0$ leads to $c_{s1}+c_{s2}=0$ for $s=\pm$ because $\xi_+$ and $\xi_-$ are orthogonal to each other. As a consequence, in order to satisfy the condition $\psi(x\rightarrow\infty)=0$, one needs $\lambda_{1}\lambda_{2}=-M/B>0$, which is exactly the condition for band inversion discussed above. Under this condition, one can further show that for $A/B<0$, $\psi(x)=c_{+1}(e^{\lambda_{1}x}-e^{\lambda_{2}x})\xi_{+}$ due to $Re(\lambda_{1,2})<0$, while $A/B>0$, $\psi(x)=c_{-1}(e^{-\lambda_{1}x}-e^{-\lambda_{2}x})\xi_{-}$. Therefore, only one kind of spin state is normalizable for the zero mode. Furthermore, we can consider the dispersion along the $k_y$ direction for the zero mode. To the leading order, if we assume the wavefunction of the zero mode does not depend on $k_y$, i.e. if the edge state wavefunction with momentum $k_y$ is given by $\phi_{k_y}(x,y)=\psi(x)e^{ik_yy}$, the energy of this state is simply
\begin{eqnarray}
	E\left(k_y\right)=\left\{
	\begin{array}{c}
		-Ak_y \qquad A/B<0\\
		Ak_y \qquad A/B>0
	\end{array}
	\right.
	\label{eq:Heff_edge}
\end{eqnarray}
The linear dispersion shows that the edge state is chiral, with a velocity determined by $A$ and $B$.
The existence of chiral edge modes directly demonstrate the non-zero quanitzed Hall conductivity in the 2D bulk system.

\subsection{Two conditions for realistic materials}
\label{sec:conditions}
In the above, we have presented a minimal model for the QAH effect and obtain the quantized Hall conductivity, as well as chiral edge modes, for this model. Based on this understanding, one may ask what are the physical criteria for the QAH effect that will guide our search for realistic QAH materials. One essential condition is inverted band structure (or $M/B<0$ for the minimal model). In realistic systems, it is known that time reversal invariant topological insulators possess this type of band structure\cite{bernevig2006c}. For topological insulators, inverted band structure is driven by strong spin-orbit coupling and time reversal symmetry is preserved. Thus, Hall conductance must be zero and its edge/surface state is of helical type, consisting of both spin. We may make an analog and require the band inversion also to occur only for one species of spin, accompanied by the breaking of time reversal symmetry. This can be achieved in ferromagnetic materials with a large spin splitting due to the so-called exchange coupling between magnetic moments and electron spin. Thus, to look for the QAH effect, we require ferromagnetic insulators with a strong exchange coupling.

This disucssion suggests two necessary conditions for the QAH effect: (1) inverted band structures and (2) ferromagnetic insulators. For the first condition, we may consider materials with electronic band structure in the inverted regime or close to the transition point when there is no magnetization. Only for these materials, it is possible for the exchange coupling to drive the system into the inverted regime. This indicates that we should consider TR invariant TIs, of which the band structure is already inverted, or other narrow gap (or zero gap) semiconductors. The second condition is not easy to satisfy since ferromagnetism usually coexists with metallic behaviors and ferromagnetic insulators are rare. The reason is that in most ferromagnetic systems, ferromagnetism originates from RKKY mechanism, of which ferromagnetic coupling is usually mediated by free carriers. For example, the QAH effect was first predicted in Mn doped HgTe quantum wells \cite{liu2009C}. However, due to the lack of free carriers, magnetic moments of Mn atoms in this system is coupled antiferromagnetically rather than ferromagnetically. This difficulty prevents a direct observation of the QAH effect in this system. There are several ferromagnetic insulators, such as EuO and GdN. However, inverted band structure is difficult to be realized in these materials due to the large band gap. This difficulty was overcome in magnetically doped Bi or Sb chalcogenides. A key insight \cite{yu2010} is that strong spin-orbit coupling can significantly enhance spin susceptibility, particularly in the inverted band regime, through the Van Vleck mechanism. With this mechanism for ferromagnetic insulators in an inverted material, the last missing piece for realizing QAH state was found, and this theoretical prediction was verified by experiments soon afterwards. In the next section, we will follow the history of this field and discuss Mn doped HgTe quantum wells first, which is also helpful for illustrating the general mechanism of relating the QAH effect with time-reversal invariant topological insulators. Then we will discuss the rapid development of the theoretical and experimental studies of this effect in magnetically doped Bi or Sb chalcogenides. Finally, we will comments on other possible materials.


\section{Introduction}

Nim forms the foundation of the mathematical study of two-player strategy games. In his landmark $1901$ paper, \textit{Nim, a game with a complete mathematical theory}, Charles L. Bouton provided a solution to the game of Nim, essentially founding the field of Combinatorial Game Theory \cite{Nim}.

Since Bouton's discovery, many extensions or variants of Nim have been explored. Some variations that come to mind are Wythoff's Game, Poker Nim and Kayles. These variations often yield winning strategies that bare little resemblance to that of Nim \cite{WinningWays}.

Interestingly, very few if any of these variations use moves predicated upon modular congruence. In this paper, we explore a modular extension to Nim, which we call $m$-Modular Nim, in which moves are indeed predicated upon modular congruence are added to the traditional Nim moves.

We start this paper with preliminaries in Section~\ref{prelim}. In Section~\ref{mModularNim}, we introduce the game of Modular Nim which is similar to Nim but in addition to Nim moves it allows players to remove a positive multiple of $m$ tokens total from the position.

Section~\ref{2HeapOdd} considers $2$ heap Modular Nim for odd modular bases. Starting with an example for $m = 3$, we prove that the number of P-positions is finite and equal to $m$ for odd $m$. In Section~\ref{2HeapEven}, we expand our result to even values of $m$ by observing a self-similar structure in the set of P-positions. In Section~\ref{sec:explicit} we describe the P-positions explicitly.

In Section~\ref{sec:manyheaps} we describe P-positions in $m$-Modular Nim for any number of heaps and odd $m$.

\section{Preliminaries}\label{prelim}

We will be investigating the broad field of Combinatorial Game Theory (CGT). Roughly speaking, CGT concerns the study of winning strategies in two-player perfect information games. Our exploration of this large topic begins with some basic yet essential definitions \cite{WinningWays}.

\begin{definition}An \emph{impartial combinatorial game} is a two-player game where each player has both the same moves available at each and every point in the game and a complete set of information about the game and the potential moves.  \end{definition}

This implies that no randomness such as rolling dice can exist.

\begin{definition}In \emph{normal play}, the first player unable to move is declared the loser. \end{definition}

\begin{definition} We will call a position a \emph{terminal position} if no moves may be made from it. \end{definition}

In general, impartial combinatorial games are analyzed using the notion of P-positions and N-positions. This system of notation allows for games to be solved from the bottom up.

\begin{definition}A \emph{P-position} is a position from which the \emph{previous} player will win given perfect play. The set of P-positions is denoted as $\mathcal{P}$. \end{definition}

We can observe that all terminal positions are P-positions.

\begin{definition}An \emph{N-position} is a position from which the \emph{next} player will win given perfect play. The set of N-positions is denoted as $\mathcal{N}$.\end{definition}

Any position in the game is either a P-position or an N-position. All the moves from any P-position lead to an N-position. On the other hand, from any N-position, there exists some move to a P-position. These observations motivate the following theorem \cite{LessonsInPlay}.

\begin{theorem}Suppose that the positions of a finite impartial game can be partitioned into disjoint sets $\mathcal{A}$ and $\mathcal{B}$ such that:

\begin{enumerate}
\item Every move from a position in $\mathcal{A}$ is to a position in $\mathcal{B}$.
\item Every move from a position in $\mathcal{B}$ has at least one move to a position in $\mathcal{A}$.
\item All terminal positions are elements of $\mathcal{A}$.
\end{enumerate}

Then $\mathcal{A} = \mathcal{P}$ and $\mathcal{B} = \mathcal{N}$.
\end{theorem}

\subsection{Nim}

Nim is the most fundamental impartial combinatorial game. It is one of the earliest \emph{take-away games} and is testament to the complexity that can arise from simple rules \cite{LessonsInPlay}.

\begin{definition}In the game of \emph{Nim}, each position consists of a set of heaps of tokens. In a move, a player must remove a positive number of tokens from a single heap. \end{definition}

In a Nim-like game, denote positions with $h$ heaps as ordered $h$-tuples. To describe the set of P-positions in Nim, we need to define the bitwise XOR operation.

\begin{definition} The \textit{bitwise XOR} of two numbers is calculated by writing both numbers in binary and adding them without carrying over. We will use the $\oplus$  symbol to denote the bitwise XOR operation.
\end{definition}

The set of P-positions in Nim is well understood and summarized by the following theorem \cite{Nim}.

\begin{theorem}[Bouton's Theorem]In Nim, $(a_1, \ldots, a_n) \in \mathcal{P}$ if and only if \[\bigoplus_{i=1}^n a_i = 0.\] \end{theorem}

\section{$m$-Modular Nim}\label{mModularNim}

We will now introduce a natural extension to Nim which we will subsequently discuss in detail. Our game involves loosening the restrictions on Nim moves with conditions based on modular congruence.

\begin{definition}In the game of \emph{$m$-Modular Nim}, each position consists of a set of heaps of tokens, like in Nim. However, we have two types of moves:

\begin{enumerate}[Type I.]
\item Remove a positive number of tokens from a single heap.
\item Remove $km$ tokens total where $k$ is a positive integer.
\end{enumerate}
\end{definition}

In our analysis of $m$-Modular Nim, we introduce some additional notation and positional functions for convenience.

\begin{definition}Let the \emph{heap-sum} of a position $A$ be the total number of tokens. We denote it as $|A|$. \end{definition}

We introduce a partial order on the set of positions to allow ourselves to speak more concisely about important concepts.

\begin{definition}If $A = \left(a_1, a_2, a_3, \ldots, a_k\right)$ and $B = \left(b_1, b_2, b_3, \ldots, b_k\right)$ are positions in $m$-Modular Nim such that $a_i \geq b_i$ for all integers $1 \leq i \leq k$, we say that $A$ \emph{dominates} $B$. Position $A$ \emph{strictly dominates} $B$ if $A$ dominates $B$ and $A$ is not equal to $B$, that is, there exists $i$ such that $a_i > b_i$. We denote domination as $A \succeq B$ and strict domination as $A \succ B$.

Moreover if all members of set $\mathcal{S}$ dominate all members of set $\mathcal{T}$, we say that $\mathcal{S}$ dominates $\mathcal{T}$ or $\mathcal{S} \succeq \mathcal{T}$. Similarly, if all members of $\mathcal{S}$ strictly dominate all members of $\mathcal{T}$, we say that $\mathcal{S}$ strictly dominates $\mathcal{T}$ or $\mathcal{S} \succ \mathcal{T}$. \end{definition}

If a P-position $A$ dominates a P-position $B$, one might expect that then there exists an optimal game when $A$ occurs as a position in the game before $B$. This is true for a $2$-heap game, but is not true for a game with more heaps.

For example, consider a 3 heap game of $4$-Modular Nim. The P-position $(1,2,2)$ dominates the P-position $(0,1,1)$, but the latter position can not be reached from the former position in any optimal play.

\begin{lemma}\label{DominateTypeII}
A Type II move from position $A$ to position $B$ exists, if and only if $|A| \equiv |B| \pmod{m}$ and $A \succ B$. \end{lemma}

\begin{proof}
If a move from position $A$ to position $B$ exists then $A \succ B$. If in addition, this is a Type II move, the total number of tokens is decreased by a multiple of $m$ implying $|A| \equiv |B| \pmod{m}$.

On the other hand, suppose $|A| \equiv |B| \pmod{m}$ and $A \succ B$. Let $A = (a_1,a_2,\ldots,a_n)$ and $B = (b_1,b_2,\ldots,b_n)$. From the $i^\text{th}$ heap in $A$, take away $a_i - b_i$ tokens. Because $A \succ B$, we have that $a_i - b_i \geq 0$ implying that the move is well-defined. Moreover, the total number of tokens removed must be divisible by $m$ as $|A| \equiv |B| \pmod{m}$.
\end{proof}

\section{$2$ Heap $m$-Modular Nim for Odd $m$}\label{2HeapOdd}

Rather than dealing with any number of heaps, we will start with $2$ heaps and odd $m$.

\subsection{An Example: $m=3$.}

\begin{example}[$m = 3$]

The P-positions of $3$-Modular Nim are the ordered pairs:

\begin{center}
$(0,0),$ \\
$(1,1),$ \\
$(2,2).$
\end{center}

Let the specified set be $\mathcal{S}$. We can manually verify that no members of $\mathcal{S}$ are connected by a legal move. Thus, it suffices to show that a move from any position $(a,b) \notin \mathcal{S}$  to an element of $\mathcal{S}$ exists.

\begin{enumerate}
\item Suppose that $\min(a,b) < 3$, then a Type I move must exist.
\item On the other hand, if $\min(a,b) \geq 3$, then $(a,b) \succ \mathcal{S}$. By Lemma \ref{DominateTypeII}, a Type II move must exist because all residue classes modulo $3$ are covered by $\mathcal{S}$.
\end{enumerate}
\end{example}

Figure~\ref{2HeapMod3} displays the P-positions in 3-Modular Nim on a coordinate grid.

\begin{figure}[htbp]
\begin{center}
\begin{tikzpicture}[scale = 1]
{\filldraw[black] (0,0) circle (2pt) node[anchor=west] {(0,0)};
\filldraw[black] (1,1) circle (2pt) node[anchor=west] {(1,1)};
\filldraw[black] (2,2) circle (2pt) node[anchor=west] {(2,2)};}
{\draw[black, thick] (0,0) -- (2,2);}
\end{tikzpicture}
\end{center}
\caption{P-positions in $3$-Modular Nim.}
\label{2HeapMod3}
\end{figure}

These patterns of P-positions for $m = 3$ suggests a similar structure exists for other integers. This motivates a generalization to all odd integers.

\begin{theorem}\label{Odd2HeapPPosition}For odd $m$, a position of $m$-Modular Nim with $2$ heaps is a P-position if and only it is of the form $(i,i)$ for integers $i$, where $0 \leq i < m$. \end{theorem}

\begin{proof}
Let the set of specified positions be $\mathcal{S}$. We first prove that no move exists between any of the specified positions. Because no distinct positions have a heap size in common, no Type I move exists. Suppose that a Type II move exists between distinct positions $(i,i)$ and $(j,j)$ where $0 \leq i,j < m$. By Lemma~\ref{DominateTypeII}, we must have that:
\[2i \equiv 2j \pmod{m}. \]

Because $\gcd(2,m) = 1$, we may divide both sides by 2:
\[i \equiv j \pmod{m}. \]

Because $0 \leq i,j < m$, we have that $i = j$, contradicting the assumption of distinctness.

Now we must prove that for any position $(a,b) \notin \mathcal{S}$ there is a move to a position in $\mathcal{S}$.

\begin{enumerate}
\item Suppose that $\min(a,b) < m$, then a Type I move must exist.
\item On the other hand, if $\min(a,b) \geq m$, then $(a,b) \succ \mathcal{S}$. Therefore, a Type II move must exist because all residue classes modulo $m$ are covered by $\mathcal{S}$.
\end{enumerate}
\end{proof}

\section{$2$ Heap $m$-Modular Nim for Any $m$}\label{2HeapEven}

\subsection{Another Example: $m=6$.}

We consider an example of $m$-Modular Nim, where $m$ is an even integer: $m = 6$.

\begin{example}[$m = 6$] \label{m=6}
We claim that the set of P-positions for $m = 6$ is the following set:

\begin{center}
$(0,0)$ \\
$(1,1)$ \\
$(2,2)$
\begin{align*}
(3,4) & \hspace{40pt} (4,3) \\
(5,6) & \hspace{40pt} (6,5) \\
(7,8) & \hspace{40pt} (8,7).
\end{align*}
\end{center}

Let the specified set be $\mathcal{S}$. As before, we can manually verify that no Type I or Type II moves connect any two members of $\mathcal{S}$.

To show that $\mathcal{S}$ is the set of all P-positions, we must now show that for any position $(a,b) \notin \mathcal{S}$ there is a move to an element of $\mathcal{S}$.

\begin{enumerate}
\item If $\min(a,b) \leq 8$, then we may reach a member of $\mathcal{S}$ by removing the necessary number of tokens from the larger heap using a Type I move.

\item On the other hand, if $\min(a,b) > 8$, we have that $(a,b) \succ \mathcal{S}$. Because $\mathcal{S}$ contains positions with each possible total heap-sum modulo $6$, there must exist a Type II move.
\end{enumerate}
\end{example}

We can see the existence of three distinct groups of P-positions for the case where $m = 6$. The first group has both heaps equal in size and is the same as P-positions for 3-Modular Nim. The second and third group can be viewed as the P-positions of 3-Modular Nim shifted by $(3,4)$ and $(4,3)$ respectively.

This idea is further elucidated by Figure~\ref{2HeapMod6}. The red sections indicate the locations of the second and third groups that are shifted replicas of Figure~\ref{2HeapMod3}. This nesting structure is essential in finding a formula for the set of P-positions.

\begin{figure}[htbp]
\begin{center}
\begin{tikzpicture}[scale = 0.5]
\scriptsize
\filldraw[black] (0,0) circle (4pt) node[anchor=west] {(0,0)};
\filldraw[black] (1,1) circle (4pt) node[anchor=west] {(1,1)};
\filldraw[black] (2,2) circle (4pt) node[anchor=west] {(2,2)};
\draw[black, thick] (0,0) -- (2,2);
{
\filldraw[black] (3,4) circle (4pt) node[anchor=west] {(3,4)};
\filldraw[black] (4,3) circle (4pt) node[anchor=west] {(4,3)};
\filldraw[black] (5,6) circle (4pt) node[anchor=west] {(5,6)};
\filldraw[black] (6,5) circle (4pt) node[anchor=west] {(6,5)};
\filldraw[black] (7,8) circle (4pt) node[anchor=west] {(7,8)};
\filldraw[black] (8,7) circle (4pt) node[anchor=west] {(8,7)};
\draw[black, thick] (2,2) -- (3,4);
\draw[black, thick] (2,2) -- (4,3);
\draw[black, thick] (3,4) -- (7,8);
\draw[black, thick] (4,3) -- (8,7);

{
\filldraw[red] (3,4) circle (4pt) node[anchor=west] {(3,4)};
\filldraw[red] (4,3) circle (4pt) node[anchor=west] {(4,3)};
\filldraw[red] (5,6) circle (4pt) node[anchor=west] {(5,6)};
\filldraw[red] (6,5) circle (4pt) node[anchor=west] {(6,5)};
\filldraw[red] (7,8) circle (4pt) node[anchor=west] {(7,8)};
\filldraw[red] (8,7) circle (4pt) node[anchor=west] {(8,7)};
\draw[red, thick] (3,4) -- (7,8);
\draw[red, thick] (4,3) -- (8,7);
}

}
\end{tikzpicture}

\end{center}
\caption {P-positions for $6$-Modular Nim.}
\label{2HeapMod6}
\end{figure}

To make this nesting pattern more clear, we display the set of P-positions in $12$-Modular Nim in Figure~\ref{2HeapMod12}. The red section indicates one of the two embedded shifted replicas of P-positions in 6-Modular Nim as shown in Figure~\ref{2HeapMod6} whereas the green section shows the doubly nested and shifted copy of P-positions in $3$-Modular Nim. Note that because we are chiefly concerned with the recursive structure, the individual labels have been removed.

\begin{figure}[htbp]
\centering
\begin{tikzpicture}[scale = 0.25]
\filldraw[black] (0,0) circle (8pt) ;
\filldraw[black] (1,1) circle (8pt) ;
\filldraw[black] (2,2) circle (8pt) ;
\filldraw[black] (3,3) circle (8pt) ;
\filldraw[black] (4,4) circle (8pt) ;
\filldraw[black] (5,5) circle (8pt) ;
\filldraw[black] (6,7) circle (8pt) ;
\filldraw[black] (7,6) circle (8pt) ;
\filldraw[black] (8,9) circle (8pt) ;
\filldraw[black] (9,8) circle (8pt) ;
\filldraw[black] (10,11) circle (8pt) ;
\filldraw[black] (11,10) circle (8pt) ;
\filldraw[black] (12,15) circle (8pt) ;
\filldraw[black] (13,14) circle (8pt) ;
\filldraw[black] (14,13) circle (8pt) ;
\filldraw[black] (15,12) circle (8pt) ;
\filldraw[black] (16,19) circle (8pt) ;
\filldraw[black] (17,18) circle (8pt) ;
\filldraw[black] (18,17) circle (8pt) ;
\filldraw[black] (19,16) circle (8pt) ;
\filldraw[black] (20,23) circle (8pt) ;
\filldraw[black] (21,22) circle (8pt) ;
\filldraw[black] (22,21) circle (8pt) ;
\filldraw[black] (23,20) circle (8pt) ;
\draw[black, thick] (0,0) -- (5,5);
\draw[black, thick] (5,5) -- (6,7);
\draw[black, thick] (5,5) -- (7,6);
\draw[black, thick] (6,7) -- (10,11);
\draw[black, thick] (7,6) -- (11,10);
\draw[black, thick] (10,11) -- (12,15);
\draw[black, thick] (10,11) -- (14,13);
\draw[black, thick] (11,10) -- (15,12);
\draw[black, thick] (11,10) -- (13,14);
\draw[black, thick] (12,15) -- (20,23);
\draw[black, thick] (13,14) -- (21,22);
\draw[black, thick] (14,13) -- (22,21);
\draw[black, thick] (15,12) -- (23,20);

{
\filldraw[red] (6,7) circle (8pt) ;
\filldraw[red] (8,9) circle (8pt) ;
\filldraw[red] (10,11) circle (8pt) ;
\filldraw[red] (12,15) circle (8pt) ;
\filldraw[red] (14,13) circle (8pt) ;
\filldraw[red] (16,19) circle (8pt) ;
\filldraw[red] (18,17) circle (8pt) ;
\filldraw[red] (20,23) circle (8pt) ;
\filldraw[red] (22,21) circle (8pt) ;
\draw[red, thick] (6,7) -- (10,11);
\draw[red, thick] (10,11) -- (12,15);
\draw[red, thick] (10,11) -- (14,13);
\draw[red, thick] (12,15) -- (20,23);
\draw[red, thick] (14,13) -- (22,21);
}

>{
\filldraw[green] (15,12) circle (8pt) ;
\filldraw[green] (19,16) circle (8pt) ;
\filldraw[green] (23,20) circle (8pt) ;
\draw[green, thick] (15,12) -- (23,20);
}

\end{tikzpicture}
 \caption{P-positions for 12-Modular Nim.}
\label{2HeapMod12}
\end{figure}

\subsection{Potential P-positions}\label{sec:potentialpositions}

In order to formalize the previous notion of nesting, we will recursively define the positions that we later prove to be P-positions for any value of $m$.

\begin{definition}
\[\mathcal{A}_m = \left \{
\begin{array}{ll}
\hspace{5pt} (i,i) \text{ where } 0 \leq i < m & \hspace{112pt} : m \text{ is odd} \\
\hspace{5pt} (i,i) \text{ where } 0 \leq i < \frac{m}{2} &  \hspace{112pt}: m \text{ is even.}
\end{array}
\right. \]
We will refer to any member of $\mathcal{A}_m$ as a \emph{trunk position}.
\end{definition}

Note that the P-positions in $6$-Modular Nim (from Example \ref{m=6}) that appear before the ``splitting'' are trunk positions.

\begin{definition}
\[ \mathcal{B}_m = \left \{
\begin{array}{ll}
\hspace{5pt} \varnothing &: m \text{ is odd} \\
\begin{array}{l}\left(2a + \frac{m}{2} + 1, 2b + \frac{m}{2}\right) \\ \left(2a + \frac{m}{2}, 2b + \frac{m}{2} + 1\right) \end{array}
\text{ where } (a,b) \in  \mathcal{B}_{\frac{m}{2}}\cup \mathcal{A}_{\frac{m}{2}} &: m \text{ is even.}
\end{array}
\right. \]
We will refer to any member of $\mathcal{B}_m$ as a \emph{branch position}.
\end{definition}

Note that the P-positions in $6$-Modular Nim (from Example \ref{m=6}) that appear after the ``splitting'' are branch positions. Also note that both positions in $\mathcal{B}_m $ generated from $(a,b)$ have the same sum of coordinates.

We seek to prove that the set of P-positions in $m$-Modular Nim is $\mathcal{Q}_{m} = \mathcal{A}_{m} \cup \mathcal{B}_{m}$.

\begin{definition}  We call the set of positions $\mathcal{Q}_{m}$ \emph{potential positions}. \end{definition}

Note that the recursion that builds $\mathcal{Q}_{2m}$ from $\mathcal{Q}_{m}$ allows us to provide a recursion for the number of elements in $\mathcal{Q}_{m}$:
\[|\mathcal{Q}_{2m}| = m+ 2|\mathcal{Q}_{m}|.
\]

Because of the recursive doubling involved in generating trunk and branch positions, it is natural to consider the number of times that this doubling can occur. In other words, we wish to count the greatest power of $2$ that divides $m$.

\begin{definition}Define the \emph{$2$-adic order} of $m$ as the highest power of $2$ that divides $m$. We denote this arithmetic function as $\nu_2(m)$ \cite{NumberTheory}.\end{definition}

Now we prove a series of preliminary results.
The following lemma shows that a number may appear in a potential position in a particular coordinate at most once.

\begin{lemma} \label{2HeapTypeI}
If $(a,b)$ and $(a,c) \in \mathcal{Q}_m$, then $b = c$. \end{lemma}

\begin{proof}
We proceed by induction on the $2$-adic order of $k$. For the base case, suppose that $k$ is odd. From the definition of $\mathcal{Q}_k$, where $k$ is odd, no two distinct positions share a heap size, because all positions are of the form $(i,i)$ where $0 \leq i < k$.

For the inductive step, suppose that the lemma is true for $\mathcal{Q}_k$. We wish to prove that if $(a,b)$ and $(a,c) \in \mathcal{Q}_{2k}$, then $b = c$. We will split this into cases:

\begin{enumerate}
\item Suppose that both $(a,b)$ and $(a,c)$ are trunk positions. No two trunk positions share a heap size because they are of the form $(i,i)$ for different $i$.
\item Suppose one of the positions $(a,b)$ and $(a,c)$ is a trunk position and the other is a branch position. Any trunk position in $\mathcal{Q}_k$ is strictly dominated by $(k,k)$ and any branch position strictly dominates $(k,k)$. Therefore, a branch position and a trunk position cannot share a coordinate.
\item Suppose that both $(a,b)$ and $(a,c)$ are elements of $\mathcal{B}_{2k}$. If $a$ has the same parity as $k$, then $a = 2a^\prime + k$. Thus, we may write $b = 2b^\prime + k + 1$ and $c = 2c^\prime + k + 1$ where $\left(a^\prime, b^\prime\right)$ and $\left(a^\prime, c^\prime\right)$  are in $\mathcal{Q}_{k}$. By the inductive hypothesis, $b^\prime = c^\prime$. Therefore, $b = c$.

A similar symmetric argument works when $a$ has the opposite parity as $k$.
\end{enumerate}
Thus, we are done by induction.
\end{proof}

We can strengthen Lemma~\ref{2HeapTypeI} by showing that the set of integers allowed to be a coordinate of a potential position only consists of consecutive numbers.

\begin{lemma}\label{2HeapPPositionCount}
Integers between $0$ and $k\left(\frac{\nu_2(k)}{2} + 1\right) -1$ inclusive form the set of all first coordinates in the set of potential positions $\mathcal{Q}_k$. By symmetry the same is true for the second coordinate.\end{lemma}

\begin{proof}
As before, we induct on the $2$-adic order of $k$. Our base case is the odd integers. In Theorem \ref{Odd2HeapPPosition}, we defined the set of potential positions for odd integers $k$ to be all $(i,i)$ where $0 \leq i < k$. This concludes the base case, as $\nu_2(k)=0$.

For the inductive step, consider $\mathcal{Q}_{2k}$. For trunk positions, the first coordinate ranges over all numbers from $0$ to $k - 1$ inclusive. For branch positions, it ranges over all numbers of the form $2a+k$ and $2a+k+1$, where $a$ is a first coordinate of an element of $\mathcal{Q}_k$.

By induction, $a$ is between 0 and $k\left(\frac{\nu_2(k)}{2} + 1\right) -1$ inclusive. Thus the new branch positions for $2k$ have first coordinates in the range between $k$ and

\begin{align*}
2\left(k\left(\frac{\nu_2(k)}{2} + 1\right) -1 \right) + k + 1 &= \\
2k\left(\frac{\nu_2(k)}{2} + 1\right) + k - 1 &= \\
2k\left(\frac{\nu_2(k)}{2} + 1 + \frac{1}{2} \right) -1 &= \\
2k\left(\frac{\nu_2(2k)}{2} + 1  \right) -1.
\end{align*}

Therefore, an integer is between $0$ and $2k\left(\frac{\nu_2(k)}{2} + 1\right) -1$ inclusive if and only if it appears as a first coordinate in the set of potential positions $\mathcal{Q}_{2k}$, completing the induction.
\end{proof}

\begin{corollary}\label{thm:totalcount}
\[|\mathcal{Q}_m| = m\left(\frac{\nu_2(m)}{2} + 1\right).\] \end{corollary}

With Lemma~\ref{2HeapTypeI} and Lemma~\ref{2HeapPPositionCount}, we can now explicitly determine the number of potential positions.

The following lemma shows that two potential positions that share a remainder modulo $m$ have the same sum.

\begin{lemma}\label{2HeapTypeII} If $(a,b)$ and $(c,d) \in \mathcal{Q}_k$ and $a + b \equiv c + d \pmod{k}$, then $a + b = c + d$. \end{lemma}

\begin{proof}
We proceed by induction on $\nu(k)$. For the base case, suppose that $k$ is odd. All potential positions are of the
from $(i,i)$ where $0 \leq i < k$. Let $(a,b) = (m,m)$ and $(c,d) = (n,n)$. We can rewrite our condition as:
\[2m \equiv 2n \pmod{k}. \]

Because $\gcd(2,k) = 1$, we may divide both sides by $2$:
\[m \equiv n \pmod{k}. \]

Because $0 \leq m,n < k$, we have that $m = n$. Thus, $a + b = c + d$ where $k$ is odd.

For the inductive hypothesis, assume the lemma is true for $k$. We seek to now prove the corresponding statement for $2k$. We will divide this into cases:

\begin{enumerate}
\item Suppose that both $(a,b)$ and $(c,d)$ are trunk positions in $\mathcal{Q}_{2k}$. Because trunk positions are of the from $(i,i)$ where $0 \leq i < k$, we have that  $(a,b) = (m,m)$ and $(c,d) = (n,n)$ where $0 \leq m,n < k$. We can rewrite our condition as:
\[2m \equiv 2n \pmod{2k}. \]

Because $\gcd(2,2k) = 2$, we may divide both sides by $2$ only if we divide the modular base by $2$ as well:
\[m \equiv n \pmod{k}. \]

Because $0 \leq m,n < k$, this implies that $m = n$, finishing this case.

\item Suppose one of $(a,b)$ and $(c,d)$ is a trunk position and the other is a branch position in $\mathcal{Q}_{2k}$. The heap-sum of a trunk position is even and the heap-sum of a branch position is odd. Therefore, they cannot have the same remainder modulo an even number $2k$.

\item Suppose that both $(a,b)$ and $(c,d)$ are branch positions in $\mathcal{Q}_{2k}$. Therefore,  $(a,b)$ is equal to $(2a^\prime + k, 2b^\prime + k + 1)$ or $(2a^\prime + k + 1, 2b^\prime + k)$ where $\left(a^\prime, b^\prime\right)$ is a position in $\mathcal{Q}_{2k}$. In either case,
\[a + b = 2a^\prime + 2b^\prime + 2k + 1.\]

Similarly,
\[c + d = 2c^\prime + 2d^\prime + 2k + 1.\]

with $\left(c^\prime, d^\prime\right)$ in $\mathcal{Q}_{2k}$. We may now rewrite our given condition
\[a + b \equiv c + d \pmod{2k}\]
as
\[2a^\prime + 2b^\prime \equiv 2c^\prime + 2d^\prime \pmod{2k}.\]

Because $\gcd(2,2k) = 2$, we may divide both sides by $2$ only if we divide the modular base by $2$ as well:
\[ a^\prime + b^\prime \equiv c^\prime + d^\prime \pmod{k}. \]

By the inductive hypothesis, $a^\prime + b^\prime = c^\prime + d^\prime$. Therefore, $a + b = c + d$.

Thus, we are done by induction.
\end{enumerate}
\end{proof}

\begin{corollary}\label{cor:trunkremainders}
Trunk positions in $\mathcal{Q}_k$ have distinct remainders modulo $k$.
\end{corollary}

We are now ready to prove that every position not in $\mathcal{Q}_m$ has a move to a position in $\mathcal{Q}_m$ in the game of $m$-Modular Nim.

\begin{lemma}\label{MovesFromNotQtoQ}Every position not in $\mathcal{Q}_m$ has a move to an element of $\mathcal{Q}_m$ in the game of $m$-modular Nim.\end{lemma}

\begin{proof}
Let there be a position $(a,b) \notin \mathcal{Q}_k$. We will as usual induct on the $2$-adic order. Our base case is when $k$ is odd. By Lemma~\ref{Odd2HeapPPosition}, $\mathcal{Q}_k$ is the complete set of P-positions in $k$-Modular Nim. Therefore, all other positions are N-positions. From the definition of an N-position, every N-position has a move to a P-position.

For the inductive step, assume that this lemma is true for $\mathcal{Q}_k$ in the game of $k$-Modular Nim. We want to show the corresponding statement for $\mathcal{Q}_{2k}$ in the game of $2k$-Modular Nim. Again we will use casework. Assume without loss of generality that $a \leq b$.

\begin{enumerate}
\item Suppose that $0 \leq a < k$. Then, we have $(a,a) \in \mathcal{Q}_{2k}$. Therefore, there is a Type I move from $(a,b)$ to $(a,a)$, a member of $\mathcal{Q}_{2k}$.
\item Suppose that $k \leq a \leq b$ and $a \equiv b \pmod{2}$. The heap-sums of trunk positions in $\mathcal{Q}_{2k}$ form a set of $k$ even integers. By Corollary~\ref{cor:trunkremainders} there exists a trunk position with any given even remainder modulo $2k$. Therefore, there exists a trunk position with the same even remainder modulo $2k$ as $a+b$. As $(a,b)$ dominates $(k,k)$, it also strictly dominates any trunk position. Therefore by Lemma \ref{DominateTypeII}, there exists a Type II move from $(a,b)$ to the trunk position with the same remainder.

\item Suppose that $k \leq a \leq b$ and $a \equiv k \pmod{2}$ while $b \equiv k + 1 \pmod{2}$. Consider the position:
\[ \left(a^\prime, b^\prime\right) = \left(\frac{a - k}{2}, \frac{b - k - 1}{2} \right) \]

Note that $\left(a^\prime, b^\prime\right)$ is not an element of $\mathcal{Q}_k$, as $\left(2a^\prime + k, 2b^\prime + k + 1\right) = (a,b)$ is not an element of $\mathcal{Q}_{2k}$. Therefore, by the inductive hypothesis there exists a position $\left(q_1,q_2\right) \in \mathcal{Q}_{k}$ that can be reached with a Type II move from $\left(a^\prime, b^\prime\right)$ in the game of $k$-Modular Nim.

Thus, there exists a Type II move from $(a,b) = \left(2a^\prime + k, 2b^\prime + k + 1\right)$ to $\left(2q_1 + k,2q_2 + k + 1\right)$ in the game of $2k$-Modular Nim. Note that $\left(2q_1 + k,2q_2 + k + 1\right)$ is an element of $\mathcal{Q}_{2k}$. This finishes this case.

\item Suppose that $k \leq a \leq b$ and $a \equiv k + 1\pmod{2}$ while $b \equiv k  \pmod{2}$. We are done by a similar symmetric argument to the previous case.
\end{enumerate}
This completes the induction. \end{proof}

These results enable us to determine a recursive definition of P-positions in $m$-Modular Nim with $2$ heaps.

\begin{theorem}If $\mathcal{P}_m$ is the set of P-positions for $m$-Modular Nim with $2$ heaps, then $\mathcal{P}_m = \mathcal{Q}_m$. \end{theorem}

\begin{proof}We must first prove that no move exists between any two elements of $\mathcal{Q}_m$. Let $(a,b)$ and $(c,d)$ be distinct elements of $\mathcal{Q}_m$. If a Type I move exists between $(a,b)$ and $(c,d)$, we may assume without loss of generality that $a = c$. Thus $(a,b)$ and $(a,d)$ are distinct elements of $\mathcal{Q}_m$. However, from Lemma \ref{2HeapTypeI}, then we have that $b = d$ which contradicts distinctness.

Now suppose that a Type II move exists. By Lemma \ref{DominateTypeII}, we have that $a + b \equiv c + d \pmod{m}$. However, from Lemma \ref{2HeapTypeII}, we have that $a + b = c + d$. However, every move strictly decreases heap-sum, so no such move exists.

Moreover by Lemma \ref{MovesFromNotQtoQ}, every position not in $\mathcal{Q}_m$ has a move to an element of $\mathcal{Q}_m$. Therefore, $\mathcal{P}_m = \mathcal{Q}_m$. \end{proof}

We are now ready to describe the set of P-positions explicitly.

\section{Explicit Description of P-positions}\label{sec:explicit}

Let $m = k \cdot 2^n$ where $k$ is odd. The P-positions in the game of $m$-Modular Nim are built recursively from the positions of $k$-Modular Nim, by using the recursion described in Section~\ref{sec:potentialpositions}. The recursion is used $n$ times.

The recursion procedure is similar for all $m$ with the same 2-adic order. Figure~\ref{fig:8mod} depicts P-positions for 8-Modular Nim to emphasize the branching.

\begin{figure}[htbp]
\begin{center}
\begin{tikzpicture}[scale = 0.25]
\scriptsize
\filldraw[black] (0,0) circle (4pt);
\filldraw[black] (1,1) circle (4pt);
\filldraw[black] (2,2) circle (4pt);
\draw[black, thick] (0,0) -- (3,3);
{
\filldraw[black] (3,3) circle (4pt);
\filldraw[black] (4,5) circle (4pt);
\filldraw[black] (5,4) circle (4pt);
\draw[black, thick] (3,3) -- (4,5);
\draw[black, thick] (3,3) -- (5,4);

\filldraw[black] (6,7) circle (4pt);
\filldraw[black] (7,6) circle (4pt);

\filldraw[black] (8,11) circle (4pt);
\filldraw[black] (9,10) circle (4pt);
\filldraw[black] (10,9) circle (4pt);
\filldraw[black] (11,8) circle (4pt);

\filldraw[black] (12,19) circle (4pt);
\filldraw[black] (13,18) circle (4pt);
\filldraw[black] (14,17) circle (4pt);
\filldraw[black] (15,16) circle (4pt);

\filldraw[black] (19,12) circle (4pt);
\filldraw[black] (18,13) circle (4pt);
\filldraw[black] (17,14) circle (4pt);
\filldraw[black] (16,15) circle (4pt);

\draw[black, thick] (5,4) -- (7,6);
\draw[black, thick] (4,5) -- (6,7);

\draw[black, thick] (7,6) -- (11,8);
\draw[black, thick] (7,6) -- (9,10);
\draw[black, thick] (6,7) -- (8,11);
\draw[black, thick] (6,7) -- (10,9);

\draw[black, thick] (8,11) -- (12,19);
\draw[black, thick] (8,11) -- (16,15);
\draw[black, thick] (11,8) -- (19,12);
\draw[black, thick] (11,8) -- (15,16);
\draw[black, thick] (9,10) -- (13,18);
\draw[black, thick] (9,10) -- (17,14);
\draw[black, thick] (10,9) -- (18,13);
\draw[black, thick] (10,9) -- (14,17);
}
\end{tikzpicture}

\end{center}
\caption {P-positions for $8$-Modular Nim.}
\label{fig:8mod}
\end{figure}

\begin{definition}\label{Branch Level}The \textit{$i$-level branch} is the set of P-positions that are generated from the trunk of $\mathcal{Q}_{k\cdot 2^{n-i}}$ using $i$ splitting procedures as described in Section~\ref{sec:potentialpositions}.\end{definition}

The P-positions in the $i$-level branch are of the form
\[(2^ia + f_im + b_1, 2^ia + f_im + b_2),\]
where $f_i$ is the coefficient by which $m$ is multiplied after completing $i$ splitting procedures. Moreover, we have $0 \leq a < k \cdot 2^{n-i-1}$, for $i < n$, and $0 \leq a < k$, for $i=n$.

Let us trace through recursion and find $f_i$. The $i$-level branch is recursively generated from the trunk of $\mathcal{Q}_{\frac{m}{2^i}}$. After the first recursion the coefficient is $\frac{m}{2^i}$. After the next recursion it is \[2\cdot \frac{m}{2^i} + \frac{m}{2^{i-1}}=2 \cdot \frac{m}{2^{i-1}}.\] Repeating again we get \[2\cdot \frac{m}{2^{i-2}} + \frac{m}{2^{i-2}}=3 \cdot \frac{m}{2^{i-2}}.\] As we continue, we see that $f_i = \frac{i}{2}$.

We also can calculate that the largest possible value for $b_1$ is $2^i-1$, and the smallest is 0. The same is true for $b_2$. Therefore, the following lemma follows.

\begin{lemma}
The P-positions within the $i$-level branch correspond to all possible values for $0 \leq b_1 < 2^i$ and $b_2 = 2^i-1-b_1$.
\end{lemma}

\begin{proof}
The P-positions corresponding to the same $a$ within the $i$-level branch have the same heap-sum. They all have different first coordinates and there are $2^i$ of them.\end{proof}

Combining all the results together, we deduce the following theorem.

\begin{theorem}
The P-positions of $m$-Modular Nim are of the form
\[\left(2^ia + \frac{im}{2} + b, 2^ia + \frac{im}{2} + 2^i-1-b\right),\]
where $0 \leq i \leq n$ and $0 \leq b < 2^i$. In addition, $0 \leq a < k \cdot 2^{n-i-1}$, for $i < n$, and $0 \leq a < k$, for $i = n$.
\end{theorem}

Note that 0-level branch corresponds to the trunk and the formula correctly produces positions of the form:
\[(a,a)\]
for $0 \leq a < k \cdot 2^n = m$. Similarly, the 1-level branch positions are of the form:
\[\left(2a + \frac{m}{2} + 1, 2a + \frac{m}{2}\right) \text{ and } \left(2a + \frac{m}{2}, 2a + \frac{m}{2}+1\right),\]
for $0 \leq a < k \cdot 2^{n-1}$.

We see that in the list of all P-positions the range of the coordinates with the higher heap-sum is above the range of the coordinates with the lower heap-sum.

\begin{corollary}
Any P-position with a greater sum dominates every P-position with a lower sum.
\end{corollary}

In addition, on level $i$ the first coordinate ranges from $f_im=i\frac{m}{2}$ to \[f_im+\frac{m}{2}-1=(i+1)\frac{m}{2}-1 = f_{i+1}m- 1.\]
This further confirms Lemma~\ref{2HeapPPositionCount} in that the numbers in the first and second coordinates of P-positions are consecutive. We can also see that the largest number in the range is $2^nk+n\frac{m}{2}-1$, which, not surprisingly, matches Corollary~\ref{thm:totalcount}.

The description above allows us to calculate the second coordinate of the P-position, given the first coordinate.

\begin{lemma}
Suppose the first coordinate of an $m$-Modular Nim P-position is $x$, then the second coordinate is
\[ 2^{\left \lfloor \frac{2x}{m}\right  \rfloor+1} \left \lfloor \frac{x-\left \lfloor \frac{2x}{m}\right  \rfloor \frac{m}{2}}{2^i}\right \rfloor + \left \lfloor \frac{2x}{m}\right \rfloor m + 2^{\left \lfloor \frac{2x}{m}\right \rfloor}-1 -x.\]
\end{lemma}

\begin{proof}
Suppose the first coordinate is $x$, then the level $i$ is defined as $\left\lfloor \frac{2x}{m}\right\rfloor$. After that we can calculate $a$ in the formula as  $\left\lfloor \frac{x-im/2}{2^i}\right\rfloor$. Therefore, the second coordinate is:
\[ 2^{i+1}a + im + 2^i-1 -x.\]\end{proof}

\section{$m$-Modular Nim for Odd $m$ and any Number of Heaps}\label{sec:manyheaps}

We will now generalize the results from  Theorem~\ref{Odd2HeapPPosition} to any number of heaps.

\begin{theorem}Let $a = (a_1,a_2,\ldots,a_n)$ be a position in $m$-Modular Nim with $n$ heaps, where $m$ is odd. It is a P-position if and only if:

\begin{enumerate}
\item $|a| < 2m$
\item $\oplus_i a_i = 0$.
\end{enumerate}\end{theorem}

\begin{proof}
Let the specified set be $\mathcal{Q}$. We will first show that no move exists between any two positions in $\mathcal{Q}$. By Bouton's Theorem, no Type I move exists.

Because the heap-sum of all members of $\mathcal{Q}$ is less than $2m$, any Type II move must subtract exactly $m$ tokens total. However, all members of $\mathcal{Q}$ have an even heap-sum as the bitwise XOR of their heap sizes is $0$. Because $m$ is odd, no two elements of $\mathcal{Q}$ can be connected by a Type II move.

We must now show that every position $p=(p_1,p_2,\ldots,p_n)$ not in $\mathcal{Q}$ has a move to an element of $\mathcal{Q}$.

\begin{enumerate}
\item Suppose that $p$ is not a P-position in Nim and there is a Nim-move from $p$ to a position $p^\prime$ such that $|p^\prime| < 2m$. Then $p^\prime \in \mathcal{Q}$ and we have found our move.
\item Suppose that $p$ is not a P-position in Nim; and there is a Nim-move from $p$ to a position $q$ such that $|q| \geq 2m$. There exists a number $k$, such that $|p| -km < 2m$ and $|p| -km$ is even. In the game of Nim, there exists an optimal play in which only $1$ token is removed per turn. That means, for any P-position $q$, there exists a P-position dominated by $q$ for every even heap-sum less than $|q|$. In particular, there exists a P-position in Nim, $p^\prime$, such that $|p^\prime| = |p| -km$ and $q$ dominates $p^\prime$. Therefore, $p$ dominates $p^\prime$, and $p^\prime$ is reachable from $p$ by a Type II move.
\item Suppose that $p$ is a P-position in Nim, and $|p| \geq 2m$. A similar argument to the one above shows that there exists a P-position of Nim $p^\prime \in \mathcal{Q}$ that is reachable from $p$ via a Type II move.
\end{enumerate} \end{proof}

\section{Acknowledgements}
The authors are grateful to the MIT-PRIMES program for supporting this research and to Prof. Aviezri Fraenkel for helpful suggestions.