\section{Introduction}
Control over the propagation of light using artificial photonic crystals\cite{PhC_John,PhC_Yablonovitch}
and electromagnetic metamaterials\cite{MetaSmith} has received enormous attention in recent decades mainly
because of its importance for many applications in the vibrant field of photonics. For example, metamaterials
such as left-handed media\cite{Veselago68,Chan18} have shown promising potential for novel
technologies\cite{Pendry00,Klimov12}. In recent decades, great progress in this field
has been often made by taking advantages of analogies with electronic systems in solid state physics.
For example, the concept of a photonic band gap material\cite{PhC_John,PhC_Yablonovitch}, a man-made system
with a periodic dielectric function, was inspired by the electronic Bloch states in a crystalline semiconductor.

More recently, there have been growing interests in using topological photonic materials\cite{Lu2014}
to manipulate the flow of light, again inspired by the recent developments of electronic topological
materials \cite{Hasan2010,Qi2011,hmwengreview}. In particular, the electronic quantum anomalous Hall (QAH) phase
is a two-dimensional (2D) bulk ferromagnetic insulator with a nonzero Chern number in the presence of spin-orbit
coupling (SOC) but in the absence of applied magnetic fields\cite{Haldane1988,QAHexp}. Its associated metallic
chiral edge states in this Chern insulator carry dissipationless unidirectional electric current. Haldane and Raghu
recently proposed\cite{photon_QAH_Haldane} to construct analogs of this intriguing QAH in photonic crystals made
of time-reversal symmetry (TRS) breaking materials to realize unidirectional optical waveguides. Subsequently,
these topological electromagnetic states in a number of gyromagnetic photonic crystals with broken TRS were
further proposed \cite{QAH_Wangzhen_theory,ochiai2009photonic,QAH_poo2011,zeromode,QAH_largeChern,Dirac_mirrorsym,Chan2018}
and observed\cite{QAH_Wangzhen_exp,QAH_poo2013,QAH_largeChern_exp,He2016}. Interestingly, photonic analog of
electronic quantum spin Hall effect in 2D topological insulators with TRS\cite{Hasan2010,Qi2011} were also
observed in bi-anisotropic photonic crystals\cite{Khanikaev2013,z2new,z2}.

Nevertheless, investigations of photonic unidirectional edge modes have mostly been limited to topological
photonic crystals periodic on the scale of the operational wavelength, and this considerably restricts the
applications of topological photonic materials. Very recently, Gao et. al.\cite{chiralhypermeta} theoretically
demonstrated topological photonic phase in chiral hyperbolic metamaterials (CHM) made of continuous TRS media
with photonic edge states robust against disorder on all length scales\cite{chiralhypermeta}.
In hyperbolic metamaterials\cite{HMM13},
which are plasmonic metamaterials, equi-frequency surfaces (EFSs) of transverse electric-field (TE) and
transverse magnetic-field (TM) modes are degenerate on the high-symmetry points in the momentum space [see Fig. 1(b)].
When the bi-anisotropic property (chirality) is introduced in hyperbolic metamaterials with TRS,
which then become CHMs, the degeneracies are broken due to the coupling between TE and TM modes\cite{Khanikaev2013}
and consequently, a nontrivial band gap is opened\cite{chiralhypermeta}. The nontrivial topology of the CHM
results from the nonzero Berry curvature due to the chirality (equivalent to the SOC in electronic topological
insulators) and broken spatial inversion symmetry in continuous medium\cite{continuousmediachern}. Furthermore,
the design principle was based on the concept of  a Floquet topological insulator (FTI)\cite{floquetrechtsman2013}
where by treating the distance $z$ propagated by a waveguide mode as a timelike coordinate, e.g.,
a honeycomb lattice of helical waveguides behaves like a FTI.

In this work, as a novel mechanism for controlling light flow, we introduce the photonic Chern insulators
made of continuous gyromagnetic hyperbolic metamaterials (GHM) with the band gap opened by the TRS-breaking
gyromagnetic effect [Figs. 1(b) and (c)]. The non-trivial topology is demonstrated by the calculated Berry
curvature and nonzero Chern number due to the broken TRS. The unidirectional backscattering-immune non-radiative
edge modes at the interface between the GHM and vacuum are uncovered by the finite-element electromagnetic simulations.
As in the photonic quantum spin-Hall insulators (QSHI) made of the CHMs\cite{chiralhypermeta}, the photonic
edge states in our Chern insulators are robust against disorder on all length scales. In contrast to the
QSHIs\cite{chiralhypermeta}, however, the light flow direction on the surface of our Chern insulators
can be easily flipped by switching the direction of applied magnetic fields. Furthermore, our Chern insulators
made of the GHM can be easily fabricated.

\section{Gyromagnetic hyperbolic metamaterial}

We consider a GHM as a hyperbolic metamaterial\cite{HMM13} with the gyromagnetic response
which is described by the constitutive relation,
\begin{equation}
\left( \begin{matrix}
\mathbf{D}  \\
\mathbf{B}  \\
\end{matrix} \right)=\left( \begin{matrix}
{{\varepsilon_0}}\hat{\varepsilon} & 0  \\
0 & {{\mu_0}}\hat{\mu}  \\
\end{matrix} \right)\left( \begin{matrix}
\mathbf{E} \\
\mathbf{H}
\end{matrix} \right),\\
\end{equation}
where
\begin{equation}
\hat{\varepsilon}=\left(
\begin{array}{ccc}
	\epsilon_{xx}  & 0 & 0 \\
	0 & \epsilon_{yy} & 0 \\
	0 & 0 & \epsilon_{zz} \\
\end{array}
\right),
\hat{\mu }=
\begin{pmatrix}
	\mu_{xx} & -i\gamma & 0 \\
	i\gamma & \mu_{yy} & 0 \\
	0 & 0 & \mu_{zz}
\end{pmatrix}
\end{equation}
are the relative permittivity and permeability tensors, respectively, and $\gamma$ is the gyromagnetic parameter,
representing the degree of TRS breaking upon application of a magnetic field in the $z$-direction (gyromagnetic effect).
For simplicity, let us set $\mu_{xx}=\mu_{yy}=\mu_{zz}=\mu=1$, $\epsilon_{xx}=\epsilon_{yy}=\alpha$ and
$\epsilon_{zz}=\beta$ such that EFS dispersion are isotropic in $X-Y$ plane.
The propagation behavior in the gyromagnetic hyperbolic metamaterials can be described by the wave equation
of the electric field $\mathbf{E}=(E_x,E_y,E_z)^T$ as
\begin{equation}
\label{eq_wave}
\mathbf{k}\times \hat{\mu}_r^{-1} \mathbf{k}\times\mathbf{E}+k_0^2 \hat{\epsilon_r} \mathbf{E}=0,
\end{equation}
where $\mathbf{k}=(k_x,k_x,k_z)$ is the wave vector, and $k_0=\omega/c$ is the wave number in vacuum. The left
side of Eq. \eqref{eq_wave} can be rewritten in the matrix form as
\begin{equation}
\label{eq_waveMatrix}
\left(
\begin{array}{ccc}
\alpha k_0^2 - \frac{k_y^2}{\mu}-Ak_z^2  & \frac{k_xk_y}{\mu}+Bk_z^2 & Ak_xk_z-Bk_yk_z \\
\frac{k_xk_y}{\mu}-Bk_z^2 & \alpha k_0^2 - \frac{k_x^2}{\mu} - Ak_z^2 & Ak_yk_z+Bk_xk_z\\
Ak_xk_z+Bk_yk_z & Ak_yk_z-Bk_xk_z & \beta k_0^2-Ak_t^2\\
\end{array}
\right)
\end{equation}
where $k_x^2+k_y^2=k_t^2$, $A=\frac{\mu}{\mu^2-\gamma^2}$, $B=\frac{\gamma}{\mu^2-\gamma^2}$. The nontrivial
solutions of $\mathbf{E}=(E_x,E_y,E_z)^T$ exist when the determinant of Eq.~\ref{eq_waveMatrix} equals to zero,
resulting in the eigen equation of
\begin{equation}
\label{eq_Det0}
\alpha ^2 \beta  k_0^2 \mu-(\alpha A\mu+\beta-A)\alpha k_t^2
+(A-2\alpha \mu A)\beta k_z^2=0.
\end{equation}

Starting with an isotropic optical medium with $\epsilon_{xx}=\epsilon_{zz}=2$ and $\gamma=0$, one has two routes
to arrive at a GHM, as illustrated in Fig. 1. The EFS dispersion of the isotropic medium is a perfect sphere,
as shown in Fig. 1(a). Furthermore, the eigenstates consisting of transverse electric mode TE ($E_x=E_y=0$, $H_z\neq0$)
and transverse magnetic modes TM ($H_x=H_y=0$, $E_z\neq0$), are degenerate and their dispersion are identical
and on the same sphere. When  the isotropic medium is transformed to a hyperbolic metamaterial
($\epsilon_{xx}=2$, $\epsilon_{zz}=-1$, and $\gamma=0$) by varying $\epsilon_{zz}$ from 2 to -1, the EFS sphere
of the TM mode splits and becomes two parabola along the $k_z$-axis while that of the TE mode remains spherical,
as depicted in Fig. 1(b) [see Fig. 2(a, b) for polarization analysis]. This is because TM "sees" $\epsilon_{xx}=\epsilon_{yy}=2$
only while, in contrast, TM "sees" $\epsilon_{x-y}=2$ and =$\epsilon_{zz}=-1$. Interestingly, the two TM parabola
touch the TE sphere at $k_z =\pm 1.0k_0$, respectively. When the gyromagnetic effect is further introduced into
the system by making $\gamma$ nonzero (e.g., $\gamma=0.8$), the system becomes a GHM. As a result, a band gap
is opened at both degenerate points $k_z =\pm1.0k_0$ [see Fig. 1(d)]. Interestingly, since the gyromagnetic effect breaks TRS,
the eigenstates near $\pm 1.0k_0$ become circularly (elliptically) polarized, while eigenstates on the upper and
lower hyperbolic bands far away from the singular points remain TM-polarized [see Fig. 2(c,d)].

\begin{figure}[htbp]
\includegraphics[width=3.2in]{figure1.pdf}
\caption{EFS evolution from (a) an isotropic medium ($\epsilon_{xx}=\epsilon_{zz}=2$ and $\gamma=0$) to (d)
a GHM ($\epsilon_{xx}=2$,$\epsilon_{zz}=-1$, and $\gamma=0.8$) via either (b) a hyperbolic metamaterial
($\epsilon_{xx}=2$,$\epsilon_{zz}=-1$, and $\gamma=0$) or (c) a gyromagnetic medium ($\epsilon_{xx}=\epsilon_{zz}=2$
and $\gamma=0.8$). TE (TM) and RCP (LCP) denote TE (TM) polarization and right-handed (left-handed) circular
polarization, respectively.}
\end{figure}

Alternatively, one can go from an isotropic medium first to a gyromagnetic medium (e.g., $\epsilon_{xx}=\epsilon_{zz}=2$
and $\gamma=0.8$) [see Fig. 1(c)] and then to the GHM. In a gyromagnetic medium, because of broken TRS, TE and TM
modes are coupled, resulting in the separation of the degenerate EFS spheres into two ellipsoids [Fig. 1(c)].
Also eigenstates becomes elliptically polarized [see Fig. 2(e,f)]. In particular, the eigenstates on the top
and bottom surfaces of the ellipsoids become nearly fully circularly polarized, while those in the vicinity
of the $k_z=0$ plane remain almost purely either (inner ellipsoid) TE or (outer ellipsoid) TM polarization.
Finally, when $\epsilon_{zz}$ is tunned to a negative value (e.g. $\epsilon_{zz}=-1$), the system becomes a GHM
and the outer ellisoid splits to form two open parabolas, as shown in Fig. 1(d).

\begin{figure}[htbp]
\includegraphics[width=3.2in]{figure2.pdf}
\caption{Polarization analysis of eigenstates. Calculated polarization of EFS dispersion for (a, b) the hyperbolic
metamaterial, (c, d) the GHM (c, d) and (e, f) the gyromagnetic medium. RCP (LCP) and TE (TM) denote
right-handed (left-handed) circular polarization and TE (TM) polarization, respectively.}
\end{figure}

\section{Topological phase transition}

To investigate the topological property of the EFS dispersions of all four kinds of optical media and metamaterials
and also to examine the topological nature of the band gaps near $k_z=\pm k_0$,
we calculate the Berry phase $\Phi =\iint{\mathbf{F}(\mathbf{k})\cdot ds}$ and hence the Chern number $\mathcal{C}=\tfrac{1}{2\pi}\Phi$
of all EFS surfaces.\cite{fukui2005chern,Chan2018} Following our previous work\cite{Chan2018},
we adopt the efficient numerical algorithm reported in Ref. [\onlinecite{fukui2005chern}] to
evaluate the Berry curvature $\mathbf{F}(\mathbf{k})$. Note that in principle, the Berry phase and the Chern number are
well defined only for a closed surface. Fortunately, our test calculations show that the Berry curvature $\mathbf{F}(\mathbf{k})$
is negligibly small when the radial wave vector $k_\rho=\sqrt{k_x^2+k_y^2}>15$.
Therefore, although the hyperbolic shaped EFS surfaces are open (Figs. 1 and 2), we find that the calculated Chern numbers would converge
well to integers as long as the surface integration is carried out from the $k_0$ up to the radial wave vector $k_\rho=\sqrt{k_x^2+k_y^2}$
being larger than $15$. The Chern numbers calculated in this way are shown in Fig. 1, and the Berry curvature distributions
are displayed in Fig.~3. We notice that the Berry curvature has the following symmetry properties:
(a) $\mathbf{F}(\mathbf{k})=\mathbf{F}(-\mathbf{k})$ if the system has the spatial inversion symmetry (IS) and
(b) $\mathbf{F}(\mathbf{k})=\mathbf{-F}(-\mathbf{k})$ if the system has the TRS symmetry. Therefore, the Berry curvature
is identically zero and hence the Chern number is zero for both isotropic medium and hyperbolic metamaterial because they have
both IS and TRS [See Figs.~1(a) and 1(b)].

\begin{figure}[htbp]
\includegraphics[width=3.2in]{figure3.pdf}
\caption{Berry curvature distribution on EFSs for the gyromagnetic hyperbolic metamaterials with
(a) $\alpha=\beta=2$, and $\gamma=0.8$, (b) $\alpha=2$, $\beta=-1$, and $\gamma=1.2$.}
\end{figure}

In contrast, the Berry curvatures become nonzero in both gyromagnetic medium and GHM since their TRS is broken,
as shown in Fig.~3. Nevertheless, in the gyromagnetic medium, the Berry curvatures of the eigenstates on each EFS ellipsoid
form an odd function of $k_z$ with respect to the out-surface element [Fig. 3(a)]. Consequently, the sum of
the Berry curvatures on the EFS surface of each ellipsoid (i.e., the Berry phase and the Chern number) remains zero [Fig. 1(c)].
This is also true for the inner ellipsoid in the GHM [Fig. 3(b)] and hence the Chern number of the inner EFS is zero [Fig. 1(d)].
However, the outer ellipsoid in the GHM now splits and hence transforms into two open parabolic surfaces with the Berry curvatures
pointing to the positive $z$-direction. Although the Berry curvatures on the upper and lower surfaces still form
an odd function of $k_z$, they are on the different surfaces. Resultantly, the upper and lower open surfaces each
acquires a Chern number of -1. Therefore, the phase transformation from the hyperbolic to gyromagnetic hyperbolic
due to the gyromagnetic effect is a topological one, and the band gaps centered at $\pm k_0$ are topological nontrivial.

\section{Nonradiative one-way photonic edge mode}
Next, we follow the method suggested in Ref. [\onlinecite{Dya88}] to calculate analytically the surface band between the dielectric and
anisotropic material. Let the interface between air ($y>0$) and GHM ($y<0$) be on the $xz$ plane ($y=0$). The wave vector in
GHM normal to the interface can then be obtained by solving Eq. \eqref{eq_Det0}

\begin{widetext}
\begin{equation}
\label{eq_ky}
k_y^{GHM} = \pm \left( \frac{\alpha  k_0^2 \left(\alpha  \mu ^2+\beta  \left(\mu
	^2-\kappa ^2\right)\right)-2 \alpha  \mu  k_x^2-\alpha
	\mu  k_z^2-\beta  \mu  k_z^2 \pm Y}{2 \alpha  \mu }\right) ^{0.5}.
\end{equation}
where
\begin{equation}
Y \equiv \sqrt{\alpha ^2 k_0^4 \left(\alpha  \mu ^2+\beta
	\left(\kappa ^2-\mu ^2\right)\right)^2+2 \alpha  k_0^2
	\mu  k_z^2 \left(-\alpha ^2 \mu ^2+\alpha  \beta
	\left(\kappa ^2+2 \mu ^2\right)+\beta ^2 \left(\kappa
	^2-\mu ^2\right)\right)+\mu ^2 (\alpha -\beta )^2 k_z^4}
\end{equation}

Then two linear independent eigen fields in GHM $\vec{E}_i(k_x,k_y,k_z)$ ($i=3,4$) can be obtained by solving the null space
of Eq. \eqref{eq_waveMatrix} with $k_y=k_y(k_x,k_z)$ from Eq. \eqref{eq_ky}. Once $\vec{E}_i$ are obtained, $\vec{H}_i$
can also be obtained via Faraday's law $\vec{H}_i=\frac{1}{Z_0}\vec{k} \times \vec{E}_i$. Likewise, two orthogonal eigen fields
in the air can be expressed as
\begin{equation}
\vec{E}_1=\left(
\begin{array}{c}
-k_z \\
0 \\
k_x \\
\end{array}
\right), \quad
\vec{E}_2=\left(
\begin{array}{c}
-k_y \\
k_x \\
0 \\
\end{array}
\right), \quad
\vec{H}_1=\frac{1}{Z_0}\left(
\begin{array}{c}
k_x k_y \\
-k_x^2-k_z^2 \\
k_y k_z \\
\end{array}
\right), \quad
\vec{H}_2=\frac{1}{Z_0}\left(
\begin{array}{c}
-k_x k_z \\
-k_y k_z \\
k_x^2+k_y^2 \\
\end{array}
\right),
\end{equation}
where $Z_0=\sqrt{\mu_0/\epsilon_0}$ is the vacuum impedance, and $k_y=i\sqrt {k_x^2 + k_z^2 - k_0^2}$. Since the tangential
components of electromagnetic fields should be continuous across the interface, we arrive at $Det(D)=0$ with
\begin{equation}
D(k_x,k_z)=
\left(
\begin{array}{cccc}
-k_z & -k_y & E_{3x} & E_{4x} \\
k_x & 0 & E_{3z} & E_{4z} \\
k_x k_y & -k_x k_z & H_{3x} & H_{4x} \\
k_y k_z & k_x^2+k_y^2 & H_{3z} & H_{4z} \\
\end{array}
\right).
\end{equation}
In order to obtain a surface mode, it is important to choose the range of the imaginary part of $k_y$ ($k_y^{GHM}$)
in the vacuum (GHM) to be positive (negative). In other words, the fields in the air and GHM must decay exponentially
along the vertical direction of the interface.
\end{widetext}

According to the bulk-edge correspondence, the gap Chern numbers of $\Delta C=-1$ implies that there is
one propagating edge state in each band gap for each air-GHM interface. Furthermore, these edge states
are chiral and topology-protected, i.e., they are reflection-free one-way edge states. To verify these
amazing predictions, we solve the Maxwell’s equations for a slab made of the GHM, as illustrated in Fig.~4(a).
The calculated EFS dispersions for this slab, displayed in Fig. 4(b), show that in addition to the bulk bands
(black curves), there is indeed one edge state for both top (red curve) and bottom (blue curve) interface
between air and the GHM slab in each band gap.

We further simulate light propagation in the air-GHM interface by putting a line source along the $z$-axis
at $x = 0$. The results are displayed in Figs. 4(c)-4(f). Our simulations show clearly that light on
each interface runs without radiation in either $+x$ or $–x$ direction but not both directions, i.e.,
light propagation on the edges is unidirectional. We notice that in Fig.~4(c), both edges at $k_z=1.6k_0$
carry light towards $–x$ direction. This is because the dispersions of both bulk and edge states are
symmetric with respect to the mirror located at $k_z=0$. And the propagation along the $–x$ direction is
caused by the negative gap Chern numbers ($\Delta C=-1$) of both band gaps. Moreover, Fig.~4(g) shows
that light can overcome the obstacles on the surface, indicating light propagation is reflection-free
due to the topological protection. Finally, the direction of light propagation depends on the sign
of gyromagnetic parameter $\gamma$ [see Figs. 4(c) and 4(d)], demonstrating that one can control light
transmission direction by switching the magnetization direction of the GHM slab.

\begin{figure
\includegraphics[width=3.2in]{figure4.pdf}
\caption{One-way propagating edge states. (a) A GHM slab in air. (b) Calculated EFS dispersions for
the GHM slab in air (a). Black, red, and blue lines denote dispersions of bulk state, top and bottom edge states,
respectively. The dashed green circle denotes light cone. (c) and (d) Simulated propagations of light emitted
by a line source at $x=0$ on the top surface at $k_z=1.6k_0$ with $\gamma=-0.8$ and $\gamma= +0.8$, respectively.
(e) The same as in (c) except that the surface is now uneven, which does not stop light propagate along the negative $x$-direction.
(f) The same as in (d) but on the bottom surface.
}
\end{figure}

\begin{figure
\includegraphics[width=3.2in]{figure5.pdf}
\caption{(a) EFS dispersions and (b) decay rate of the electromagnetic waves for the GHM slab in the air.
In (a), the solid black and dashed blue curves denote the dispersions of the bulk states
and the air, respectively. In (a) and (b), the yellow regions denote the overlap regions
of the air dispersion and bulk gap regions, while the grey regions denote the complete gap regions.
(c) and (d) Electric field distributions of the edge states with $k_z=1.06$ and $k_z=0.88$, respectively.
Note that in (c), the electric field is completely confined to the edge of the GHM slab and
hence there is no radiative energy loss, because the $k_z$ is outside the light EPS sphere.
In contrast, in (d), because the $k_z$ is within the light EPS sphere, the electromagnetic wave
of the edge state propagates strongly into the air, resulting in a large radiative energy loss.
}
\end{figure}

A GHM is based on a hyperbolic metamaterial, which is generally a
bilayer superlattice consisting of a metal layer and a dielectric layer stacked
along the normal direction (see next Section for details). Thus, the GHMs
possess the property of surface plasmon polaritons (SPPs). Compared with
the electromagnetic waves in the air, the SPPs are slower waves and hence
have larger wave vectors than that of light in the air. As a result, the edge
states on a GHM may inherit this larger wave vector propagation, as shown in
Fig. 5. Therefore, there may be no intersection between the air light line and
edge states, i.e., the edge states are out of the light line. Consequently, the
edge state would propagate on the interface without radiation [Fig. 5(c)] because the
edge states cannot couple to the electromagnetic waves in the air. In strong
contrast, most of photonic topological insulators are made of photonic
crystals and consequently, the one way edge states on the interfaces of these insulators
suffer from severe radiation loss and a metal wall serving as a cladding perfect
electric conductor or a perfect magnetic conductor has to be added to the
surrounding edges to stop radiation loss into the air \cite{QAH_Wangzhen_exp}.
In practice, there is no lossless perfect conductor and even good metals such as copper
used as the cladding walls absorb the radiation significantly. This certainly hinders
their applications. In contrast, in a GHM, the air serves as the good insulator for
stopping the radiation of the edge states into the air [Fig. 5(c)].

\section{Negative refraction of edge mode}
Another interesting property of the topological surface states of the GHM
is negative refraction, similar to that of the CHM reported recently in Ref. [\onlinecite{Klimov18}].
As shown in Fig. 6(b), the $x$ component of the group velocity (energy flow) of all the surface states
of the GHM with $\gamma=-0.8$ (the blue line in the bottom panel) is positive ($v_{v,x}>0$),
while the phase velocity ($v_{ph,x} = \omega/k_x$) changes sign across $k_z=0$. As a result,
in the $k_x<0$ region the $x$ components of the group velocity and phase velocity have opposite signs
($v_{g,x}>0$ and $v_{ph,x}<0$) [see Fig. 4(b)]. This is the signature
of negative refraction \cite{HMM13,Klimov18}. On the contrary, the surface states on the GHM with $\gamma=0.8$
in the $k_x<0$ region have $v_{g,x}<0$ and $v_{ph,x}<0$ [see upper panel of Fig. 6(b)].
Following the method proposed in Ref. [\onlinecite{He18}],
we show in Fig. \ref{figNR} that by drawing a vertical line in the EFS diagram (i.e., conservation of the momentum
parallel to the interface), the group velocities of the three points intersecting the surface modes will
give the propagation directions of the surface waves [Fig. \ref{figNR}(b)]. Consequently, as shown in Fig. \ref{figNR}(c),
the incident (wave 1) and refracted (wave 3) surface waves lie on the same side of the normal line, indicating a negative refraction.
The existence of negative refraction of the topological surface states would lead to many interesting
effects\cite{Veselago68,Pendry00,Chan18,He18,Fang05,Liu07,Smolyaninov07,Lin08}
such as superlens effect\cite{Pendry00,Fang05,Liu07,Smolyaninov07} and suppression of reflected waves \cite{He18}.

\begin{figure
\includegraphics[width=3.2in]{figure6.pdf}
\caption{\label{figNR}
Schematic illustration of the negative refraction of the TPSS of the GHM.
(a) Side view and (b) the EFS of the system. (c) The incident,
reflected, and transmitted waves at the interface between the two different GHM insulators.
whose group velocities are denoted as  1, 2, and 3 in (b), respectively.
Here $v_{vg,x}$ and $v_{ph,x}$ denote the $x$ components of the group and phase velocities, respectively.}
\end{figure}

\section{Realization of the proposed gyromagnetic hyperbolic metamaterials}

The hyperbolic metamaterials are a highly anisotropic material with real parts of the principal components
of its permittivity tensor having opposite signs \cite{HMM13} and have been intensively investigated
because of this unique property. Here we consider a bilayer superlattice composed of a metal layer with
relative permittivity $\epsilon_m$ and thickness $d_m$ and a dielectric layer with relative
permittivity$\epsilon_d$ and thickness $d_d$ stacked along the $z$-direction. The effective  dielectric constant
of this metamaterial is given by
\begin{equation}
\epsilon_{eff}=\left[
\begin{array}{ccc}
\epsilon_{x-y} & 0 & 0 \\
0 & \epsilon_{x-y} & 0 \\
0 & 0 & \epsilon_{z} \\
\end{array}
\right],
\end{equation}
where
\begin{subequations}
\begin{equation} \epsilon_{x-y}=\frac{\epsilon_md_m+\epsilon_dd_d}{d_m+d_d}, \end{equation}
\begin{equation} \frac{1}{\epsilon_{z}}=\frac{d_m/\epsilon_m+d_d/\epsilon_d}{d_m+d_d}. \end{equation}
\end{subequations}

If the dielectric layer is replaced by a gyromagnetic medium layer, the superlattice becomes a gyromagnetic
hyperbolic metamaterial with its relative permeability tensors
$\begin{bmatrix}
\mu_{g} & -i\gamma & 0 \\
i\gamma & \mu_{g} & 0 \\
0 & 0 & \mu_{g} \\
\end{bmatrix}$
and
$\begin{bmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 \\
\end{bmatrix}$
for the gyromagnetic and metal layers, respectively. The bilayer superlattice with relative
permittivity $\epsilon_g= \epsilon_d$ and thickness $d_g= d_d$ for the gyromagnetic layer, is schematically shown in Fig. 7.

\begin{figure
\includegraphics[width=3.2in]{figure7.pdf}
\caption{Proposed gyromagnetic hyperbolic metamaterial as a bilayer
superlattice composed of a gyromagnetic slab and a metal slab.}
\end{figure}

For simplicity, let us write the fields in the metal as $(\mathbf{E}^g, \mathbf{H}^g)=(E_x^g, E_y^g, E_z^g, H_x^g, H_y^g, H_z^g)$
and in the gyromagnetic medium as $(\mathbf{E}^m, \mathbf{H}^m)=(E_x^m, E_y^m, E_z^m, H_x^m, H_y^m, H_z^m)$ .
We assume that the wavelength is much larger than the thickness. Using the boundary conditions
that $D_z$, $E_x$, $E_y$, $B_z$, $H_x$ and $H_y$ must be continuous across an interface, we obtain
$\epsilon_mE_z^m=\epsilon_gE_z^g$, $E_x^m=E_x^g=E_x$, $E_y^m=E_y^g=E_y$, $H_z^m=\mu_gH_z^g$,
$H_x^m=H_x^g=H_x$, $H_y^m=H_y^g=H_y$.
Consequently, the average $D$-field and $B$-field are
\begin{subequations}
\begin{equation} D_x^{eff}=\frac{d_g\epsilon_g+d_m\epsilon_m}{d_g+d_m}E_x   \end{equation}
\begin{equation} D_y^{eff}=\frac{d_g\epsilon_g+d_m\epsilon_m}{d_g+d_m}E_y   \end{equation}
\begin{equation} D_z^{eff}=\epsilon_mE_z^m=\epsilon_gE_z^m                  \end{equation}
\begin{equation} B_x^{eff}=\frac{d_g \mu_g+d_m}{d_g+d_m}H_x -\frac{id_g\gamma}{d_g+d_m}H_y \end{equation}
\begin{equation} B_y^{eff}=\frac{d_g \mu_g+d_m}{d_g+d_m}H_y +\frac{id_g\gamma}{d_g+d_m}H_x  \end{equation}
\begin{equation} B_z^{eff}=\mu_0H_z^m=\mu_gH_z^m. \end{equation}
\end{subequations}
Therefore, the effective relative permittivity and permeability tensors are given by
\begin{subequations}
\begin{equation}
\epsilon_{eff}=\left[
\begin{array}{ccc}
\frac{d_g\epsilon_g+d_m\epsilon_m}{d_g+d_m} & 0 & 0 \\
0 & \frac{d_g\epsilon_g+d_m\epsilon_m}{d_g+d_m} & 0 \\
0 & 0 & \frac{(d_g+d_m)\epsilon_g\epsilon_m}{d_g\epsilon_m+d_m\epsilon_g} \\
\end{array}
\right],
\end{equation}

\begin{equation}
\mu_{eff}=\left[
\begin{array}{ccc}
\frac{d_g\mu_g+d_m}{d_g+d_m} & \frac{-id_g\gamma}{d_g+d_m} & 0 \\
\frac{id_g\gamma}{d_g+d_m} & \frac{d_g\mu_g+d_m}{d_g+d_m} & 0 \\
0 & 0 & \frac{(d_g+d_m)\mu_g}{d_g\epsilon_m+d_m\mu_g} \\
\end{array}
\right].
\end{equation}

\end{subequations}

We can easily find suitable materials as the gyromagnetic medium and the metal slab to
construct a GHM described above. For example, yttrium iron garnet (YIG) under an applied magnetic
field of 1600 Gauss has the effective relative permittivity $\epsilon_g=15$ and permeability $\mu_g=1.12$
and $\gamma=0.124$ at 1.94 THz~\cite{QAH_Wangzhen_theory}, and thus can be used as the gyromagnetic medium. InSb has the
effective relative permittivity $\epsilon_m$ = -10.78 at 1.94 THz~\cite{Wu2014} and thus can serve as the metal slab.
In the bilayer multilayer, the thicknesses of the gyromagnetic medium and metal layers are taken to be
the same and much less than 1.546$\times$ $10^{-4}$ m (wavelength of 1.94 THz). In this case,
the effective relative permittivity and permeability tensors can be written as

\begin{subequations}
\begin{equation}
\epsilon_{eff}=\left[
\begin{array}{ccc}
2.11 & 0 & 0 \\
0 & 2.11 & 0 \\
0 & 0 & -76.73 \\
\end{array}
\right],
\end{equation}

\begin{equation}
\mu_{eff}=\left[
\begin{array}{ccc}
1.06 & 0.062i & 0 \\
-0.062i & 1.06 & 0 \\
0 & 0 & 1.06 \\
\end{array}
\right].
\end{equation}
\end{subequations}

The EFS dispersions obtained using the effective parameters are also displayed in Fig. 8. Clearly,
the shape of the EFS dispersions is similar to that of the EFS presented in the last section,
and furthermore, the topological properties such as the Chern numbers of the bands and the topological nature
of the band gaps are identical to those shown in the last section.

\begin{figure}[htbp]
\includegraphics[width=3.2in]{figure8.pdf}
\caption{Calculated EFS of the bilayer superlattice shown in Fig. 7
made of a YIG slab as the gyromagnetic medium and an InSb layer as the metal slab.}
\end{figure}

\begin{figure}[htbp]
\includegraphics[width=3.2in]{figure9.pdf}
\caption{Calculated band gap map on the $\gamma-(\epsilon_{zz}/\epsilon_{xx})$ plane of the bilayer superlattice
shown in Fig. 7, which is made of a YIG slab as the gyromagnetic medium and an InSb layer as the metal slab.}
\end{figure}

To see the robustness of the topological band gap against the variations of the material parameters,
we show the calculated band gap as a function of both $\gamma$ and $\epsilon_{zz}/\epsilon_{xx}$ in Fig. 9.
Clearly, the band gap is significant $(>0.5)$ in a wide range of the parameters. For example,
when $\epsilon_{zz}/\epsilon_{xx}$ is in between -0.1 and -1.0, the band gap is significant for $\gamma$
being from 0.4 to 2.0. Interestingly, Fig. 9 shows that for $\epsilon_{zz}/\epsilon_{xx}$ is in
between -1.0 and -2.0, one tunes the system from the gapless to the Chern insulating state
and then back to gapless state.
Finally, it should be pointed out that the present photonic Chern insulator would be much simpler
to design and fabricate than many other photonic insulators. For example, the photonic Floquet
topological insulators reported in Ref. [\onlinecite{floquetrechtsman2013}] are composed of an array of
micrometer-scale helical waveguides. The photonic chiral hyperbolic topological insulators proposed
in Ref. [\onlinecite{chiralhypermeta}] would be based on an array of hyperbolic coils which serve as chiral resonators.
In contrast, as discussed above, the photonic Chern insulators would consist of a simple YIG/InSb bilayer
superlattice, and thus would be much easier to fabricate.

\section{Discussion and conclusions}

In the QAH effect, vacuum is an insulator for electrons and thus the chiral edge states are localized
at the surface of the Chern insulator. In contrast, vacuum to photons is often like a free electron
metal and thus light propagation on the unidirectional edge modes would generally
suffer from radiation loss. Consequently, a metal film such as copper is usually inserted between air
and the topological photonic insulator to suppress the radiation leakage\cite{QAH_Wangzhen_exp,QAH_poo2013,QAH_largeChern_exp}.
However, metals like copper could incur significant Ohmic loss. In strong contrast, the chiral edge modes
in our GMHs are not only reflection-free but also non-radiative, as demonstrated
by our finite element electromagnetic simulations (see Figs. 4 and 5). This is simply
because the chiral edge dispersions are located outside light cone
[see Fig. 4(b)]. This is an important advantage of the GMHs
over the gyromagnetic photonic crystals\cite{QAH_Wangzhen_exp,QAH_poo2013}.

In conclusion, by both theoretical analysis and electromagnetic simulations,
we have demonstrated that gyromagnetic hyperbolic metamaterials (GHM)
are photonic Chern insulators with fascinating properties.
We further show that the large topogical band gaps with gap Chern number of one in these metamaterials,
result from the simultaneous presence of the hyperbolicity and also the gyromagnetic effect,
which breaks the time-reversal symmetry and thus gives rise to nonzero Berry curvatures on the EFSs.
Remarkably, unlike many other photonic Chern insulators, the GHM Chern insulators possess
non-radiative chiral edge modes on their surfaces, and thus allow to fabricate unidirectional waveguides
without cladding metals which generally incurr considerable Ohmic loss.
Furthermore, the photonic edge states in the proposed Chern insulators are robust against disorder
on a wide range of length scales, in strong contrast to crystalline topological insulators,
and the light flow direction on the surface of the Chern insulators can be easily
flipped by switching the direction of an applied magnetic field.
We also uncover negative refraction of the topological surface wave at the boundary between
the GHMs with the opposite signs of gyromagnetic parameters. Finally, we show that
compared with other photonic topological materials such as chiral hyperbolic materials\cite{chiralhypermeta}, the present
GHM Chern insulators can be much easier to fabricate.

\begin{acknowledgments}
This work is supported by the Ministry of Science and Technology, the National Center for Theoretical Sciences
and the Thematic Research Program (AS-TP-106-M07) of the Academia Sinica of the R. O. C. (Taiwan).
\end{acknowledgments}


\section{Introduction}

  Ultracold Bosons in optical lattices (OLs) \cite{bloch_08} have been the
subject of intense research since the experimental realization of Bose
Einstein condensate (BEC) in OLs
\cite{anderson_98, greiner_01, greiner_02,lewenstein_07}.
The rapid developments in the control of BECs in OLs have made
these systems an elegant experimental tool to explore the foundations of
strongly correlated quantum many-body systems. Till recently, many of these
were limited to the realm of theoretical studies. The near ideal, defect free,
experimental realizations of OLs make these excellent proxies to explore
quantum many-body effects in condensed matter systems. One remarkable recent
experimental development is the introduction of synthetic magnetic fields in
OLs\cite{aidelsburger_13, miyake_13}. This makes the physics of
quantum Hall (QH) effect accessible to OLs. An important outcome is the
study of Harper-Hofstadter model \cite{harper_55,hofstadter_76} and observation
of fractal spectrum \cite{dean_13} for interacting Bosons in OL with synthetic
magnetic field \cite{jaksch_03}. The scope to investigate the interplay of
lattice geometry, synthetic magnetic fields, and strong interactions have made
these systems an excellent platform to explore exotic quantum many-body phases.

In the quantum description, the energies of electrons in an external magnetic
field are quantized in to Landau levels. These have large degeneracy and in
a lattice these correspond to the Bloch bands \cite{harper_55} and is
sensitive to applied magnetic field. The key point is the geometrical phase
an electron acquires when completing a loop in the cyclotron motion. Thus,
neutral atoms in OLs can mimic the physics of electrons in magnetic fields
if a geometrical phase can be induced to the atoms. This is achieved through
the generation of a synthetic magnetic field in OLs through an artificial gauge
potential \cite{dalibard_11, lin_11, lin_09a, garcia_12} using lasers. Then,
an atom hopping around a single unit cell in the OL, also called a plaquette,
acquires Peierls' phase \cite{peierls_33} of $\Phi = 2\pi\alpha$. Where
$\alpha$ is the flux quanta per plaquette and it is related to the strength of
the synthetic magnetic field. In the condensed matter systems realizing
high $\alpha$ require magnetic fields $\approx 10^3$ Tesla.
In this respect, the OLs have the advantage that by suitable choice of
external as well as internal parameters, various topological states such as
fractional QH (FQH) states are obtainable with the current experimental
realizations \cite{palmer_06, sorensen_05} and a variety of FQH-like states
can be expected to emerge from these systems.

A paradigmatic model which describes BECs in OLs is the Bose Hubbard
model (BHM) \cite{fisher_89,jaksch_98}. In this model the kinetic energy of
the bosons competes with the on site interaction and drives a quantum phase
transition (QPT) from  superfluid (SF) to bosonic Mott insulator
(MI) phase \cite{greiner_02, stoferle_04}. Various theoretical methods, such
as mean field theory \cite{fisher_89}, strong coupling expansion
\cite{freericks_96, niemeyer_99, freericks_09, wang_18}, quantum Monte Carlo
\cite{wessel_04}, density matrix renormalization group \cite{peotta_14} have
been used to study the role of quantum fluctuation and short range on site
interaction on QPT. The SF phase is compressible with finite SF order parameter
and phase coherent; MI phase on the other hand is incompressible
with zero order parameter and shows integer commensurate filling per
lattice site. In contrast to these two phases, the QH states are
incompressible states with zero order parameter and have incommensurate
filling. Several previous works
\cite{sorensen_05,palmer_06,hafezi_07,umucalilar_07, palmer_08,umucalilar_10,
natu_16, hugel_17,kuno_17,gerster_17,bai_18} have theoretically explored the
existence of FQH states in OLs using BHM with synthetic magnetic fields, the
bosonic counterpart of the Harper-Hofstadter model
\cite{harper_55,hofstadter_76}. These theoretical works have also examined the
possible signatures of the FQH states. One of the possibilities is the
measurement of two-point correlation function in the bulk and in the edge of
the lattice \cite{he_17}. Such measurements may be experimentally possible
using the concepts from quantum information theory \cite{elliott_16,streif_16}.

In the present work we use cluster Gutzwiller mean field (CGMF) theory to
provide a better description of the atom-atom correlations. It is proven to
be more accurate than single site Gutzwiller theory and previous works
have also used CGMF theory to study both the integer quantum Hall (IQH) and
FQH states. In ref \cite{natu_16} the appearance of incompressible QH ground
state in the hard-core limit with stripe order for $\alpha$ = 1/5 and
$\nu = 1/2$ is reported. Similarly, using reciprocal cluster mean field (RCMF)
analysis, a competing FQH state is reported for $\alpha = 1/4$ in the recent
study by H$\ddot{o}$gel {\it et al.} \cite{hugel_17}. Motivated by the above
observation of QH states, we explore the possible QH states for $\alpha = 1/3$
with different cluster size in the hard core limit and demonstrate the
dependence of the QH state geometry on the cluster size. For the present
studies we have considered $3 \times 2$ and $3 \times 3$ clusters and
report the improvement in the description of the QH states with the
increase of cluster sizes. This stems from more accurate accounting of the
correlation effect with larger cluster sizes. We have also performed a
comparative study between the obtained QH and SF states with both the cluster
sizes and detect the emergence of various patterns.

\section{Theory}

We study a system of spinless bosonic atoms at T$ = 0$K, confined
in a two-dimensional (2D) optical lattice of square geometry under the
influence of artificial gauge field
\cite{jaksch_03,aidelsburger_11,aidelsburger_13,miyake_13}.
In the Landau gauge $\mathbf{A} = (A_x, 0, 0)$ with $A_x = 2\pi\alpha q$, the
system is described by the following Hamiltonian
~\cite{jaksch_98,jaksch_03,sorensen_05,palmer_06,palmer_08},
BHM Hamiltonian where Peierls substitution is incorporated in the
nearest-neighbor (NN) hopping ~\cite{peierls_33,harper_55,hofstadter_76},
\begin{equation}
 \hat{H} = -\sum_{\langle jk\rangle}\left({\rm e}^{i2\pi\alpha q}J_x
           + J_y \right ) \hat{b}_j^{\dagger}\hat{b}_k
           + \sum_j\hat{n}_j\left[\frac{U}{2}(\hat{n}_j-1)
           - \mu )\right],
\label{bhm}
\end{equation}
$j \equiv (p,q)$ corresponds to the lattice site index where $q$ is the index
of the lattice site along $y$-axis, $\hat{b}_{j}$ ( $\hat{b}^{\dagger}_{j}$)
are the bosonic annihilation (creation) operators, $\hat{n}_{j}$ is the
occupation number operator at $j^{th}$ lattice site, $J_x$ ( $J_y$) are the
hopping strengths between two neighbouring sites along $x$ ($y$) direction,
$U$ corresponds to the on-site interaction and $\mu$ is the chemical potential.
Based on the experimental realizations, we consider isotropic hopping
$J_x = J_y = J$, and repulsive on-site interaction energy ($U>0$).
In presence of synthetic magnetic field, the atoms acquire a
phase $2\pi\alpha$ upon hopping around a plaquette, where, $\alpha$ is
the number of flux quanta per plaquette, and it has values
$0 \leq\alpha\leq 1/2$. In the absence of synthetic magnetic field
($\alpha$ = 0), the Hamiltonian (\ref{bhm}) reduces to the familiar
BHM Hamiltonian which admits two possible phases ---
MI and SF~\cite{fisher_89,jaksch_98,greiner_02}. The MI phase appears in the
strongly interacting regime $(J/U\ll 1)$, while SF phase occurs in the
limit $(J/U\gg 1)$. In homogeneous system, where the OL does not include any
background potential, the phase-boundary between MI and SF forms lobes of
different fillings and in presence of magnetic field the MI-lobes
are enhanced \cite{oktel_07}. We employ the single site Gutzwiller mean-field
method (SGMF) and CGMF to analyze the system in presence of synthetic magnetic
field and obtain the QH states.

\subsection{Mean Field Theory And Gutzwiller Approximation}

  Following the mean-field \cite{sheshadri_93} calculations of
the BHM, we decompose the creation and annihilation operators
in Eq.~(\ref{bhm}) into mean field and fluctuation around the
mean-field, that is, $\hat{b}_{j} = \phi_{j} + \delta \hat{b}_{j}$,
with $\phi_{j} = \langle\hat{b}_{j}\rangle$, and similarly,
$\hat{b}^{\dagger}_{j} = \phi^{*}_{j} + \delta \hat{b}^{\dagger}_{j}$,
with $\phi^{*}_{j} = \langle\hat{b}^{\dagger}_{j}\rangle$. Neglecting the
term quadratic in the fluctuations, the mean-field Hamiltonian is
\begin{eqnarray}
\hat{H}^{\rm MF} &=&- \sum_{\langle jk\rangle}
                  \left[(J_x {\rm e}^{i2\pi\alpha q} + J_y)
                  \left(\hat{b}_{j}^{\dagger}\phi_{k}
                  + \phi^{*}_{j}\hat{b}_{k}
                  -\phi^{*}_{j}\phi_{k}\right)
                  \right. \nonumber \\
                  && + \left.  {\rm h.c.}\right]
                  + \sum_{j} \left[\frac{U}{2}\hat{n}_{j}
                  (\hat{n}_{j}-1)
                  - \mu\hat{n}_{j}\right].
\label{mf_hamil}
\end{eqnarray}
We can, therefore, express the total Hamiltonian as the sum of single site
mean field Hamiltonians
\begin{eqnarray}
\hat{h}_{j} =&-&\left [(J_x{\rm e}^{i2\pi\alpha q} + J_y)
                  \left(\phi_{k}^*\hat{b}_{j} - \phi^{*}_{k}\phi_{j} \right)
                  + {\rm H.c.}\right ]\\ \nonumber
                  &+& \frac{U}{2}\hat{n}_{j}(\hat{n}_{j}-1)
                  - \mu \hat{n}_{j}.
   \label{ham_ss}
\end{eqnarray}
The next step is to diagonalize the Hamiltonian (\ref{mf_hamil}) for each
site separately. For this, we consider the Gutzwiller ansatz, that is the
ground  state of the entire lattice is the direct product of the ground
states of all the individual sites, and can be written in the Fock basis as
\begin{eqnarray}
 |\Psi_{\rm GW}\rangle = \prod_j|\psi\rangle_j
                       = \prod_j \sum_{n = 0}^{N_{\rm b}}c^{(j)}_n
                         |n\rangle_j,
 \label{gw_state}
 \end{eqnarray}
with the normalization condition $\sum_{n} |c^{(j)}_n|^2$ = 1.
Here, $N_b$ is the occupation number state maximum number of particles at a
site and  $c^{(j)}_n$ corresponds the complex co-efficients for the
ground state $|\psi\rangle_{j}$ at the $j$th site. Gutzwiller ansatz is
the exact solution of the system in the strongly interacting regime ($J\ll U$).

For the numerical computations, we consider $N_{\rm b} = 10$ and choose an
initial guess of $\phi$. Then, we diagonalize the Hamiltonian for each site
and retain the ground state as the state $\ket{\psi}_j$ in
$|\Psi_{\rm GW}\rangle$. Then, using this $|\psi\rangle_j$, we
calculate new $\phi$ for the next iteration and this cycle is
continued till convergence is reached. To distinguish the different phases, we
compute the SF order parameter at each site, and  for the $j$th lattice site
SF order parameter is
\begin{equation}
\phi_j = \langle\Psi_{\rm GW}|\hat{b}_j|\Psi_{\rm GW}\rangle
            = \sum_{n = 1}^{N_{\rm b}}\sqrt{n}
              {c^{*(j)}_{n-1}}c^{(j)}_{n}.
\label{gw_phi}
\end{equation}
From the above expression it is evident that $\phi_{j}$ is zero in the MI
phase of the system since only one of the co-efficients in Eq.~(\ref{gw_state})
is non-zero, and it is finite for the SF phase due to the different $c_n$s
contribution. We also compute the average lattice occupancy or density at
each of the lattice site as
\begin{equation}
\rho_j = \langle\Psi_{\rm GW}|\hat{n}_j|\Psi_{\rm GW}\rangle
            = \sum_{n = 0}^{N_{\rm b}}n |c^{(j)}_n|^2.
\label{aver_n}
\end{equation}
These are the essence of SGMF theory.

\subsection{Theory of CGMF}

In the SGMF Hamiltonian Eq.~(\ref{ham_ss}), we decouple the hopping terms
between two neighbouring sites by considering the mean field or SF order
parameter $\phi$. Thus, we could write the Hamiltonian of the entire
system as the sum of Hamiltonians of individual sites and implement it as a
site wise computations.  However, this approximation is inadequate
to incorporate the correlation effects arising from the NN hopping. To remedy
this short coming, which assumes great importance to describe strongly
correlated states like QH states, previous works have relied on
CGMF \cite{natu_16}. To derive the CGMF Hamiltonian, we consider the entire
lattice size as $K \times L$, which we divide into $W$ clusters ($C$) of
size $M \times N$, i.e., $W = (K\times L)/(M\times N)$. The case of
$M = N = 1$ corresponds to the SGMF theory. In CGMF Hamiltonian, the hopping
term is decomposed into two parts. First is the actual hopping term in the
internal link of the cluster($\delta$C) and the second term takes care of the
boundary via mean fields. Our recent study \cite{bai_18} describes the
decomposition of hopping term and CGMF method more clearly. After
decomposition, the Hamiltonian for a single cluster($C$) is expressed in
the following way
\begin{eqnarray}
 \hat{H}_C  &=& -\sum_{p, q \in C}\left[\left({\rm e}^{i2\pi\alpha q}J_x
              \hat{b}_{p + 1, q}^{\dagger}\hat{b}_{p, q} + {\rm H.c.}\right)
               \right. \nonumber\\
              &&+ \left. \left(J_y \hat{b}_{p, q - 1}^{\dagger}\hat{b}_{p, q}
              + {\rm H.c.}\right)\right]
              \nonumber\\
              &&-\sum_{p, q\in \delta C}
              \left[\left({\rm e}^{i2\pi\alpha q}J_x
              \langle a_{p,q}^{*} \rangle \hat{b}_{p, q} + {\rm H.c.}\right)
              \right. \nonumber \\
              &&+ \left. \left(J_y \langle a_{p,q}^{*} \rangle\hat{b}_{p, q}
              + {\rm H.c.}\right)\right]
              \nonumber\\
              &&+\sum_{p, q \in C}
              \left[\frac{U}{2}\hat{n}_{p, q}(\hat{n}_{p, q}-1) -
              \mu\hat{n}_{p, q}\right],
\label{cg_hamil}
\end{eqnarray}
where $\langle a_{p,q} \rangle  = \sum_{p^{'},q^{'}\in \not C} \langle
b_{p^{'},q^{'}}\rangle$. Then, we use Gutzwiller ansatz and the local
cluster wavefunction in a Fock basis can be expressed as
\begin{equation}
\ket{\psi_c} = \sum_{n_1,n_2,\ldots,n_{MN}}
  C_{n_1,n_2..,n_{MN}}\ket{n_1,n_2,\ldots,n_{MN}},
\end{equation}
with $n_i$ being the index of the occupation number state of $i$th lattice
site within the cluster, and $C_{n_1,n_2,\ldots,n_{MN}}$ is the amplitude
of the cluster Fock state $|n_1,n_2,\ldots, n_{MN}\rangle$.
Here also, the total Hamiltonian of the system can be written as the sum of
all the individual cluster Hamiltonians \cite{luhmann_13}.
The SF order parameter $\phi$ is computed for each cluster in the similar way
as discussed in SGMF method. The next step is to find the ground state and we
adopt the similar process as is described in SGMF theory.
We take the initial solution for $\phi$, construct the Hamiltonian matrix
elements for a single cluster and diagonalize it. After diagonalization,
we consider the lowest ground state $\ket{\psi_c}$ for the cluster and
calculate the new $\phi$ and repeat the cycle until we get the converged
solution. Here, it is worth to be mention that in case of CGMF, the
convergence is very sensitive to the initial
conditions \cite{natu_16,kuno_17} and we use the method of
successive over-relaxation for the better convergence \cite{barrett_94}.

\begin{figure}[t]
  \includegraphics[width=8.3cm]{compr_alq_1b3_3b2and3b3_cluster.pdf}
  \caption{The variation in the number density $\rho$ for $\alpha = 1/3$
           as function of $\mu$. The states in SF phase are compressible and
           have non-zero superfluid order parameter $\phi$. As a result,
           $\rho$ varies linearly with $\mu$ and the green curve represents
           the SF states. For specific values of filling factor $\nu$ there are
           states with constant $\rho$, represented by the blue lines, and
           these correspond to the QH states.
           (a) Results from $3\times2$ cluster, the plateaus or
           the constant $\rho$ values correspond to $\nu = n/2, n =1,2,..,5$
           and the corresponding $\rho$ values are $n/6$.
           (b) Results from $3\times3$, the platues correspond to
           $\nu = n/3, n = 1,2,..,8$ and the corresponding $\rho$ values
           are $n/9$.}
  \label{kappa_alpha1b3}
\end{figure}
\begin{figure}[t]
  \includegraphics[width=8.3cm]{hallstate_3b2_3b3_nu1_comb.pdf}
  \caption{The IQH state for $\alpha = 1/3$ and $\nu = 1$ with (a) $3\times2$
           cluster and (b) with $3\times3$ cluster. The IQH state switches from
           stripe to checkerboard geometry with the mentioned cluster sizes.}
  \label{hs_nu_1_3b2_3b3}
\end{figure}

\section{Results and Discussions}

We start our computations by considering a single cluster for
$\alpha = 1/3$, with $3\times2$ and $3\times3$ clusters. We choose
these clusters as Hamiltonian becomes periodic over a $3 \times 1$
magnetic unit cell in the Landau gauge for this flux value of $\alpha$.
We obtain QH states from the CGMF and characterize them based on the
compressibility $\kappa = \partial \rho /\partial \mu$, where the
density for the cluster is
$\rho = \sum_j\bra{\psi_c} \hat{n}_j \ket{\psi_c}/(K\times L)$.
As the QH states are incompressible $\kappa=0$ for these states,
and $\kappa$ is finite for the compressible SF states.
Therefore, $\rho(\mu)$ of QH states has plateaus for different fillings
$\nu$ and $\rho(\mu)$ is linear for the SF phase. In the
Fig.~\ref{kappa_alpha1b3}, the plateaus corresponding to constant $\rho$
indicate the existence for the QH states. Our computations, as mentioned
earlier, are in the hard-core boson limit where $\rho < 1$. We obtain the QH
states at $\nu = n/2$, with $n =1,2,..,5$ by taking the $3\times2$ cluster,
and at $\nu = n/3$, with $n = 1,2,..,8$ by taking the $3\times3$ cluster. The
QH states are enhanced with the larger cluster size as mentioned above.
Here, $3\times3$ cluster is close to exact diagonalization (ED) as the central
lattice site has exact hopping contributions from the nearest neighbor sites.
And, indeed, the diagonalization of the cluster can be transformed into ED with
minor modifications in the computations of the Hamiltonian matrix elements. One
main reason for enhancement in the QH states with $3\times3$ cluster is that
it describes correlations effects more accurately compared to the $3\times2$
cluster and hence, the results are more accurate. Further, we show the density
plots for the QH and SF states for the larger lattice system. For this, we take
$12 \times 12$ lattice sites and $J/U = 0.01$. This system size,
and hopping energies are kept the same for all the QH and SF states discussed
in the rest of the manuscript. The IQH state at filling $\nu = 1$,
$\mu/U = -0.008$ with density $\rho = 1/3$ is shown in
Fig.~\ref{hs_nu_1_3b2_3b3}. The IQH state has stripe pattern with $3\times2$
cluster and transforms into checkerboard pattern with $3\times3$ cluster. The
transformation from the stripe to checkerboard pattern is observed for the
other IQH state of $\nu = 2$ as well. We observe that all the
IQH and FQH states have stripe pattern except the $\nu = 3/2$ state, which
has homogeneous density with $\rho = 0.5$ with $3\times2$ cluster. And the
corresponding SF states have a zigzag pattern in the density $\rho$ and
in the SF order parameter $\phi$. One of the FQH state for $\nu = 5/2$ and
corresponding SF state is shown in the Fig.~\ref{nu_5b2_3b2}.
\begin{figure}[t]
  \includegraphics[width=8.3cm]{nu_5b2_3b2_hs_sf.pdf}
  \caption{(a) The FQH state with $\alpha = 1/3$
           and $\nu = 5/2$ with $3\times2$ cluster, it has a stripe pattern
           in the density with vanishing SF order parameter.
           (b) The analogous SF state with zigzag pattern in
           the density as well as in the SF order parameter.}
  \label{nu_5b2_3b2}
\end{figure}
\begin{figure}[t]
  \includegraphics[width=8.3cm]{nu_1b3_3b3_hs_sf.pdf}
  \caption{(a) The represents FQH state with $\alpha = 1/3$
           and $\nu = 1/3$ with $3\times3$ cluster with a checkerboard pattern
           in the density and vanishing SF order parameter.
           (b) The analogous SF state with diagonal stripe pattern
           in the density as well as in the SF order parameter.}
  \label{nu_1b3_3b3}
\end{figure}
As discussed earlier, we do not observe the half integer FQH states with
$3\times3$ cluster, but do observe the FQH states at the one third fillings.
One of the FQH and SF state with $3\times3$ cluster for $\nu = 1/3$ is
shown in the Fig.~\ref{nu_1b3_3b3}. Here, with $3\times3$ cluster, we observe
all the FQH states have checkerboard pattern and all the SF states have the
diagonal stripe pattern. We find that in all the cases SF states is the
ground state and QH state is metastable state.

\section{Conclusion}

We obtain QH states by considering the two cluster sizes as $3 \times 2$
and $3 \times 3$ in the CGMF theory. With the larger cluster $3 \times 3$, we
approach ED as the CGMF provides exact description of the hopping term for
the central lattice site. We obtain QH states with fillings
$\nu = n/2$, $n = 1,2,..,5$ and corresponding density $\rho = n/6$ with
$3\times 2$ cluster size. However, with $3\times 3$ cluster, we obtain
a larger set of QH states with fillings $\nu = n/3$, $n = 1,2,..,8$
and corresponding density $\rho = n/9$. We also observe the competing SF
states corresponding to all the QH states. We find that the SF state is
the ground state and QH state is metastable state in all the cases.
We have demonstrated that the QH states change geometry from stripe to
checkerboard by switching $3\times 2$ to $3\times 3$ cluster size. On the
other hand, SF state has zigzag pattern with both the cluster sizes, specially
with $3\times3$ the zigzag pattern is equivalent to diagonal stripe pattern.
Thus, we have established that to obtain correct density pattern of the
QH states, it is essential to consider larger cluster sizes in the CGMF
theory.

\begin{acknowledgments}
The results presented in the paper are based on the computations
using Vikram-100, the 100TFLOP HPC Cluster at Physical Research Laboratory,
Ahmedabad, India. We thank Arko Roy, S. Gautam and S. A. Silotri for valuable
discussions.
\end{acknowledgments}