\section{Introduction}
Recently, an exactly solvable many-body fermionic model with an infinite-range interaction, i.e. the Hatsugai-Kohmoto (HK) model,\cite{Hatsugai1992,Baskaran1991,Hatsugai1996} and its extension has been hotly studied.\cite{Phillips2018,Yeo2019,Phillips2020,Yang2021,Zhu2021,Zhao2022,Setty2021,Mai2022,Huang2022,Li2022,Setty2020,Setty2021b,Zhong2022}
In contrast with celebrated and more familiar Sachdev-Ye-Kitaev, Kitaev's toric code and honeycomb model with either quenched disorder or local $Z_{2}$ gauge symmetry,\cite{Sachdev,Maldacena,Chowdhury,Kitaev1,Kitaev2,Prosko,Zhong2013,Smith2017} the original HK model has translation invariance with topologically trivial nature, but surprisingly, it provides a strictly exact playground for non-Fermi liquid (NFL) and featureless Mott insulator in any spatial dimension, which is rare in statistical mechanics and condensed matter physics. (See Fig.~\ref{fig:fig_0})


The solvability of HK model results from its locality in momentum space and one can diagonalize HK Hamiltonian (just diagonal $4\times4$-matrix) for each momentum. The current studies have mainly focused on an interesting extension of HK model, i.e. the superconducting instability from the intrinsic NFL state in HK model,\cite{Phillips2020} which is inspired by ubiquitous NFL behaviors and their link to unconventional superconductivity in cuprate, iron-based superconductors (SC) and many heavy fermion compounds. Unexpected properties such as topological $s$-wave pairing and two-stage superconductivity have been discovered.\cite{Zhu2021,Zhao2022,Li2022}
\begin{figure}
\includegraphics[width=0.75\linewidth]{fig_0.pdf}
\caption{\label{fig:fig_0} (a) The Hatsugai-Kohmoto (HK) model with hopping strength $t$ and interaction $U$ in one spatial dimension and (b) on a square lattice. (c) The exact ground-state phase diagram for HK model exhibits a
Mott insulator and a non-Fermi-liquid-like metal. ($U_{c}=W$ with $W$ being the band-width) The transition from metallic state to
gapped Mott insulating phase belongs to the universality of the continuous Lifshitz transition.}
\end{figure}


At the same time, it is well-known that observable (like magnetization and resistivity) in metals show periodic oscillation under external magnetic field.\cite{Shoenberg1984} Such phenomena is called quantum oscillation and is generally believed to result from the oscillation of electron's density of state at Fermi energy when Landau level crosses Fermi energy periodically. The underlying microscopic description is captured by the famous Lifshitz-Kosevich (LK) formula, which provides a direct tool to extract effective mass and Fermi surface from period and amplitude of oscillation. Furthermore, LK formula is still valid if electron-electron interaction effect does not alter the Fermi liquid nature of the considered system.\cite{Luttinger1961,Wasserman1996}


However, a natural question arises when the system is not described by Fermi liquid, and we may ask how about its quantum oscillation and the fate of LK formula. This question is partially motivated by recent experiments of quantum oscillation in many strongly correlated electron systems, including cuprate and iron-based superconductor,\cite{Sebastian2015,Carrington2011} heavy fermion compounds,\cite{McCollam2005} topological Kondo insulator,\cite{Tan2015} excitonic insulator and twisted-bilayer graphene.\cite{Wang2021,Cao2016} In literature, electron correlation effect in quantum oscillation has been widely studied by mean-field approximation,\cite{Doh1998} slave-particle theory and dynamic mean-field theory.\cite{Galitski2009,Acheche2017,Vucicevic2021} It should be emphasized that all these mentioned treatments involve either uncontrollable approximation or infinite-dimension limit, thus our understanding on quantum oscillation in interacting many-electron system is obviously not complete.


In this work, inspired by the NFL state in solvable HK model, we study its possible quantum oscillation. With the Luttinger approximation, the
NFL state indeed shows quantum oscillation and its zero-temperature behavior can be captured by LK-like formula, which reflects the existence of two Fermi surface structure of (non-Landau) quasiparticle in NFL. We also perform an exact diagonolization calculation to check the validity of Luttinger approximation, but due to the limited size of the system, although numerical results exhibit certain oscillation behavior, it is hard to extract its oscillation period and amplitude. Therefore, more work (particularly the large-scale numerical simulation) on this interesting issue is highly desirable and we expect the current study on HK model will be helpful to understand generic quantum oscillation in correlated electron materials.

\section{The Hatsugai-Kohmoto model}\label{sec1}
\subsection{Quick review of Hatsugai-Kohmoto model}
The Hatsugai-Kohmoto model we study has the following form,
\begin{eqnarray}
\hat{H}&=&-\sum_{i,j,\sigma}t_{ij}\hat{c}_{i\sigma}^{\dag}\hat{c}_{j\sigma}-\mu\sum_{j\sigma}\hat{c}_{j\sigma}^{\dag}\hat{c}_{j\sigma}\nonumber\\
&+&\frac{U}{N_{s}}\sum_{j_{1},j_{2},j_{3},j_{4}}\delta_{j_{1}+j_{3}=j_{2}+j_{4}}
\hat{c}_{j_{1}\uparrow}^{\dag}\hat{c}_{j_{2}\uparrow}\hat{c}_{j_{3}\downarrow}^{\dag}\hat{c}_{j_{4}\downarrow}\label{eq1}.
\end{eqnarray}
Here, $\hat{c}_{j\sigma}^{\dag}$ is the creation operator of conduction electron (called $c$-electron for simplicity) at site $j$. $t_{ij}$ are hopping integral between $i,j$ sites. Furthermore, chemical potential $\mu$ has been added to fix electron's density. $N_{s}$ is the number of sites. The last term of $\hat{H}$ is the HK interaction,\cite{Hatsugai1992} unlike the usual Hubbard interaction $U\sum_{j}\hat{c}_{j\uparrow}^{\dag}\hat{c}_{j\uparrow}\hat{c}_{j\downarrow}^{\dag}\hat{c}_{j\downarrow}$, HK interaction is an infinite-range interaction between four electrons but preserves the center of motion for $c$-electron due to the constraint of $\delta$ function. This interaction plays a fundamental role in solving this model as we will see later.



Importantly, Eq.~\ref{eq1} is local in momentum space after Fourier transformation and the resultant Hamiltonian reads as $\hat{H}=\sum_{k}\hat{H}_{k}$,
\begin{eqnarray}
\hat{H}_{k}&=&\sum_{\sigma}(\varepsilon_{k}-\mu)\hat{c}_{k\sigma}^{\dag}\hat{c}_{k\sigma}+U
\hat{f}_{k\uparrow}^{\dag}\hat{f}_{k\uparrow}\hat{f}_{k\downarrow}^{\dag}\hat{f}_{k\downarrow}\label{eq2},
\end{eqnarray}
where $\varepsilon_{k}$ are dispersion of electrons. It is emphasized that the locality of above Hamiltonian stems from infinite-range HK interaction preserving center of motion. In contrast, the Hubbard interaction in momentum space is rather nonlocal as $U\sum_{k,k',q}\hat{c}_{k+q\uparrow}^{\dag}\hat{c}_{k\uparrow}\hat{c}_{k'-q\downarrow}^{\dag}\hat{c}_{k'\downarrow}$, thus it cannot lead to simple formalism in our model.


Now, if we choose Fock state
\begin{equation}
|n_{1},n_{2}\rangle\equiv
(\hat{c}_{k\uparrow}^{\dag})^{n_{1}}|0\rangle(\hat{c}_{k\downarrow}^{\dag})^{n_{2}}|0\rangle\label{eq3}
\end{equation}
with $n_{i}=0,1$ as basis, $\hat{H}_{k}$ can be written as a diagonal $4\times4$ matrix, whose eigen-energy is $0,\varepsilon_{k}-\mu,\varepsilon_{k}-\mu,2(\varepsilon_{k}-\mu)+U$ and the corresponding eigen-state is $|0\rangle_{k}\equiv|00\rangle,|\sigma=\uparrow\rangle_{k}\equiv|10\rangle,|\sigma=\downarrow\rangle_{k}\equiv|01\rangle,|\uparrow\downarrow\rangle_{k}\equiv|11\rangle$, which means states are empty, single occupation with spin-up and spin-down, and double occupation.


Therefore, the many-body ground-state of $\hat{H}$ is just the direct-product state of each $\hat{H}_{k}$'s ground-state, i.e. $|\Psi_{g}\rangle=\prod_{k\in\Omega_{0}}|0\rangle_{k}\prod_{k\in\Omega_{1}}|\sigma\rangle_{k}\prod_{k\in\Omega_{2}}|\uparrow\downarrow\rangle_{k}$.($\Omega_{0},\Omega_{1},\Omega_{2}$ are the momentum range for different occupation) If $\Omega_{0}=\Omega_{2}=0$, the system is a Mott insulator otherwise, we obtain a metallic state with non-Fermi liquid properties (e.g. absence of Fermi surface, violation of Luttinger theorem). Similarly, excited states and their energy are easy to be constructed, so $\hat{H}$ (Eq.~\ref{eq1}) has been solved since all eigen-states and eigen-energy are found.


It is noted that the solvability of our model involves only locality in momentum space, therefore Eq.~\ref{eq1} is solvable for arbitrary lattice geometry, spatial dimension and electron filling. Furthermore, HK model does not rely on disorder average and large-$N$ limit, which are crucial in classic Sherrington-Kirkpatrick's spin glass model and the more recent Sachdev-Ye-Kitaev model.\cite{Sherrington1975,Sachdev,Maldacena,Chowdhury} In addition, (mean-field-based) pairing term like $\hat{c}_{k\uparrow}^{\dag}\hat{c}_{-k\downarrow}^{\dag}$ has been studied in Ref.~\cite{Phillips2020,Zhu2021,Zhao2022}.


For our purpose, it is useful to present the single-particle Green's function and some ground-state or thermodynamic quantities for HK model. For example, the single-particle Green's function can be obtained in terms of equation of motion, which is read as (See Appendix.\ref{ap_A})
\begin{eqnarray}
G_{\sigma}(k,\omega)&=&\frac{1+\frac{U\langle\hat{n}_{k\bar{\sigma}}\rangle}{\omega-(\varepsilon_{k}-\mu+U)}}{\omega-(\varepsilon_{k}-\mu)}\nonumber\\
&=&\frac{1-\langle\hat{n}_{k\bar{\sigma}}\rangle}{\omega-(\varepsilon_{k}-\mu)}+\frac{\langle\hat{n}_{k\bar{\sigma}}\rangle}{\omega-(\varepsilon_{k}-\mu+U)}\label{eq4}
\end{eqnarray}
where $\langle\hat{n}_{k\bar{\sigma}}\rangle$ is the expectation value of electron number operator $\hat{n}_{k\bar{\sigma}}=\hat{c}^{\dag}_{k\bar{\sigma}}\hat{c}_{k\bar{\sigma}}$ with spin $\bar{\sigma}=-\sigma$. The pole of above Green's function tells us that there exist two quasiparticle (it is not Landau quasiparticle since adiabatical continuity into non-interacting limit $U=0$ does not work) energy bands $\varepsilon_{k}-\mu,\varepsilon_{k}-\mu+U$.


Next, at finite-$T$, the thermodynamics of HK model is determined by
its free energy density $f$, which is related to partition function $\mathcal{Z}$ as
\begin{eqnarray}
f=-\frac{T}{N_{s}}\ln\mathcal{Z},\mathcal{Z}=\mathrm{Tr }e^{-\beta \hat{H}}=\prod_{k}\mathrm{Tr} e^{-\beta \hat{H}_{k}}=\prod_{k}f_{k}\label{eq5}
\end{eqnarray}
Here, one notes that the partition function is easy to calculate since each $k$-state contributes independently. We have defined $f_{k}=1+2z_{k}+z_{k}^{2}e^{-\beta U}$ and $z_{k}=e^{-\beta(\varepsilon_{k}-\mu)}$.Then, the typical thermodynamic quantity, i.e. the heat capacity, is calculated by standard thermodynamic relation $C_{V}=-T\frac{\partial^{2}f}{\partial T^{2}}$. In addition to $C_{V}$, one can also calculate spin susceptibility $\chi_{s}$ if inserting Zeeman energy term $\hat{H}_{h}=-B(\hat{c}_{k\uparrow}^{\dag}\hat{c}_{k\uparrow}-\hat{c}_{k\downarrow}^{\dag}\hat{c}_{k\downarrow})$ into Hamiltonian $\hat{H}_{k}$. Then, it follows that the magnetization $M=-\frac{\partial f}{\partial B}$ and $\chi_{s}=\frac{\partial M}{\partial B}=-\frac{\partial^{2}f}{\partial B^{2}}$.

At zero temperature, the free energy density reduces into the ground-state energy density, which has very simple expression,
\begin{equation}
e_{g}=\frac{1}{N_{s}}\sum_{k}[(\varepsilon_{k}-\mu)\theta(\mu-\varepsilon_{k})+(\varepsilon_{k}-\mu+U)\theta(\mu-\varepsilon_{k}-U)],\label{eq6}
\end{equation}
where $\theta(x)$ is the standard unit-step function ($\theta(x)=1$ for $x>0$ and $\theta(x)=0$ if $x<0$). Therefore, the electron density at $T=0$ is found to be
\begin{equation}
n=\frac{1}{N_{s}}\sum_{k}[\theta(\mu-\varepsilon_{k})+\theta(\mu-\varepsilon_{k}-U)].\label{eq7}
\end{equation}
which indicates two Fermi surface located at $\varepsilon_{k}=\mu$ and $\varepsilon_{k}=\mu-U$.
\subsection{Adding orbital magnetic field}
Because we are interested in possible quantum oscillation of HK model, orbital magnetic field should be included. This can be realized easily via Pierls substitution as $t_{ij}\rightarrow t_{ij}e^{ia_{ij}}$, where phase factor $a_{ij}=\frac{e}{\hbar}\int_{r_{j}}^{r_{i}}d\vec{r}'\cdot \vec{A}(\vec{r}')$ and $\vec{A}$ is the vector potential of external magnetic field. Then, we can write HK model as $\hat{H}=\hat{H}_{0}+\hat{H}_{\mu}+\hat{H}_{U}$
\begin{eqnarray}
&&\hat{H}_{0}=-\sum_{i,j,\sigma}t_{ij}e^{ia_{ij}}\hat{c}_{i\sigma}^{\dag}\hat{c}_{j\sigma},\nonumber\\
&&\hat{H}_{\mu}=-\mu\sum_{j\sigma}\hat{c}_{j\sigma}^{\dag}\hat{c}_{j\sigma},~~~~\nonumber\\
&&\hat{H}_{U}=\frac{U}{N_{s}}\sum_{j_{1},j_{2},j_{3},j_{4}}\delta_{j_{1}+j_{3}=j_{2}+j_{4}}
\hat{c}_{j_{1}\uparrow}^{\dag}\hat{c}_{j_{2}\uparrow}\hat{c}_{j_{3}\downarrow}^{\dag}\hat{c}_{j_{4}\downarrow}.\label{eq8}
\end{eqnarray}
Frankly speaking, the above Hamiltonian is hard to solve since adding external magnetic field breaks the translation invariance, which leads to the violation of solvability of original HK model (Eq.~\ref{eq1}).

Instead, following the original treatment of Luttinger,\cite{Luttinger1961,Wasserman1996} whose approximation is able to capture quantum oscillation in interacting electron systems such as Fermi liquid, we can still use results of last section but only need to replace single-particle energy $\varepsilon_{k}$ with eigen-energy of $\hat{H}_{0}$ (denoted as $\{E_{m}\}$). In the language of many-body physics, such replacement is a statement of that the external magnetic field only modifies the single-particle spectrum but not the expression of self-energy function.


Therefore, the ground-state energy density under external magnetic field ($\vec{B}=\nabla\times\vec{A}$) has the new formalism,
\begin{equation}
e_{g}=\frac{1}{N_{s}}\sum_{m}[(E_{m}-\mu)\theta(\mu-E_{m})+(E_{m}-\mu+U)\theta(\mu-E_{m}-U)],\label{eq9}
\end{equation}
and the related electron density is
\begin{equation}
n=\frac{1}{N_{s}}\sum_{m}[\theta(\mu-E_{m})+\theta(\mu-E_{m}-U)].\label{eq10}
\end{equation}
Thus, the magnetization and susceptibility in ground-state can be found by $M=-\frac{\partial e_{g}}{\partial B},\chi_{s}=\frac{\partial M}{\partial B}$. For thermodynamics at finite $T$, the only modification is $f_{k}\rightarrow f_{m}=1+2z_{m}+z_{m}^{2}e^{-\beta U}$ with $z_{m}=e^{-\beta(\varepsilon_{m}-\mu)}$ in Eq.~\ref{eq5}.


To be specific, let us consider the square lattice on $xy$-plane with nearest-neighbor-hopping ($t_{ij}\rightarrow t$) and the magnetic field is along the $z$-axis $\vec{B}=B\vec{e}_{z}$. In Landau gauge, the corresponding vector potential is $\vec{A}=Bx \vec{e}_{y}$, so $a_{i,i+x}=0$ and $a_{i,i+y}=\frac{e}{\hbar}Ba^{2}i_{x}$, ($a$ is the lattice constant and will be set to unit) which gives rise to
\begin{eqnarray}
\hat{H}_{0}=-t\sum_{i,\sigma}(\hat{c}_{i\sigma}^{\dag}\hat{c}_{i+x,\sigma}+e^{i\frac{e}{\hbar}Bi_{x}}\hat{c}_{i\sigma}^{\dag}\hat{c}_{i+y,\sigma}).\label{eq11}
\end{eqnarray}
and phase factor part can be parameterized as $e^{i\frac{e}{\hbar}Bi_{x}}=e^{i\frac{e}{\hbar}\Phi i_{x}}=e^{i2\pi w i_{x}}$ with flux $\Phi=2\pi\frac{\hbar}{e}w$ on each plaquette. In fact $\hat{H}_{0}$ is just the famous Hofstadter model with spin degree of freedom and has been widely studied in recent years due to its realization in cold atom experiments and relevance to moir\'{e} superlattice systems.\cite{Hofstadter1976,Aidelsburger2013,Miyake2013,Kennedy2015,Dean2013,Hunt2013,Spanton2018}


After diagonalization of $\hat{H}_{0}$, we obtain its eigen-energy $\{E_{m}\}$ for given magnetic field $B$, thus interested quantities can be calculated straightforwardly.
\section{Quantum oscillation in non-Fermi liquid: The case study on Hatsugai-Kohmoto model}\label{sec2}
\subsection{Lattice calculation}

In this section, we study the HK model under magnetic field on square lattice whose single-particle Hamiltonian is $\hat{H}_{0}$ (Eq.~\ref{eq11}). For parameters $t=e=\hbar=1$, the single-particle spectrum $\{E_{m}\}$ has been shown in Fig.~\ref{fig:1}, in which we consider a $60\times60$ system with periodic boundary condition. It is clear to see the well-known butterfly spectrum of $\hat{H}_{0}$ for $E_{m}$ versus flux $\Phi=2\pi w=B$.\cite{Hofstadter1976}
\begin{figure}
\includegraphics[width=0.75\linewidth]{fig_1.pdf}
\caption{\label{fig:1} The butterfly spectrum of $\hat{H}_{0}$ for $E_{m}$ versus flux $\Phi=2\pi w=B$.}
\end{figure}

Next, using Eq.~\ref{eq10}, we have plotted particle density $n$ versus $\Phi$ and $1/\Phi$ for different interaction strength $U$.
\begin{figure}
\includegraphics[width=0.85\linewidth]{fig_2.pdf}
\caption{\label{fig:2} The particle density $n$ versus $\Phi$ (a) and $1/\Phi$ (b) for different $U/t$.}
\end{figure}
To have regular data for our model on square lattice, we have used periodic boundary condition along $y$-direction while an open boundary condition is assumed for $x$-direction. Thus, a clear quantum oscillation of particle density appears in Fig.~\ref{fig:2}. (setting $\mu=0$) Similarly, in Fig.~\ref{fig:3} we have shown the data of magnetization, which also exhibits the quantum oscillation signal though with irregular error introduced by numerical differentiation. Therefore, we finds that the HK model indeed shows quantum oscillation phenomena in terms of the Luttinger's approximation.
\begin{figure}
\includegraphics[width=1.05\linewidth]{fig_3.pdf}
\caption{\label{fig:3} The magnetization $M$ versus $\Phi$ for $U/t=0.5,1.0,1.5,2.0$.}
\end{figure}

\subsection{Continuity limit}
After presenting the results of direct lattice calculation, in this subsection, we try to understand the origin of quantum oscillation from the continuity limit.


In continuity limit, one can approximate electron's dispersion $\varepsilon_{k}$ with free particle energy $\frac{k^{2}}{2m_{e}}$. At the same time,  $\varepsilon_{m}=\hbar\omega_{c}(m+1/2)$ if we consider usual Landau gauge. ($\omega_{c}=eB/m_{e}$) Then, the electron's density reads
\begin{equation}
n=\frac{1}{N_{\Lambda}}\sum_{k_{y}}\sum_{m=0}^{\Lambda}[\theta(\mu-\hbar\omega_{c}(m+1/2))+\theta(\mu-\hbar\omega_{c}(m+1/2)-U)].\label{eq12}
\end{equation}
Here, $\Lambda$ is the cutoff of Landau level index, which acts as a band-width. The factor $N_{\Lambda}$ has been introduced to give sensible definition of electron's density. The summation over $k_{y}$ gives the degeneracy of Landau level ($D$).

To proceed, we can rewrite $n$ as follows,
\begin{eqnarray}
n&=&\frac{D}{N_{\Lambda}}\lim_{\beta\rightarrow\infty}\int_{0}^{\infty}d\omega\frac{\sum_{m}\delta(\omega-\varepsilon_{m})+\sum_{m}\delta(\omega-\varepsilon_{m}-U)}{e^{\beta(\omega-\mu)}+1}\nonumber\\
&=&\frac{D}{N_{\Lambda}}\lim_{\beta\rightarrow\infty}\int_{0}^{\infty}d\omega\frac{1}{e^{\beta(\omega-\mu)}+1}\nonumber\\
&\times&\mathrm{Re}\int_{0}^{\infty}d\tau\sum_{m}(e^{i(\omega-\varepsilon_{m})\tau}
+e^{i(\omega-\varepsilon_{m}-U)\tau}).\label{eq13}
\end{eqnarray}
Then, using Poisson summation formula
$\sum_{m}e^{-im\hbar\omega_{c}\tau}=\frac{2\pi}{\hbar\omega_{c}}\sum_{p=0,\pm1,\pm2...}\delta(\tau-\frac{2\pi p}{\hbar\omega_{c}})$, and integrate over $\tau$, one finds
\begin{eqnarray}
n&=&\frac{2\pi D}{N_{\Lambda}\hbar\omega_{c}}\lim_{\beta\rightarrow\infty}\int_{0}^{\infty}d\omega\left(\frac{1}{e^{\beta(\omega-\mu)}+1}+\frac{1}{e^{\beta(\omega-\mu+U)}+1}\right)\nonumber\\
&\times&\mathrm{Re}\sum_{p}e^{i\frac{2\pi p}{\hbar\omega_{c}}(\omega-\frac{\hbar\omega_{c}}{2})}\nonumber\\
&\propto&\lim_{\beta\rightarrow\infty}\frac{2\pi}{\beta\hbar\omega_{c}}\mathrm{Re}\sum_{\omega_{n}>0}\sum_{p>0}(-1)^{p}\left(e^{i\frac{2\pi p}{\hbar\omega_{c}}(\mu+i\omega_{n})}+e^{i\frac{2\pi p}{\hbar\omega_{c}}(\mu-U+i\omega_{n})}\right)\nonumber\\
&=&\lim_{\beta\rightarrow\infty}\sum_{p>0}\frac{(-1)^{p}}{2\pi p}\frac{X}{\sinh X}\left(\cos\frac{2\pi p\mu}{\hbar\omega_{c}}+\cos\frac{2\pi p(\mu-U)}{\hbar\omega_{c}}\right).\label{eq14}
\end{eqnarray}
Here, $X=\frac{2\pi^{2}pT}{\hbar\omega_{c}}$ and usually one defines $R_{T}=\frac{X}{\sinh X}$ as damping factor for temperature. It is amusing to note that the above Eq.~\ref{eq14} is just the well-known LK formula with oscillation frequency $\mu/(\hbar\omega_{c})$ and $(\mu-U)/(\hbar\omega_{c})$. When $T\rightarrow0$ or $\beta\rightarrow\infty$, we obtain $n\propto\sum_{p>0}\frac{(-1)^{p}}{2\pi p}\left(\cos\frac{2\pi p\mu}{\hbar\omega_{c}}+\cos\frac{2\pi p(\mu-U)}{\hbar\omega_{c}}\right)$. Thus, it seems that under the approximation of Luttinger, despite the non-Fermi liquid nature, the quantum oscillation exists in the ground-state of HK model and it has two oscillation period determined by chemical potential $\mu$ and $\mu-U$.

\section{Exact calculation based on exact diagonolization}
Until now, our calculation is based on Luttinger's approximation, which is in principle not known whether it works or not.


Instead, one may use (numerical) exact diagonolization to extract reliable information from our model (Eq.\ref{eq8} and \ref{eq11}). Since $y$-direction of our model has been chosen to have periodic boundary condition under Landau gauge, we can use Fourier transformation to write
\begin{equation}
\hat{c}_{j\sigma}=\frac{1}{\sqrt{L_{y}}}\sum_{k_{y}}e^{ik_{y}j_{y}}\hat{c}_{j_{x}\sigma}(k_{y}),\label{eq15}
\end{equation}
where $x-$direction remains to be intact. Then, the Hamiltonian reads as $\hat{H}=\sum_{k_{y}}\hat{H}(k_{y})$,
\begin{eqnarray}
\hat{H}(k_{y})&=&-t\sum_{j_{x}\sigma}(\hat{c}_{j_{x}\sigma}^{\dag}(k_{y})\hat{c}_{j_{x}+1\sigma}(k_{y})+\hat{c}_{j_{x}+1\sigma}^{\dag}(k_{y})\hat{c}_{j_{x}\sigma}(k_{y}))\nonumber\\
&+&\sum_{j_{x}\sigma}\left(-2t\cos(\Phi j_{x}+k_{y})-\mu\right)\hat{c}_{j_{x}\sigma}^{\dag}(k_{y})\hat{c}_{j_{x}\sigma}(k_{y})\nonumber\\
&+&\frac{U}{L_{x}}\sum_{j_{1x},j_{2x},j_{3x},j_{4x}}\delta_{j_{1x}+j_{3x}=j_{2x}+j_{4x}}
\hat{c}_{j_{1x}\uparrow}^{\dag}(k_{y})\hat{c}_{j_{2x}\uparrow}(k_{y})\nonumber\\
&\times&\hat{c}_{j_{3x}\downarrow}^{\dag}(k_{y})\hat{c}_{j_{4x}\downarrow}(k_{y}).\label{eq16}
\end{eqnarray}
The above Hamiltonian means that for given momentum $k_{y}$, one can just study the reduced object $\hat{H}(k_{y})$, which is able to be solved in terms of exact diagonolization up to $L_{x}$ sites along $x-$direction (i.e. $j_{x}=1,2...Lx$).


In Fig.~\ref{fig:4}, we have shown the particle density $n$ versus $\Phi$ for $L_{x}=6$ and $L_{y}=16$. However, due to the limited size, particularly the smallness of $L_{x}$, although certain oscillation behavior appears, it is hard to detect clear quantum oscillation and extract its oscillation period and amplitude. Therefore, much larger size of $L_{x}$ should be used to
explore the true signature of quantum oscillation of HK model. (system with $L_x$ up to $16$ may be calculated in present workstation and it will be explored in our future work)
\begin{figure}
\includegraphics[width=0.85\linewidth]{fig_4.pdf}
\caption{\label{fig:4} The particle density $n$ versus $\Phi$ for $U/t=0,0.5,1.0,1.5,2.0$ from exact diagonolization.}
\end{figure}
\section{Conclusion and Future direction}\label{sec3}
In conclusion, we have taken Hatsugai-Kohmoto model as an example to study the quantum oscillation behavior in strongly correlated system. With Luttinger approximation, it is found that the quantum oscillation indeed appears in non-Fermi liquid state and one can use Lifshitz-Kosevich-like formula to extract its basic properties. We have also performed a small size exact diagonolization calculation, its results exhibits certain oscillation behavior but we cannot extract exact value of oscillation period and amplitude. Therefore, we expect that future work on this interesting issue is highly desirable, particularly the large scale numerical simulation.


\emph{Note added}: After completing this work, we have noticed the paper of  Leeb and Knolle,\cite{Leeb2023} whose results generally agree with ours though their calculation is based on the projection into Landau level.

