\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.


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

In longitudinal MRI studies, repeated MRI scans are performed in specific clinical scenarios, such as follow up of patients \cite{Angelini12},  therapy response assessment \cite{Weizman14}, and large-scale longitudinal studies  \cite{Schilling22}. They constitute one of the most efficient tools to track pathology changes and to evaluate treatment efficacy in diseases. In such studies, patients are scanned on a regular basis, which calls for efficient longitudinal scan protocols.

Conventional longitudinal MRI workflows consist of multiple steps. First, complex-valued image data is acquired (in k-space) at different time points, from which (also complex-valued) images are individually reconstructed. Next, the magnitudes of the resulting images are mutually geometrically aligned. Finally, local structural deformations are detected by computing deformation vector fields, allowing for voxel-wise deformation analysis \cite{Dufresne20}.

Unfortunately, conventional longitudinal MRI workflows are characterized by a low time efficiency. To accelerate imaging, techniques based on compressed sensing (CS) theory have been applied that exploit sparsity in some spatial transform domain and reconstruct images from substantially less data (through k-space sub-sampling) \cite{Lustig07}. Furthermore, since changes between two subsequent images in a longitudinal MRI study are generally small (apart from global transformations), temporal redundancy can be additionally exploited to shorten image acquisition time.  Indeed, CS-based methods have successfully been applied to also accelerate longitudinal MRI scanning by extending spatial regularization with temporal regularization \cite{Weizman15, Ye19, asif2013motion, jung2009k}. Such methods, however, require careful tuning of multiple regularization parameters to balance sparsity constraints (both in spatial and time domain) and data consistency, leading to impractical methods. Alternatively, deep learning methods have been investigated to automatically detect and exploit the distribution of expected images in both the spatial \cite{Tezcan19} and temporal \cite{yoo2021time, bustin2020compressed} domain. However, most data driven methods strongly rely on extensive training with a large number of datasets, which is especially non-trivial to obtain for longitudinal MRI applications.

In this paper, we propose DELTA-MRI, a novel framework to estimate longitudinal image deformations {\it directly} from a reference image and subsequently acquired, strongly sub-sampled MRI k-space data. In contrast to state-of-the-art longitudinal MR image reconstruction methods, which rely on the image reconstruction followed by deformation vector field estimation, our method avoids this multi-step process. Instead, the set of follow-up images, along with motion parameters and deformation vector fields that describe their relation to the reference image, are estimated in one go. We will show that DELTA-MRI performs significantly better compared to state-of-the-art CS based longitudinal MRI methods in terms of the normalized reconstructed error.

\section{methods}
Let $\bm{x}_1=\{x_{1,k}\}, \bm{x}_2=\{x_{2,k}\} \in \mathbb{C}^N$ be two complex valued 3D MRI images with $N$ voxels in the image domain, of which the elements can be expressed in polar form as
\begin{equation}
    x_{j,k} = r_{j,k} e^{i\phi_{j,k}} \quad,\ \textrm{with}\ j = 1, 2 \ ;\ k=1,...,N
\end{equation}
and with $\bm{r}_j=\{r_{j,k}\} \in \mathbb{R}^N$ and $\bm{\phi}_j=\{\phi_{j,k}\} \in \mathbb{R}^N$ the magnitude and phase of $\bm{x}_j$, respectively.
Furthermore, let $\bm{d}_1, \bm{d}_2 \in \mathbb{C}^N$ be the k-space (Fourier) representations of $\bm{x}_1, \bm{x}_2$, respectively: $\bm{d}_j = \bm F \bm{x}_j$, with

$\bm F \in \mathbb{C}^{N \times N}$
the discrete Fourier transform operator.

Finally, let  $\bm{v} \in \mathbb{R}^{N \times 3}$ be a deformation vector field (DVF), comprised of one 3D vector per voxel, describing the geometrical deformation from $\bm{r}_1$ to $\bm{r}_2$, and let $ W(\cdot, \bm{v})$ denote the image warping operator \cite{renders2021adjoint} along
$\bm{v}$ that (approximately) yields $\bm{r}_2$ when applied to $\bm{r}_1$:
\begin{equation}
    \label{eq:real_equation}
     W(\bm{r_1}, \bm{v}) = \bm{r}_2 \quad .
\end{equation}
We propose the estimation of $\bm{v}$ directly from k-space measurements of $\bm{x}_1$ and $\bm{x}_2$. These measurements, which will be denoted by $\{\tilde{\bm d_j}\}_{j=1}^2 \in \mathbb{C}^{n_j}$,
are noise disturbed and can be modelled as:
\begin{equation}
    \tilde{\bm d_j} =\bm{S}_j \bm{F} \bm{x}_j + \bm \varepsilon_j, \quad \textrm{with}\ j = 1, 2 \
\end{equation}
where $\bm \varepsilon_j \in \mathbb{C}^{n_j}$ is an additive noise contribution, modeled as a zero-mean, complex-valued Gaussian random variable, and $\bm{S}_{j}\in\{0,1\}^{n_j\times N }$, with $n_j \leq N$, is  a sub-sampling operator that selects the k-space points acquired for the $j$-th image.
For fully sampled k-space data, denoted by $\bm {\underline{\tilde d}}_j$, the operator $\bm{S}_j$ corresponds with the identity matrix $\bm I \in \mathbb{R}^{N \times N}$.

In what follows, it is assumed that sufficient k-space data $\tilde{\bm d_1}$ is acquired to allow a high quality reconstruction $\hat{\bm r}_1$ of the reference magnitude image $\bm r_1$. If fully sampled data $\bm {\underline{\tilde d}}_1$ is available, such an estimate can be calculated as
$\bm {\hat{r}}_1=|\bm{F}^{-1} \underline{\tilde {\bm d}}
_1|$, with $|\cdot|$ the point-wise modulus operator.
Furthermore, it is assumed that the phase $\bm{\phi}_2$ is slowly varying and that the low frequencies in k-space are sufficiently sampled to obtain a good estimate of $\bm{\phi}_2$ by calculating $\hat{\bm \phi}_2 =\angle (\bm F^{-1} \tilde{\bm{d}}_2)$, where $\tilde{\bm{d}}_2$ is zero-filled at the non-sampled indices and $\angle(\cdot)$ denotes the phase operator.
Under these assumptions, $\bm{v}$ can be directly estimated
by minimizing a regularized least-squares functional
without requiring the reconstruction of $\bm{x}_2$:
\begin{equation}
\label{eqn:estimator_v}
   \hat{\bm{v}}=  \arg \min_{\bm{v}} \left( \lVert \tilde{\bm{d}_2}-\bm{S}_2 \bm{F} W(\bm{\hat{r}_1}, \bm{v})e^{i\hat{\bm{\phi}}_2} \rVert_2^2 + \lambda R(\bm{v})\right),

\end{equation}
with
$\lambda R(\bm{v})$ a regularization term with weight $\lambda$ that
imposes prior knowledge about
$\bm{v}$.
In this work,
$R(\bm{v})$ is chosen equal to $\|\nabla \bm v\|_2$, with $\nabla$ the gradient operator. This choice
reflects the assumption that $\bm{r}_2$ only differs from $\bm{r}_1$ by a local, smooth deformation, for which this regularization term is small.

In practice, $\bm{x}_2$ will likely also differ from $\bm{x}_1$ by a rigid transformation, on top of the local deformation. This transformation can be included in the DVF $\bm{v}$, but then $\bm{v}$ will no longer be mostly zero. Alternatively, the rigid transformation can be isolated from the DVF, and represented by rotation and translation parameters $\bm{\theta}, \bm{t} \in \mathbb{R}^3$.
\refeq{eqn:estimator_v}
can then be extended as follows:
\begin{eqnarray} \label{eqn:joint_estimator}
   \{\bm{\hat  v, \hat \theta, \hat t}\}&=&\arg \min_{\bm{\theta}, \bm{t}, \bm{v}} ( \lVert  \tilde{\bm{d}}_2 - \bm{SF} W(\hat{\bm{r}}_1, \bm{\theta}, \bm{t, \bm{v}})e^{i\hat{\bm{\phi}}_2} \rVert_2^2 \nonumber \\
    &&+ \lambda R(\bm{v})) \quad
    ,
\end{eqnarray}
where $W(\cdot, \bm{\theta}, \bm{t, \bm{v}})$ is an image warping operator that rotates the image according to $\bm{\theta}$, translates the image according to $\bm{t}$ and then applies the local deformation described by $\bm{v}$. The optimization problem of \refeq{eqn:joint_estimator} can be solved using the cyclic Block Coordinate Descent (cBCD) method \cite{Beck13}, where the first block consists of the parameters $\bm{\theta}$ and $\bm t$ and the second block consist of $\bm{v}$.  In what follows, the initial estimates of  $\bm{v}$, $\bm t$ and $\bm \theta$ will be denoted as $\bm{v}_0$, $\bm{t}_0$ and $\bm{\theta}_0$, respectively.

\section{Experiments}
\label{sec:experiments}
High-resolution, fully-sampled k-space data was acquired using a Cartesian sampling scheme of a male C57BL/6 wild-type mouse using a 9.4T Biospec 94/20 USR horizontal MR system (Bruker Biospin MRI, Ettlingen, Germany), with a mouse head 2$\times$2 array cryo-coil (Bruker Biospin MRI, Ettlingen, Germany), using a 3D T2-weighted Turbo-RARE sequence with FOV: $20\times 15\times 10$ mm$^3$,
image acquisition matrix: $256\times 192\times 128$, 0.078 mm isotropic voxels, TE/TR: 67.2/1800 ms, RARE factor: 12. The fully-sampled multi-coil k-space data was then averaged across the coils to obtain
$\underline{\tilde{\bm{d}}}_2$.

Based on the experimentally acquired dataset, simulation experiments were setup to validate DELTA-MRI. First, a complex valued  image ${\bm x}_2=\bm{F}^{-1} \underline{\tilde{\bm{d}}}_2$ was reconstructed, with corresponding magnitude ${\bm r}_2$ and phase $\bm{\phi}_2$.
From $\bm{r}_2$, a predefined rigid transform ($\bm{\theta}$, $\bm t$), and a predefined DVF $\bm{v}$, a magnitude image $\bm{r}_1$ was constructed
using an inverse image warp, such that
$W(\bm r_1, \bm \theta, \bm t, \bm v) = \bm r_2$.
 \reffig{fig:GT_DVF_color} shows a cross section of the ground truth DVF $\bm v$.   The ground truth parameter
  values of $\bm{\theta}$ and $\bm t$ were given by $(2.9^\circ, 4.0^\circ, 5.7^\circ)$ and $(-6, -5, -4.5)$ pixels, respectively.
\reffig{fig:r1}, \ref{fig:r1aligned}, and \ref{fig:r2} show cross sections of $\bm{r}_1$, its fully aligned version $W(\bm {r}_1, \bm \theta, \bm t, \bm 0)$, and $\bm{r}_2$, respectively. Since DELTA-MRI does not depend on the phase of the reference image, its evaluation does not require the simulation of $\bm{\phi}_1$. However, as will be described below, DELTA-MRI was benchmarked against a state-of-the-art method that does depend on $\bm{\phi}_1$ and assumes it to be similar to $\bm{\phi}_2$.
Therefore, a complex valued image ${\bm x}_1$ with magnitude $\bm{r}_1$ was constructed, by setting $\bm {\phi}_1 = \bm{\phi}_2$.
This complex image was polluted with complex, Gaussian distributed white noise with standard deviation equal to 4\% of the average foreground value of $\bm{r}_1$. We denote this noise disturbed image by $\hat{\bm {x}}_1$ and its magnitude and phase by $\hat{\bm {r}}_1$ and $\hat{\bm{\phi}}_1$, respectively.
Next, a sub-sampling operator $\bm S_2$ was applied that under-samples $\underline{\tilde{\bm {d}}}_2$ by selecting a given percentage of the k-space data points. In particular, a small cubic central area was complemented by k-space points that were randomly drawn from a central Gaussian distribution with a standard deviation that corresponds with a quarter of the maximum frequency.
Illustrative examples of the thus designed sampling schemes are shown in \reffig{fig:normal1} and \reffig{fig:normal2}.
Finally,  $\bm v$, $\bm{\theta}$ and $\bm t$, were estimated from the under-sampled k-space data $\tilde{\bm {d}}_2$ using the estimator described by \refeq{eqn:joint_estimator}. The
experiments were repeated for different sub-sampling percentages.

\begin{figure}[hbt]
    \centering
    \begin{subfigure}[t]{2.8cm}
    \includegraphics[width=2.8cm]{images/gt_color.png}
    \caption{Ground truth DVF}
    \label{fig:GT_DVF_color}
    \end{subfigure}
    \begin{subfigure}[t]{2.8cm}
    \includegraphics[width=2.8cm]{images/s_normal1.png}
    \caption{Gaussian 1\%}
    \label{fig:normal1}
    \end{subfigure}
    \begin{subfigure}[t]{2.8cm}
    \includegraphics[width=2.8cm]{images/s_normal5.png}
    \caption{Gaussian 5\%}
    \label{fig:normal2}
    \end{subfigure}
    \caption{Ground truth DVF (a) Sampling schemes (b) and (c).}
    \label{fig:sampling_schemes}
\end{figure}

\begin{figure}[htb]
    \centering
    \begin{subfigure}[t]{2.8cm}
    \includegraphics[width=2.8cm]{images/r1.png}
    \caption{$\bm{r}_1$}
    \label{fig:r1}
    \end{subfigure}
    \begin{subfigure}[t]{2.8cm}
    \includegraphics[width=2.8cm]{images/aligned_r1_arrow.png}
    \caption{$\bm{r}_1$ aligned}
    \label{fig:r1aligned}
    \end{subfigure}
    \begin{subfigure}[t]{2.8cm}
    \includegraphics[width=2.8cm]{images/r2_arrow.png}
    \caption{$\bm{r}_2$}
    \label{fig:r2}
    \end{subfigure}
    \caption{Reconstructions from fully sampled k-space data.}
    \label{fig:r1_and_r2}
\end{figure}

\textbf{Implementation details:}
The optimization task of \refeq{eqn:joint_estimator} was carried out using a cBCD approach with
$\bm{v}_0$, $\bm{t}_0$ and $\bm{\theta}_0$ all
zero. The value of $\lambda$ was set to $\lambda_0 \cdot 5 \cdot 10^{-5}$, where $\lambda_0 = \lVert \tilde{\bm{d}_2} \rVert_2^2$ maintains the regularization strength when varying the sampling scheme. This value was obtained by empirically minimizing the reconstruction error at 1\% Gaussian sub-sampling. In the first block, $\bm{\theta}$ and $\bm t$ were optimized with the \mbox{L-BFGS-B} method of SciPy 1.9.3 \cite{2020SciPy-NMeth}, with default stopping criteria and bounds $\bm t \in [-20, 20]^3, \bm{\theta} \in [-0.3 \text{ rad}, 0.3 \text{ rad}]^3$. In the second block,  $\bm{v}$ was optimized with the Barzilai-Borwein gradient method \cite{barzilai1988two} with 2000 iterations to ensure convergence. The derivatives of the objective function towards the deformation parameters were computed using ImWIP \cite{renders2022imwip, renders2021adjoint}. Because the deformation is local, the initial estimate
$\bm{v}_0$ was close enough to give an accurate estimation of $\bm{\theta}$ and $\bm t$. Therefore, each block was optimized only once. It was verified that adding more cycles did not change the results significantly.

\textbf{Benchmarking:} DELTA-MRI was benchmarked against
the Temporal Compressed Sensing MRI  (TCS-MRI) method \cite{Weizman15},
which regularizes the reconstruction of the complex image $\bm{x}_2$ in the temporal and wavelet domain, by solving
\begin{eqnarray}
\label{eqn:TCS-MRI}
\bm{\hat x}_2 &=&  \arg \min_{\bm{x}_2} ( \lVert \tilde{\bm{d}_2}-\bm{S}_2 \bm{F} \bm{x}_2 \rVert_2^2   \nonumber \\
&& + \lambda_1 \lVert \text{diag}(\bm{w}) (\hat{\bm{x}}_1 - \bm{x}_2) \rVert_1 + \lambda_2 \lVert \bm{\Psi} \bm{x}_2 \rVert_1 ),
\end{eqnarray}
with $\bm{\Psi}$ the Daubechies-4 wavelet transform and $\bm{w}$ a weight vector that controls the demand for similarity between $\bm{x}_1$ and $\bm{x}_2$, enforcing sparsity only in regions where $\bm{x}_1$ and $\bm{x}_2$ are similar. In our experiments, $\bm{w}$ was fixed to the support of the ground truth DVF.

It should be noted that TCS-MRI assumes that the rigid transformation aligning $\bm{x}_1$ and $\bm{x}_2$ is known.
DELTA-MRI
does not make this assumption, and instead estimates the rigid transformation, along with $\bm v$,
from the available, sub-sampled data. In each experiment, the rigid transformation estimated by DELTA-MRI was used to align the data before using TCS-MRI. Furthermore, the regularization weights $\lambda_1$ and $\lambda_2$ in  \refeq{eqn:TCS-MRI} were set to $\lambda_0 \cdot 5.85 \cdot 10^{-9}$ and $\lambda_0 \cdot 3.63 \cdot 10^{-10}$, respectively, where $\lambda_0 = \lVert \tilde{\bm{d}_2} \rVert_2^2$ maintains the regularization strength when varying the sampling scheme. These values were obtained by empirical tuning to optimize reconstruction quality at 1\% Gaussian sampling.

To evaluate how much information on $\bm{r}_2$ was present in $\tilde{\bm{d}_2}$, and how much information was obtained from $\hat{\bm{x}}_1$, the results of TCS-MRI and DELTA-MRI were compared against a reconstruction
obtained from $\tilde{\bm{d}_2}$ alone by applying an inverse DFT
to the zero
filled data $\tilde{\bm{d}_2}$. This method is denoted Z-IDFT.

\textbf{Performance criterion:}
As a quantitative performance criterion, the normalized reconstruction error of the reconstructed image
$\hat{\bm{r}}_2$ was calculated as:
\begin{equation} \label{eqn:norm reconstruction error}
\epsilon(\hat{\bm r}_2) = \|\hat{\bm r}_2-\bm {r}_2 \|_2 / \|\hat{\bm {r}}_1-\bm{r}_2\|_2
    \quad ,
\end{equation}
where $\hat{\bm{r}}_2 := W(\hat{\bm r}_1, \hat{\bm \theta}, \hat{\bm t}, \hat{\bm v})$ for DELTA-MRI.

\section{Results \& Discussion}
\label{sec:results_discussion}
Figure \ref{fig:rec_normal}
shows reconstruction results obtained from Gaussian sub-sampled k-space data
with sub-sampling percentages equal to 5\% and 1\%.
The left column shows
(a cross-section of) the Z-IDFT reconstruction of $\bm{r}_2$.
The middle and right column show the reconstructions of $\bm{r}_2$ obtained by TCS-MRI and
DELTA-MRI, respectively.
Figure \ref{fig:quantitative comparison} shows the normalized reconstruction error (\ref{eqn:norm reconstruction error}) as a function of the sub-sampling percentage for Z-IDFT, TCS-MRI and DELTA-MRI.

\begin{figure}[ht]
    \centering
    \begin{subfigure}[t]{2.8cm}
    \includegraphics[width=2.8cm]{images/normal5/direct_arrow.png}
    \caption{Z-IDFT}
    \label{fig:Direct_5proc_norm}
    \end{subfigure}
    \begin{subfigure}[t]{2.8cm}
    \includegraphics[width=2.8cm]{images/normal5/TCS_MRI_r2_arrow.png}
    \caption{TCS-MRI}
    \label{fig:TCS-MRI_5proc_norm}
    \end{subfigure}
    \begin{subfigure}[t]{2.8cm}
    \includegraphics[width=2.8cm]{images/normal5/warped_r1_arrow.png}
    \caption{DELTA-MRI}
    \label{fig:Delta_MRI_5proc_norm}
    \end{subfigure}
    \newline
    \begin{subfigure}[t]{2.8cm}
    \includegraphics[width=2.8cm]{images/normal1/direct_arrow.png}
    \caption{Z-IDFT}
    \label{fig:Direct_1proc_norm}
    \end{subfigure}
    \begin{subfigure}[t]{2.8cm}
    \includegraphics[width=2.8cm]{images/normal1/TCS_MRI_r2_arrow.png}
    \caption{TCS-MRI}
    \label{fig:TCS_MRI_1porc_norm}
    \end{subfigure}
    \begin{subfigure}[t]{2.8cm}
    \includegraphics[width=2.8cm]{images/normal1/warped_r1_arrow.png}
    \caption{DELTA-MRI}
    \label{fig:Delta_MRI_1proc_norm}
    \end{subfigure}
    \caption{Reconstructions $\hat{\bm{r}}_2$ from 5\% (top) and 1\% (bottom) sub-sampled data.}
    \label{fig:rec_normal}
\end{figure}
\begin{figure}[!htb]
    \centering
    \includegraphics[width=8cm]{images/error_plot.pdf}
    \caption{
    Normalized reconstruction error of Z-IDFT, TCS-MRI and DELTA-MRI as a function of the sub-sampling percentage of a Gaussian sub-sampling scheme.
    }
    \label{fig:quantitative comparison}
\end{figure}
The results show that for sub-sampling percentages up to 20\%, DELTA-MRI outperforms Z-IDFT and TCS-MRI
. Furthermore, \reffig{fig:quantitative comparison} shows that for a sub-sampling percentage below around 13\%, TCS-MRI outperforms Z-IDF, whereas  for higher sub-sampling percentages Z-IDF outperforms TCS-MRI. This
can be explained as follows. In our simulation study, the
weights of the regularization terms and data-misfit term that together constitute the cost function of \refeq{eqn:TCS-MRI} were kept
proportional to each other
and were tuned to be optimal for a sub-sampling percentage of 1\%. For increasing sub-sampling percentages, these regularization settings
lead to over-regularization, causing TCS-MRI to be outperformed by Z-IDF beyond a certain turning point.
Alternatively, the regularization settings of  TCS-MRI and Delta-MRI could be optimized for each sub-sampling percentage individually, resulting in reduced reconstruction errors for all sub-sampling percentages excluding 1\%. However, this approach would not be realistic, as in practice ground truth images are not available and regularization settings are typically sub-optimal.

Additionally, some remarks can be made about the comparison between DELTA-MRI and TCS-MRI. First, unlike DELTA-MRI, TCS-MRI assumes the phases of $\bm{x}_1$ and $\bm{x}_2$ to be similar. This assumption is reflected by the first regularization term of \refeq{eqn:TCS-MRI}, which enforces similarity of $\bm{x}_1$ and $\bm{x}_2$ in the complex domain. In the simulation experiments, the ground truth phases $\bm{\phi}_1$ and $\bm{\phi}_2$ were chosen equal to meet this assumption, providing optimal conditions for  TCS-MRI.
In practice, however, the assumption of similar phase images may be
violated, which may negatively affect the performance of TCS-MRI. Second, as mentioned above, unlike DELTA-MRI, TCS-MRI does not include an alignment procedure, assuming that the rigid transformation aligning the images $\bm{x}_1$ and $\bm{x}_2$ is known beforehand. To allow a fair comparison that accounts for this difference, the rigid transformation estimated by DELTA-MRI was used to align the data before applying TCS-MRI. TCS-MRI could have been combined with alternative alignment procedures
, but this is not expected to change the main conclusions of the comparison.

The current work has some limitations that will be addressed in future work. First,
while the current version of DELTA-MRI was designed for single-coil data, the method can be extended to account for multi-coil data, taking coil sensitivities into account.
Furthermore, to facilitate experiments with prospectively under-sampled longitudinal data, the Gaussian distributed sampling schemes used to evaluate DELTA-MRI in this work can be replaced by sampling schemes that are practically more feasible.

Finally, it is worthwhile mentioning that the problem of determining a DVF that deforms a source image to a target image is underdetermined. In other words, many DVFs may transform the source image to the same target image. This is known as the aperture problem \cite{Marr82,Horn93}. For this reason, in simulation experiments, the quality of a DVF estimate is not appropriately measured in terms of its deviation from the ground truth DVF, but in terms of the difference between the ground truth target image and the image that is obtained by geometrically deforming the source image according to the estimated DVF. This motivates the use of the reconstruction error \refeq{eqn:norm reconstruction error} as a performance criterion in this work.

\section{Conclusion}
\label{sec:conclusion}
We proposed DELTA-MRI, a novel framework that estimates longitudinal image changes directly from a reference image and strongly sub-sampled k-space data of a follow-up image.
The DVF and motion parameters that describe the relation between the reference image and the follow-up image are estimated in one-go.
For sub-sampling percentages up to 20\%, DELTA-MRI was shown to outperform a state-of-the-art longitudinal CS based method in terms of reconstruction error.
\section{Compliance with ethical standards}
\label{sec:ethics}
This study was performed in line with the principles of the Declaration of Helsinki. All applicable institutional and/or national guidelines for the care and use of animals were followed. Procedures were performed in accordance with the European Directive 2010/63/EU on the protection of animals used for scientific purposes. The protocols were approved by the Committee on Animal Care and Use at the University of Antwerp, Belgium (permit number: 2016-10) and all efforts were made to minimize animal suffering. Experiments have been reported in compliance with the ARRIVE guidelines.
\section{Acknowledgments}
\label{sec:acknowledgments}
The authors thank J. van Rijswijk and J. Van Audekerke (Bio-Imaging Lab, Univ. Antwerp, Belgium) for the acquisition of the experimental MR data and valuable discussions. Jens Renders is an SB PhD fellow at the Research Foundation - Flanders (FWO), grant no. 1SA2920N. This work is partially supported by FWO project S007219N.

\begingroup
\raggedright

\bibliographystyle{IEEEbib-abbrev}
{