\section{Conclusion}
This paper proposes a general network slimming pipeline named ClassPruning, which could speed up existing learning-based image restoration models by approximately 1.67 times while maintaining performance.
The key idea is first to use the Classifier Module to determine the restoration difficulties of the input images.
And then, we perform dynamic $N$:$M$ fine-grained structured pruning on the restoration model based on the predicted restoration difficulty, thus obtaining the corresponding sparse subnetworks to process the input images.
Further, we propose a novel training strategy and two additional loss terms to stabilize training and improve performance.
Extensive experiments on various tasks and datasets demonstrate that our ClassPruning could save approximately 40\% FLOPs while introducing only a few parameters, indicating the effectiveness and generalization of ClassPruning.

\section{Experiment}
We evaluate our proposed ClassPruning on four different image restoration tasks: (\textbf{a}) image denoising, (\textbf{b}) image deraining, (\textbf{c}) image deblurring, and (\textbf{d}) image dehazing.
\subsection{Datasets}

\noindent
\textbf{Image Denoising.} For gaussian image denoising, we follow ~\cite{chen2016trainable,zhang2017beyond} to use 400 images of size 180 $\times$ 180 for training and BSD68~\cite{roth2009fields} for testing.
Noisy images are generated by adding additive white Gaussian noise with noise level $\sigma$ to clean images, and we set $\sigma$ to be 15, 25, and 50 for both training and testing.
And for real image denoising, we use 320 high-resolution images of the SIDD dataset~\cite{abdelhamed2018high} to train our model and perform the evaluation on 1280 patches from the SIDD validation set.

\noindent
\textbf{Image Deraining.} Following ~\cite{li2018recurrent}, we select 700 images from Rain800~\cite{zhang2019image} to train our model and perform the evaluation on the remaining 100 images, named Test100, for testing. Also, we use Rain100H~\cite{yang2017deep} as the testing dataset to evaluate our model.
Note that the PSNR/SSIM scores are computed on the RGB channels instead of the Y channel.

\noindent
\textbf{Image Deblurring.} Consistent with prior works~\cite{nah2017deep,tao2018scale}, we use the GoPro~\cite{nah2017deep} dataset for training our model. It contains 2,103 blurry sharp image pairs for training and 1,111 for validation. We also use Kohler dataset~\cite{kohler2012recording} for evaluation, which consists of 4 images, and each image is blurred with 12 different kernels.

\noindent
\textbf{Image Dehazing.} Prior work~\cite{li2018benchmarking} proposes an image dehazing benchmark RESIDE, which consists of large-scale training and testing hazy image pairs synthesized from depth and stereo datasets.
We leverage the RESIDE's Indoor Training Set for training, which contains 1399 clean images and 13990 hazy images generated by corresponding clean images.
And we evaluate the models' performance on Synthetic Objective Testing Set (SOTS), which contains 500 indoor images and 500 outdoor ones.

\subsection{Comparison with Base Models}

\begin{table*}[t]
  \caption{Comparison between the base models and ClassPruning networks on denoising, deraining, deblurring, and dehazing tasks. The CP-$Model$ means the model is obtained using our proposed ClassPruning, where $Model$ corresponds to the specific model name.}
    \centering
    \resizebox{\linewidth}{!}{
    \small
    \begin{tabular}{lccccccccc}
        \hline
        \multirow{2}*{{Methods}}&\multicolumn{4}{c}{{\hspace{-0.75cm} BSD68}} & &\multicolumn{4}{c}{\hspace{-0.75cm} {SIDD}}\\
        \cline{2-5}
        \cline{7-10}
        & {$\sigma=15$} ~~~~ {$\sigma=25$}&{$\sigma=50$} & {FLOPs} & \#Params & &{PSNR}&~~~~~~~{SSIM}& {FLOPs} & \#Params\\
        \hline
        DnCNN & 31.63 ~~~~~~~ 29.18 & 26.24 &87.28G (100\%) &0.668M & &38.43 &~~~~~~~0.910 &87.58G (100\%)& 0.671M\\
        CP-DnCNN& 31.59 ~~~~~~~ 29.17 &26.26 &50.88G (58.3\%) &0.754M & &38.38 &~~~~~~~0.908 &53.33G (60.9\%) &0.758M \\
        \bottomrule
        \toprule
        \multirow{2}*{{Methods}}&\multicolumn{4}{c}{\hspace{-0.75cm}{Test100}} & &\multicolumn{4}{c}{\hspace{-0.75cm}{Rain100H}}\\
        \cline{2-5}
        \cline{7-10}
        & \hspace{-1.3cm} {PSNR}&\hspace{-1cm} {SSIM} & {FLOPs} & \#Params & &{PSNR}&~~~~~~~{SSIM} & {FLOPs} & \#Params\\
        \hline
        RESCAN & \hspace{-1.3cm}24.46 & \hspace{-1cm}0.838 & 64.50G (100\%)&0.589M & &25.71 &~~~~~~~0.832 &64.50G (100\%)&0.589M\\
        CP-RESCAN& \hspace{-1.3cm}24.38 &\hspace{-1cm}0.838 &36.18G (56.1\%)&0.672M & &25.60 &~~~~~~~0.830 &34.31G (53.2\%)&0.672M\\
        \bottomrule
        \toprule
        \multirow{2}*{{Methods}}&\multicolumn{4}{c}{\hspace{-0.75cm}{GoPro}} & &\multicolumn{4}{c}{\hspace{-0.75cm} {Kohler}}\\
        \cline{2-5}
        \cline{7-10}
        & \hspace{-1.3cm} {PSNR}& \hspace{-1cm}{SSIM} & {FLOPs}& \#Params & &{PSNR}&~~~~~~~{SSIM} & {FLOPs}& \#Params\\
        \hline
        SRN & \hspace{-1.3cm} 30.21 & \hspace{-1cm}0.933& 217.2 (100\%)&10.25M & &26.70 &~~~~~~~0.835 &217.2 (100\%)&10.25M\\
        CP-SRN& \hspace{-1.3cm} 30.18 & \hspace{-1cm}0.932& 131.8 (60.7\%)&10.33M & & 26.71 &~~~~~~~0.834 &132.7 (61.1\%)&10.33M\\
        \bottomrule
        \toprule
        \multirow{2}*{{Methods}}&\multicolumn{4}{c}{\hspace{-0.75cm}{RESIDE-IN}} & &\multicolumn{4}{c}{\hspace{-0.75cm}{RESIDE-OUT}}\\
        \cline{2-5}
        \cline{7-10}
        & \hspace{-1.3cm} {PSNR}& \hspace{-1cm}{SSIM} & {FLOPs}& \#Params & &{PSNR}&~~~~~~~{SSIM} & {FLOPs}& \#Params\\
        \hline
        GCANet & \hspace{-1.3cm} 30.21 & \hspace{-1cm} 0.979& 36.82 (100\%)&0.702M & & 30.36  &~~~~~~~0.963 &36.82 (100\%)&0.702M\\
        CP-GCANet& \hspace{-1.3cm} 30.17 & \hspace{-1cm}0.977& 22.24 (60.4\%)&0.785M & & 30.35 &~~~~~~~0.959 &21.68 (58.9\%)&0.785M\\
        \bottomrule
    \end{tabular}
    }
    \label{tab:Compare}
\end{table*}

To demonstrate the effectiveness and generalization of the method, we use ClassPruning to speed up some representative networks for various image restoration tasks.
Specifically, we use DnCNN~\cite{zhang2017beyond}, RESCAN~\cite{li2018recurrent}, SRN~\cite{tao2018scale} and GCANet~\cite{chen2019gated} for denoising, deraining, deblurring and dehazing, respectively.
Training and testing follow the strategy described above, and results are summarized in Table~\ref{tab:Compare}.
In summary, ClassPruning helps models obtain comparable performance to the base networks and accelerate all base models approximately 1.67 times (40\% FLOPs reduction) by introducing only a few parameters, mainly belonging to the Classifier Module.
Also, we find that the reduction of FLOPs is correlated with the restoration task and the test dataset.
With comparable performance, the Kohler dataset in the deblurring task is the most difficult, reducing only 39\% FLOPs, while the Rain100H dataset in the deraining task is the easiest, reducing 45\% FLOPs.

\subsection{Comparison with Other Dynamic Networks}

\begin{table}[t]
  \caption{Comparison with different dynamic approaches. Multi-Width and Multi-Depth are formed by multiple independent networks with different widths and depths. Channel-Skip and Layer-Skip are obtained by skipping layers and channels of the network.}
  \centering
  \resizebox{\linewidth}{!}{
  \begin{tabular}{lccccc}
    \hline
      \multirow{2}*{{Method}}&\multicolumn{5}{c}{BSD68}\\
      \cline{2-6}
      &{$\sigma=15$} & {$\sigma=25$} & {$\sigma=50$}& {FLOPs} & {\#Params}\\
    \hline
    Multi-Width~\cite{eigen2013learning} & 31.47& 29.15& 26.23 & 50.77G& 1.251M\\
    Multi-Depth~\cite{park2015big} & 31.60 & 29.17 & 26.27 & 52.11G&1.272M\\
    Channel-Skip~\cite{li2021dynamic}  &31.38 & 29.03 & 26.21 &50.71G&0.754M\\
    Layer-Skip~\cite{wu2018blockdrop} &31.42 & 29.08 & 26.24 &51.58G&0.754M\\
    Ours &31.59 &29.17 &26.26 &50.88G& 0.754M\\
    \bottomrule
  \end{tabular}
  }
  \label{tab:Dynamic Ablation}
\end{table}

We compare ClassPruning with other dynamic networks on the denoising task to demonstrate the effectiveness of our method.
And here, we adopted four different dynamic networks: Multi-Width~\cite{eigen2013learning}, Multi-Depth~\cite{park2015big}, Channel-Skip~\cite{li2021dynamic} and Layer-Skip~\cite{wu2018blockdrop}.
The first two selectively execute one of the multiple independent candidate networks with different depths and widths, respectively.
And the last two skip the execution of intermediate layers and channels, respectively.
The results are shown in Tabel~\ref{tab:Dynamic Ablation}, and we can see that our method outperforms most dynamic approaches, especially the channel-skip, which performs channel-wise structured pruning.
We also found that dynamic networks ensembling all the subnetworks in one model are more challenging to optimize than independent subnetworks, and our training strategy could be a good solution.
Although our method is slightly inferior to the multi-depth dynamic approach, it is worth noting that ClassPruning introduces fewer parameters, \emph{i.e.} more lightweight, and the $N$:$M$ sparsity pattern makes our method more friendly to hardware.

\subsection{Analysis}

\begin{figure*}[t]
  \centering
  \includegraphics[width=1.0\textwidth]{Figs/Results.pdf}
  \centering
  \caption{Visual samples obtained by sparse networks with different capabilities. The purple outer box corresponds to the difficulty class predicted by our Classifier Module.}
  \label{Fig.Result}
\end{figure*}

\noindent
\textbf{Visualization.}
Figure~\ref{Fig.Result} shows some samples of results obtained by sparse models with different capabilities, and the purple box indicates the branch selected by our Classifier Module.
We can see that lightly degraded images can be restored well by sparse networks with high pruning ratios, and only minor improvements are achieved with more complex models.
However, when we employ sparse networks with high pruning ratios to handle highly degraded images, the performance drops dramatically and unsatisfying results are produced.
Also, it can be seen that our method appropriately assigns networks with different capacities to images with varying restoration difficulties, thus reducing the average computational cost while maintaining performance.

\begin{table}[t]
  \caption{Comparison between training from scratch ($^\dagger $) and training with our proposed strategy. Note that we train from scratch multiple times and select the best results due to its instability.}
  \centering
  \resizebox{\linewidth}{!}{
  \large
  \begin{tabular}{lccc}
      \hline
      \multirow{2}*{{Methods}}&\multicolumn{3}{c}{\hspace{-1.3cm}{BSD68}}\\
      \cline{2-4}
      & {$\sigma=15$} ~~ {$\sigma=25$}&{$\sigma=50$} & {FLOPs}\\
      \hline
      CP-DnCNN$^\dagger $&31.43 ~~~~ 29.01&26.08 &51.06G (58.5\%)  \\
      CP-DnCNN& 31.59 ~~~~ 29.17 &26.26 &50.88G (58.3\%) \\
      \bottomrule
      \toprule
      \multirow{2}*{{Methods}}&\multicolumn{3}{c}{\hspace{-1.3cm}{Test100}}\\
      \cline{2-4}
      & \hspace{-1.3cm}{PSNR}& \hspace{-1.3cm} {SSIM} & {FLOPs}\\
      \hline
      CP-RESCAN$^\dagger $& \hspace{-1.3cm} 24.30 &\hspace{-1.3cm} 0.835&37.47G (58.1\%)\\
      CP-RESCAN& \hspace{-1.3cm} 24.38 &\hspace{-1.3cm} 0.838 &36.18G (56.1\%)\\
      \bottomrule
      \toprule
      \multirow{2}*{{Methods}}&\multicolumn{3}{c}{\hspace{-1.3cm} {GoPro}}\\
      \cline{2-4}
      & \hspace{-1.3cm} {PSNR}&\hspace{-1.3cm}  {SSIM} & {FLOPs}\\
      \hline
      CP-SRN$^\dagger $& \hspace{-1.3cm} 30.05 & \hspace{-1.3cm} 0.928& 130.01 (59.9\%) \\
      CP-SRN& \hspace{-1.3cm} 30.18 &\hspace{-1.3cm}  0.932& 131.83 (60.7\%)\\
      \bottomrule
      \toprule
      \multirow{2}*{{Methods}}&\multicolumn{3}{c}{\hspace{-1.3cm} {Indoor}}\\
      \cline{2-4}
      & \hspace{-1.3cm} {PSNR}&\hspace{-1.3cm}  {SSIM} & {FLOPs}\\
      \hline
      CP-GCANet$^\dagger $&\hspace{-1.3cm} 30.02 & \hspace{-1.3cm} 0.970 &21.94 (59.6\%)\\
      CP-GCANet&\hspace{-1.3cm} 30.17 &\hspace{-1.3cm}  0.977& 22.24 (60.4\%)\\
      \bottomrule
  \end{tabular}
  }
  \label{tab:Strategy}
\end{table}

\noindent
\textbf{Training Strategy.}
We obtained the restoration models by training with our proposed strategy and from scratch, respectively, to explore the impact of the training strategy on various tasks.
It is worth noting that training from scratch is quite unstable and sensitive to the weight ratio of $L_C$ and $L_W$, leading the Classifier Module to classify all inputs into one restoration difficulty easily; here, we train from scratch multiple times and select the best results.
Table~\ref{tab:Strategy} shows that our strategy can trade off performance and efficiency better, especially in denoising and deraining tasks.
Also, the novel strategy enables our training to be more stable and alleviates the dilemma of being sensitive to losses.

\begin{table}[t]
  \caption{Ablation experiments with two optimization methods, variable BatchNorm and SR-STE, on the BSD68 dataset.}
  \centering
  \resizebox{\linewidth}{!}{
  \begin{tabular}{ccccc}
    \hline
    \multirow{2}*{{Method}}&\multicolumn{4}{c}{{BSD68}}\\
    \cline{2-5}
    &{$\sigma=15$} &{$\sigma=25$} &{$\sigma=50$} & {FLOPs}\\
    \hline
    vBN $\rightarrow$ BN& 30.02& 28.14& 25.58& 59.69G (68.4\%)\\
    SR-STE $\rightarrow$ STE&31.40& 28.91& 25.94&62.23G (71.3\%)\\
    Ours  &31.59&29.17&26.26&50.88G (58.3\%)\\
    \hline
  \end{tabular}
  }
  \label{tab:tricks}
\end{table}

\noindent
\textbf{Variable BN \& SR-STE.}
As mentioned earlier, we employ variable BatchNorm and SR-STE, and here we conducted ablation experiments to investigate their effects as shown in Table~\ref{tab:tricks}.
The performance drops severely with the same BatchNorm layers for all sparse networks.
We believe this is because pruning leads to different statistics of features and variable BatchNorm provides a solution.
Also, STE is inferior to SR-STE, which only updates non-zero weights that would harm other sparse networks.

\begin{figure*}[t]
  \centering
  \includegraphics[width=1.0\textwidth]{Figs/prob2.pdf}
  \centering
  \caption{Training probability curves comparison of ClassPruning with various losses on the BSD68 dataset.}
  \label{Fig.Prob}
\end{figure*}

\noindent
\textbf{Loss Function.}
We further investigate the effect of the cost loss $L_C$ and entropy loss $L_E$ by removing them from the loss function.
For simplicity, we present the restoration difficulty classification output by the Classifier Module under different losses on BSD68 in Figure~\ref{Fig.Prob}.
We can see that, without $L_C$, the Classifier Module converges quickly; however, it assigns all degraded images to the hard class, which is a bad local minimum for optimization.
And when we remove the $L_E$, the Classifier Module still classifies all inputs into one class, but the classified class is sensitive to the weight ratios of $L_W$ and $L_C$.
Here we present the probability curve by directly removing the $L_E$ from the baseline.
We also only use $L_W$ to train our model and find that the probability curve is similar to the loss without $L_C$ but converges more slowly, indicating our proposed $L_E$ could help the model converge faster.
And only when the three losses are used together, the Classifier Module can classify input images appropriately in a more stable way.

\begin{table}[t]
  \caption{Ablation experiments of performance and computational cost with different network widths on the denoising task.}
  \centering
  \resizebox{\linewidth}{!}{
  \huge
  \begin{tabular}{lcccccc}
    \hline
    \multirow{2}*{{Width}}&\multirow{2}*{{Model}}&\multicolumn{2}{c}{{BSD68}}&&\multicolumn{2}{c}{{SIDD}}\\
    \cline{3-4}
    \cline{6-7}
    & &{PSNR} & {FLOPs} & &{PSNR} &{FLOPs}\\
    \hline
    \multirow{2}*{16}&{DnCNN} & 27.77 & 5.520G(100\%)  &  &37.08 &5.580G(100\%)\\
    &{CP-DnCNN}& 27.68 & 3.320G (60.2\%)& & 35.62&3.930G (69.1\%)\\
    \hline
    \multirow{2}*{32}&{DnCNN} & 28.31 & 21.90G (100\%) & &37.54 &22.04G (100\%)\\
    &{CP-DnCNN}& 28.24& 13.01G (59.4\%)& &37.06  &14.37G (63.8\%)\\
    \hline
    \multirow{2}*{64}&{DnCNN} & 29.01 & 87.28G (100\%)  & & 38.43  &87.58G (100\%)\\
    &{CP-DnCNN}& 29.00 & 50.88G (58.3\%)&  & 38.38  &53.33G (60.9\%)\\
    \hline
    \multirow{2}*{128}&{DnCNN} & 29.32& 348.5G (100\%) & &38.77  &349.1G (100\%)\\
    &{CP-DnCNN}& 29.34& 202.5G (58.1\%)& &38.78  &208.8G (59.8\%)\\
    \hline
    \multirow{2}*{256}&{DnCNN} & 29.52& 1,393G (100\%) & &38.90  &1,394G (100\%)\\
    &{CP-DnCNN}& 29.55& 806.4G (57.9\%)& &38.87 &812.7G (58.3\%)\\
    \hline
  \end{tabular}
  }
  \label{tab:channelAblation}
\end{table}

\noindent
\textbf{Network Width.}
Considering that we perform $N$:$M$ pruning on the convolution kernel to obtain the different capacity sparse networks, it is essential to evaluate the impact of the number of the convolution kernel \emph{i.e.}, network channel.
We conducted the network width ablation study on the denoising task, and Table~\ref{tab:channelAblation} shows the performance and computational cost comparisons of models with different widths.
When the network width is large enough to be over-parameterized, our method can filter the network redundancy sufficiently, reducing more than 40\% computational cost while maintaining or exceeding the original DnCNN's performance.
Such a property is quite valuable since many models nowadays are incredibly complex, and our method can effectively speed them up.

\begin{table}[t]
  \caption{Ablation experiments with different restoration difficulty classifications on the BSD68 dataset. We use \textbf{1, 2, 4} to represent the $1$:$4$, $2$:$4$ and $4$:$4$ fine-grained structured pruning respectively.}
  \centering
  \resizebox{\linewidth}{!}{
  \normalsize
  \begin{tabular}{lccccc}
    \hline
      \multirow{2}*{{Type}}&\multicolumn{5}{c}{{BSD68}}\\
      \cline{2-6}
      &{$\sigma=15$} & {$\sigma=25$} & {$\sigma=50$}& Avg& {FLOPs}\\
    \hline
    \textbf{1} & 31.21& 28.69& 25.83& 28.58 &21.82G (25.0\%)\\
    \textbf{2} & 31.43& 28.85& 25.98& 28.75 &43.64G (50.0\%)\\
    \textbf{4} & 31.63& 29.18& 26.24& 29.02 &87.28G (100\%)\\
    \textbf{1\&2} & 31.28 & 28.74 & 25.88 &28.63&29.33G (33.6\%)\\
    \textbf{1\&4} & 31.39 & 28.89 & 26.00 & 28.76 &49.84G (57.1\%)\\
    \textbf{2\&4} & 31.62 & 29.18 & 26.26 &29.02&65.54G (75.1\%)\\
    \textbf{1\&2\&4} &31.59 &29.17 &26.26 &29.01&50.88G (58.3\%)\\
    \bottomrule
  \end{tabular}
  }
  \label{tab:numberAblation}
\end{table}

\noindent
\textbf{Classification Number.}
An essential concern in ClassPruning is choosing the appropriate classification number to optimally achieve a balance between performance and efficiency.
And we conducted experiments to compare several different classification types.
Here we name the \textbf{1, 2, 4} to represent the $1$:$4$, $2$:$4$ and $4$:$4$ fine-grained structured pruning, and the results are shown in table~\ref{tab:numberAblation}.
We find that when the pruning ratio is relatively high or only one kind of pruning exists, the restoration models' performance drops severely (\emph{e.g.}, \textbf{1}, \textbf{2} and \textbf{1\&2}).
However, when we use \textbf{2\&4} as our classification type, the model performs well, but the computational cost reduction is not attractive enough.
Therefore, we choose the classification type, which divides the restoration difficulty into three classes to trade off the computational cost and performance.

\section{Introduction}

\begin{figure}[t]
    \centering
    \includegraphics[width = 0.47\textwidth]{Figs/TestFrontImg.pdf}
    \centering
    \caption{Our ClassPruning can achieve comparable performance and reduce approximately 40\% computational costs under different network widths for various image restoration tasks.}
    \label{Fig.Front}
    \vspace{-0.22cm}
\end{figure}

Image restoration aims to reconstruct high-quality clean images from low-quality inputs corrupted by various degradations.
Due to the ill-posed nature, obtaining the optimal result in the infinite solution space is highly challenging.
Recently, deep neural networks (DNNs) have developed rapidly and gradually become a preferred alternative to conventional restoration methods.
Its strong capability to restore degraded images has been witnessed in various image restoration tasks, including image denoising~\cite{wang2022uformer,zamir2022restormer}, deraining~\cite{purohit2021spatially,zamir2021multi}, deblurring~\cite{zhang2019deep,cho2021rethinking}, dehazing~\cite{qin2020ffa,ye2021perceiving}, etc.

Despite the remarkable success, high-capacity image restoration networks tend to be quite complex, resulting in significant inference costs, especially when used on mobile devices with limited computational resources.
Therefore, in addition to performance, the computational cost is also a crucial enabling factor~\cite{szegedy2016rethinking}.
And researchers have investigate various methods such as knowledge distillation~\cite{hinton2015distilling}, quantization~\cite{jacob2018quantization,zhou2017incremental}, designing efficient model architecture~\cite{howard2017mobilenets} and pruning~\cite{wen2016learning,li2016pruning} to speed up DNNs.
However, these methods usually obtain a static network applied to all input samples, which limits the model representation capability, as the diverse demands for network parameters and capacity from different instances are neglected.
Dynamic networks, as opposed to static ones, are proposed to adapt their structures or parameters to the input during inference to enjoy favorable properties absent in static models.
It uses lightweight networks to deal with easy samples and complex networks to process hard ones, thus reducing the average computational cost.
Generally, dynamic networks build a model ensemble through a cascaded~\cite{viola2004robust} or parallel~\cite{jacobs1991adaptive} structure and selectively activate the models conditioned on the input, allocating network with appropriate computational cost to each sample.
However, the vast network capacity and low inference cost of dynamic networks come at the cost of increased parameters~\cite{kong2021classsr,park2015big}, which makes it challenging to deploy to edge devices.

We note that pruning is a common technique that discards redundant components to obtain a slimmer network with comparable performance~\cite{frankle2018lottery,gao2018dynamic,oh2022attentive}, which can be
categorized into unstructured pruning and structured pruning.
Unstructured pruning removes individual weights at any location and can achieve a high compression ratio~\cite{han2015learning,guo2016dynamic}.
However, the irregular sparsity patterns of weight tensors lead to poor utilization of current vector processing architectures~\cite{wen2016learning}.
By contrast, structured pruning removes predetermined structures, such as block~\cite{wang2019non} and kernel~\cite{tan2020pcnn,wen2016learning}, which can better utilize processor resources but fail to maintain accuracy.
Recently, $N$:$M$ fine-grained structured pruning~\cite{mishra2021accelerating,zhou2021learning} has emerged as a better alternative to combine both the advantage of unstructured pruning and structured pruning.
Here, the $N$:$M$ sparsity means that at most $N$ weights are non-zero for every continuous $M$ weight.
Such constraint allows the pruned network to be accelerated on modern hardware\cite{mishra2021accelerating}.
Moreover, we observe that high pruning ratio networks can still well restore lightly degraded images while highly degraded images require low pruning ratio networks to process, as shown in Figure~\ref{Fig.Observation}.
Therefore, we believe $N$:$M$ fine-grained structured pruning is a sweet solution to marry dynamic networks.

In this paper, we perform dynamic $N$:$M$ structured pruning on networks based on the input, aiming to reduce the average computational cost while introducing negligible parameters and maintaining performance.
We propose a new image restoration pipeline named ClassPruning, which assigns networks with different pruning ratios to images with varying restoration difficulties.
The pipeline consists of two parts: Classifier Module and Restoration Module.
The Classifier Module is a lightweight classification network that classifies the degraded image into a specific class based on its restoration difficulty.
And the Restoration Module is an image restoration network on which we perform dynamic $N$:$M$ structured pruning to obtain the sparse networks with different capacities.
However, training a dynamic pruning network is not a trivial task, so we propose an improved training strategy.
We first introduce entropy loss and cost loss to enable the network to trade off performance and efficiency.
And we also use variable BatchNorm~\cite{ioffe2015batch} and SR-STE~\cite{zhou2021learning} to stabilize training and improve performance.
Finally, we propose a three-stage training scheme: we first pre-train the Classifier Module using pseudo labels, then fix the Classifier Module and only optimize the Restoration Module, and finally finetune the two parts simultaneously end-to-end until convergence.
We conducted experiments on several representative image restoration networks, \emph{i.e.,} DnCNN~\cite{zhang2017beyond} for denoising, RESCAN~\cite{li2018recurrent} for deraining, SRN~\cite{tao2018scale} for deblurring and GCANet~\cite{chen2019gated} for dehazing. As shown in Figure~\ref{Fig.Front}, ClassPruning could help these restoration networks save approximately 40\% of computational cost while maintaining comparable performance.
\section{Method}
\subsection{Observation}

We first illustrate our observations on images with different restoration difficulties.
Specifically, we perform $N$:$M$ pruning~\cite{mishra2021accelerating} on DnCNN~\cite{zhang2017beyond} with different sparsity ratios to obtain varying capacity models for blind Gaussian denoising tasks.
From the results shown in Figure~\ref{Fig.Observation}, we find that images with low noise levels can still be restored well by DnCNN with a 75\% sparsity ratio.
In contrast, only the unpruned model yields promising restored results for images with high noise levels.
Such a finding inspires us to preserve the whole ability of the network and perform $N$:$M$ pruning on the neural network dynamically according to restoration difficulties of input images, thus achieving a satisfying efficiency and performance tradeoff.
Note that such a method introduces only a few parameters so as not to lead to a burdensome model; instead, the $N$:$M$ fine-grained structured pruning makes the model hardware friendly.

\subsection{Network Architecture}
Our proposed ClassPruning is a novel image restoration pipeline, as shown in Figure~\ref{Fig.Arch}, which consists of two parts: Classifier Module and Restoration Module.
The Classifier Module classifies the input into $L$ classes, while the Restoration Module can be an arbitrary image restoration network on which we perform $N$:$M$ fine-grained structured pruning with different sparsity ratios to obtain $L$ branches $\{f_{R}^{i}\}_{i=1}^L$ for dealing with different inputs.
Specifically, we first feed the degraded image $y$ into the Classifier Module to generate the probability vector $[p_1(y),..., p_L(y)]$ of restoration difficulty.
Then we determine which sparse subnetwork to be sampled based on the index of the maximum probability value $i = \arg{\max}_i  p_i(y)$.
Finally, we handle the input using the $i$-th sparse network of the Restoration Module.

\begin{figure}[t]
    \centering
    \includegraphics[width=1.0\linewidth]{Figs/Observation.pdf}
    \centering
    \caption{PSNR and qualitative results of various noise level images obtained by sparse DnCNNs under different $N$:$M$ pruning ratios. Here we set the $N$:$M$ to be $1$:$4$, $2$:$4$ and $4$:$4$, as well as the noise levels $\sigma$ to 15, 25 and 50.
        }
    \vspace{0.5cm}
    \label{Fig.Observation}
\end{figure}

\begin{figure*}[t]
    \centering
    \includegraphics[width=0.9\textwidth]{Figs/Arch.pdf}
    \centering
    \caption{The overview of the proposed ClassPruning, when the number of image restoration difficulty class $M$ = 3. Classifier Module: aims to generate the probability vector, Restoration Module: aims to restore the degraded images with corresponding sparse subnetworks.}
    \label{Fig.Arch}
\end{figure*}

\noindent
\textbf{Classifier Module.}
The Classifier Module aims to identify the restoration difficulty of the input degraded images efficiently.
Thus, we design a lightweight classification network to minimize the additional introduced computational cost and parameters, which contains four convolutional blocks (Conv-BatchNorm~\cite{ioffe2015batch}-ReLU~\cite{glorot2011deep}), an global average pooling layer~\cite{lin2013network}, and a fully-connected layer.
The convolutional blocks are responsible for image local features extraction, the global average pooling layer aggregates global information, and the fully-connected layer outputs the probability vector.
And our experiments show that such a simple structure can yield satisfying results.

\noindent
\textbf{Restoration Module.}
The Restoration Module is designed to restore degraded images by sparse subnetworks with different capabilities based on the image restoration difficulty.
Unlike previous works~\cite{gao2018dynamic,lin2017runtime} of judging the importance of each convolutional kernel and conducting channel-wise pruning, we would like to prune networks at a fine-grained level to keep models' performance.
Based on the prediction of the Classifier Module, we dynamically perform $N$:$M$ structured pruning with different sparsity ratios to obtain varying complexity subnetworks.
We set the candidate pruning ratios to be $2$:$4$ and $1$:$4$, indicating 50\% and 75\% reduction in computational costs, respectively.
The two pruned networks, together with the base restoration network, form the three branches of the Restoration Module.
In this way, we can restore lightly degraded images with high pruning ratio networks and highly degraded images with low pruning ratio networks, thus speeding the image restoration model up while maintaining performance.

\subsection{Training Strategy}

\noindent
\textbf{Classifier Module.}
Most datasets for image restoration do not label the degraded images with restoration difficulty, preventing us from directly training the Classifier Module.
Therefore, we attempt to label the degraded images with the restoration difficulty class.
Considering that the degradation processes of tasks are different, we employ different methods to synthesize degraded images with restoration difficulty classes.
For image deblurring, we average the adjacent images $y_j$ and then determine the restoration difficulty class based on the number of fused images:
\begin{align}
    y_i = \sum_{j=1}^{N_i} y_{j},
\end{align}
where $N_i$ is the number of fused images for the $i$-th class of restoration difficulty.
And for denoising, deraining, and dehazing tasks, we first subtract the clean image $x$ from the degraded image $y$ to obtain the degraded components $R = y - x$.
Then we categorize all the $R$ into three classes according to their variance and map each $R$ to the median values of the three classes to satisfy the fact that the same scenario often has different levels of degradation:
\begin{align}
    R_i = \frac{R \sqrt{Var(R_{i}^{m})}}{\sqrt{Var[R]}},
\end{align}
where $i$ denotes the class and $m$ means the median value.
Finally, we add the new degraded component $R_i$ to the clean image to obtain the labeled degraded image $y_i$:
\begin{align}
    y_i = x + R_i.
\end{align}
So far, we have obtained degraded images with different restoration difficulties, and we can use Cross Entropy (CE) loss to optimize the Classifier Module.
However, apart from the degradation level, the restoration difficulty is also related to the content of the degraded images, so the labels obtained with the above method may be biased.
Therefore, we use an implicit training approach to finetune the Classifier Module.
Specifically, we propose entropy loss $L_E$ and cost loss $L_C$ to train the Classifier Module combined with the Restoration Module.
The entropy loss helps the Classifier to have much higher confidence in class with the maximum probability than others, which can accelerate the model convergence and ensure the stability of training:
\begin{align}
    L_E &= \sum_{i=1}^{L} - p_i(y) log(p_i(y)).
\end{align}
Besides, we aim to trade off the effectiveness and efficiency of restoration models. However, when only $L_E$ and reconstruction loss exist, all inputs will be classified as hard samples even though the Classifier Module attains a good initialization by CE loss.
To avoid such a bad case, we assign an overhead value $C_i$ to each branch and multiply it with the corresponding probability output by the Classifier Module to form the branch overhead.
We sum all the overhead and minimize it to force the input to be classified into the most appropriate restoration difficulty class:
\begin{align}
    L_C &= \sum_{i=1}^{L} p_i(y) C_i.
\end{align}
Through subsequent training, we enable the Classifier Module to adaptively classify the images with different restoration difficulties and achieve satisfactory results.

\noindent
\textbf{Restoration Module.}
Our restoration network contains a base network and its sparse versions obtained by dynamic $N$:$M$ pruning.
During testing, the input image only passes through one branch, which is determined by the index of the maximum probability value.
In contrast, during training, we allow all degraded images with different restoration difficulties to pass through all $L$ networks and adopt the weighted-$L_1$ loss $L_W$ to guide the models to reconstruct clean images, where the loss weights of each network are the probability output by the Classifier Module:
\begin{align}
    L_W = \sum_{i=1}^{L} p_i(y) ||f_R^i(y) - x ||_1.
\end{align}
And we weight $L_W$, $L_E$ and $L_C$ by $\omega_1$, $\omega_2$ and $\omega_3$ to form the final loss $L_F$:
\begin{align}
    L_F = \omega_1 L_W + \omega_2 L_E + \omega_3 L_C.
\end{align}
Also, we use some strategies to stabilize training and improve performance.

\textit{1) Variable BatchNorm.} The Restoration Module handles images of different restoration difficulties using sparse networks of various capacities.
However, when the network consists of BatchNorm layers, using the same BatchNorm layers for all sparse networks leads to severe performance drops due to different pruning ratios altering feature mean and variance.
And we solve this problem by using the variable BatchNorm~\cite{li2021dynamic}, which assigns individual BatchNorm layers for different sparse networks, which are selectively activated depending on the predicted class.

\textit{2) SR-STE.}
The pruning operation is non-differentiable, and researchers usually employ Straight-through Estimator (STE)~\cite{bengio2013estimating} to compensate for sparse neural networks' backpropagation.
However, STE only updates the weights that are not set to zero, which might harm our networks' performance since all our subnetworks are ensembled in one model, and individual weight updates are inappropriate.
So we use SR-STE~\cite{zhou2021learning}, which updates all parameters by adding a regularization term to decrease the magnitude of the pruned weights.
Experiments show that SR-STE can further improve the model performance.

\noindent
\textbf{Training Pipeline.}
The most intuitive training method is to train the Classifier Module and Restoration Module together from scratch.
However, we find this strategy is quite unstable, and even if this model converges, it does not produce promising results.
Thus, we propose to train the model in three stages:
\begin{itemize}
    \item \textbf{Stage 1}: Pretrain the Classifier Module with CE loss using the degraded images with pseudo labels;
    \item \textbf{Stage 2}: Fix the parameters of the Classifier Module and train the Restoration Module using $L_W$;
    \item \textbf{Stage 3}: Relax all parameters and finetune the whole model with the final loss $L_F$.
\end{itemize}
Experimental results show that this scheme can stabilize training and prevent the network from falling into trivial solutions, and somehow boost the model's performance.

\subsection{Implementation Details}
Considering that we obtain the model in three stages, we first train the Classifier Module using the degraded images with pseudo labels.
The mini-batch size is set to 32, and CE Loss is adopted with Adam optimizer~\cite{kingma2014adam}.
And we use the cosine annealing learning strategy~\cite{loshchilov2016sgdr} to adjust the learning rate.
The total epoch is $200$ with the learning rate from $10^{-3}$ to $10^{-6}$.
Then, we fix the Classifier Module and optimize the Restoration Module using $L_W$ with Adam optimizer.
The training epoch is $500$ with the learning rate from $10^{-4}$ to $10^{-6}$.
In the last stage, we relax all the parameters and train the whole model using $L_F$ with Adam optimizer, and the learning rate decays from $10^{-5}$ to $10^{-7}$ for $500$ epochs.
For fair comparisons, we train all the base models for $1,000$ epochs with the learning rate from $10^{-4}$ to $10^{-7}$.
More details on training protocols are presented in the supplementary material.

\section{Related Work}

\noindent
\textbf{Image Restoration.}
In recent years, learning-based methods~\cite{zhang2017beyond,zamir2020learning,zamir2021multi,zhang2020residual,zamir2020cycleisp} are rapidly developing and achieving more promising results in various image restoration tasks.
DnCNN~\cite{zhang2017beyond} pioneered using convolution neural networks to handle different levels of Gaussian noise.
RESCAN~\cite{li2018recurrent} builds a multi-stage architecture to obtain excellent deraining results.
SRN~\cite{tao2018scale} explores a multi-scale network for boosting image deblurring performance.
And GCANet~\cite{chen2019gated} proposed a gated contextual aggregation network to predict the haze-free image.
Despite the remarkable success, most models are complex and require huge computational resources, thus limiting their applications.
So we do not aim to design a complicated network but propose a novel image restoration pipeline that reduces models' computational cost while introducing negligible parameters and maintaining performance via dynamic $N$:$M$ structured pruning.

\noindent
\textbf{Dynamic Networks.}
Compared to static models, which fix computational graphs and parameters once trained, dynamic networks adjust their structures or parameters to different inputs~\cite{wang2021not,yang2020resolution,yang2020mutualnet,li2021dynamic,chen2022arm}.
Due to its property of assigning different capacity networks to varying inputs, researchers have built dynamic networks through dynamic depth~\cite{graves2016adaptive,wang2018skipnet}, dynamic width~\cite{li2021dynamic,yuan2020s2dnas} and dynamic routing~\cite{kong2021classsr,liu2018dynamic} to reduce the models' average computational cost.
Although these methods accelerate models, they often introduce many parameters, making models burdensome and unfriendly to devices with limited memory storage.
Therefore, in this paper, we are looking for a dynamic network that adjusts its network capacity based on the image restoration difficulties, which introduces negligible parameters and is hardware friendly.

\noindent
\textbf{Network Pruning.}
Removing unnecessary components of deep neural networks is a promising direction to compress and accelerate models~\cite{he2017channel,liu2019metapruning,you2019gate,he2021cap}.
In addition to structured pruning~\cite{han2015learning,frankle2018lottery,ramanujan2020s} and unstructured pruning~\cite{he2017channel,li2016pruning,liu2017learning}, recent works~\cite{mishra2021accelerating,zhou2021learning,oh2022attentive} have proposed $N$:$M$ fine-grained structured pruning as a better alternative that has high efficiency and lossless performance as well as is hardware friendly.
Nvidia~\cite{mishra2021accelerating} first presents the design and behavior of Sparse Tensor Cores supported by Ampere GPU architecture to exploit a $2$:$4$ sparsity pattern and describe a simple workflow for training sparse networks to maintain accuracy.
Then, SR-STE~\cite{zhou2021learning} is proposed to be a simple yet universal recipe for training a $N$:$M$ sparse neural network from scratch.
And SLS~\cite{oh2022attentive} introduce a layer-wise $N$:$M$ sparsity search framework to obtain an efficient image restoration network.
Unlike previous works that proposed a training strategy to obtain a high-performance static sparse network, this paper attempts to perform dynamic $N$:$M$ structured pruning on base restoration models to obtain sparse subnetworks with different capabilities for processing images with varying restoration difficulties.


\section{Introduction}
Population-imbalanced ultracold fermionic systems remain an active field of research since past two decades. While it is known that Feshbach-coupled two-species fermionic gases can demonstrate a BCS-BEC crossover as the effective interaction is varied \cite{BCS-BEC1,BCS-BEC2,BCS-BEC3,BCS-BEC4,BCS-BEC5}, the situation changes if the system resembles a spin-polarized state i.e., one species is more populated than the other. In this case, instead of the weakly interacting homogeneous BCS state (that demands both species to be equally populated), one would have novel pairing structures that can accommodate the unpaired fermions.  Several theoretically proposed phases exist for such systems, including (i) Breached pair (BP) \cite{bp_reference1,bp_reference2,bp_reference3,bp_reference4,bp_reference5,bp_reference6,bp_new2} (ii) Phase separation (PS) \cite{ps_reference1,ps_reference2,ps_reference3} and (iii) Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) states \cite{FFLO_reference1,FFLO_reference2,FFLO_reference3}.

The tunability associated with the ultracold atomic systems went a step ahead with the advent of optical lattices \cite{optical_lattice_reference1,optical_lattice_reference2,optical_lattice_reference3,optical_lattice_reference4} as now the geometry as well as the dimension of the lattice could be easily controlled by using the laser beams. As  Fermionic superfluidity could now be realized in the optical lattices \cite{BCS-BEC5,fermionic_superfluid,hopping6}, it  opened up the possibility of studying  the non-BCS exotic pairing phases in the presence of the lattice potentials. For example, it was suggested that the noise correlation of the two component Fermi gas in an optical lattices would show different signatures for BP and FFLO phases \cite{sarma1}. How the phase diagram evolves with a varying filling fraction and a varying strength of lattice parameters has been studied in the context of BP \cite{sarma3, ps3, sarma4}, FFLO \cite{ ps3, sarma4,fflo1, fflo3}, Phase Separation \cite{ps2, ps3, ps4} , and also for specifically 1D systems \cite{ps,ps1}. Other important features, like, the collective modes in FFLO \cite{fflo2} and the stability properties of FFLO \cite{fflo_reference4,fflo_reference5,fflo_reference6} have also been studied in the presence of the optical lattices. However, it is extremely difficult to find direct experimental signatures of  the theoretically predicted exotic phases (except the PS state \cite{ps_reference3}), and a few indirect methods has been proposed in the past \cite{ alt, lus, pec, anna, singh, bp_diagram, kaj, zhou, bp_new2,continuous}.

In this paper, we study the mean-field level fluctuation dynamics of the population-imbalanced fermionic condensate in an optical lattice for 1D, 2D and 3D. We observe that the dynamics is periodic, and the oscillation frequency  changes almost linearly with a varying Feshbach detuning. We find that the information about both the amount of population imbalance and the pairing phase can be extracted from the slope and intercept of this straight line. We propose that this periodic dynamics can be used as an indirect probe to indicate the FFLO state, as it gives a distinct signature as compared to the other phases both in 2D and 3D. Considering the fact that currently several aspects about this FFLO phase is being earnestly pursued theoretically\cite{fflo_new1, fflo_new2, fflo_new3}, and yet, obtaining its experimental signature has proven to be somewhat tricky :  we believe that such an indirect method might turn out to be extremely useful.

The paper is organized as follows. The basic theoretical model and the dynamical equations are described in Sec. \ref{2nd section}. In Sec. \ref{3rd section},  possible pairing phases are discussed and the relevant parameters are set up in different dimensions. The frequency vs. detuning plots are presented in Sec. \ref{4th section}, and it is shown that the oscillation dynamics is periodic. How the slope of the frequency vs. detuning curve varies with the amount of population imbalance is discussed in Sec. \ref{5th section}, and the same for the intercept is done in Sec. \ref{6th section}. In Sec. \ref{7th section}, we provide the analytical reasoning for the numerical results of the previous two sections. The results are summarized  in Sec. \ref{8th section}.

\section{Model Hamiltonian and dynamics}\label{2nd section}
Here we consider a system of the ultracold fermionic atoms in an optical lattice, that can form a bosonic molecule by means of Feshbach resonances. The system can be described by an extended version of the Hubbard Hamiltonian
\begin{equation}
\begin{split}
H= & \sum_{\langle \vec j.\vec l\rangle}(t-\mu_\uparrow\delta_{\vec l,\vec j}) a^\dagger_{\vec j\uparrow} a_{\vec l\uparrow} +\sum_{\langle \vec j.\vec l\rangle}(t-\mu_\downarrow\delta_{\vec l,\vec j}) a^\dagger_{\vec j\downarrow} a_{\vec l\downarrow}\\ & +g_1 \sum_{\vec j}{n_{\vec j\uparrow}n_{\vec j\downarrow}} +g_2\Big( \sum_{\vec j} b^\dagger_{\vec j}  a_{\vec j\downarrow} a_{\vec j\uparrow} +\sum_{\vec j} a^\dagger_{\vec j\uparrow} a^\dagger_{\vec j\downarrow} b_{\vec j}\Big)\\
& +[\epsilon_b-(\mu_\uparrow+\mu_ \downarrow)]\sum_{\vec j}{b^\dagger_{\vec j}b_{\vec j}}
\end{split}
\end{equation}
Here, up-spin $(\uparrow)$  and down-spin $(\downarrow)$  represent either two different atomic species, or two different hyperfine states of the same atom. $\mu_{\uparrow}$ and $\mu_{\downarrow}$ are their respective chemical potentials.  The hopping amplitude  (assumed to be same for both the species) is $t$. Also, $a^\dagger_{\vec j}$ and $a_{\vec j}$ are respectively  the fermionic atom creation and annihilation operators for site $\vec j =(j_x, j_y, j_z)$, and their bosonic counterparts for molecule creation and destruction are given by $b^\dagger_{\vec j}$ and $b_{\vec j}$.
The number operator for up-spin and down-spin fermions at $\vec j$ th site  are  $n_{\Vec{j}\uparrow}$ and  $n_{\Vec{j}\downarrow}$ respectively.
The on-site interaction strength is $g_1$, while $g_2$ is the additional interaction strength of the feshbach variety which couples two fermionic atoms to form a bosonic molecule.The threshold energy of the composite bosonic molecule energy band is $\epsilon_b$, and it is the pivotal adjustable ``detuning" parameter in this model.

This Hamiltonian takes the following form in momentum space :
\begin{equation}
\begin{split}
 H=& \sum_{\vec k}\epsilon_{\vec k\uparrow}a^\dagger_{\vec k\uparrow} a_{\vec k\uparrow}+\sum_{\vec k}\epsilon_{-\vec k\downarrow}a^\dagger_{-\vec k\downarrow} a_{-\vec k\downarrow}
\\ &+g_1\sum_{\vec k,\vec k'} a^\dagger_{\vec k\uparrow} a^\dagger_{-\vec k+\vec q\downarrow} a_{-\vec k'+\vec q \downarrow} a_{\vec k'\uparrow}\\  &+g_2(\sum_{\vec q,\vec k'} b^\dagger_{\vec q} a_{-\vec k'+\vec q \downarrow} a_{\vec k' \uparrow}+\sum_{\vec k,\vec q'}a^\dagger_{\vec k\uparrow} a^\dagger_{-\vec k+\vec q'\downarrow}b_{\vec q'})\\
 & +[\epsilon_b-(\mu_\uparrow+\mu_ \downarrow)]\sum_{\vec q} b^\dagger_{\vec q} b_{\vec q}
 \end{split}
\end{equation}
$\epsilon_{\Vec{k}\uparrow}$ and $\epsilon_{-\Vec{k}\uparrow}$ are the energy-momentum dispersion relation for up-spin and down-spin respectively. Also, $a^\dagger_{\vec k\uparrow}$, $a_{\vec k\uparrow}$ are the creation and annihilation operators corresponding to momentum $k$ for one fermionic species, and $a^\dagger_{\vec k \downarrow}$, $a_{\vec k \downarrow}$ are the corresponding operators for the other fermionic species. The bosonic  molecule creation and annihilation operators in the momentum space is given by  $b^\dagger_{\vec q}$ and $b_{\vec q}$ respectively. Here $\vec q$ is the momentum of the composite bosons, and has  magnitude zero for BCS-type of pairing, and a non-zero magnitude if an exotic pairing structure like FFLO is involved. More details about these pairing phases are  discussed in \ref{different phases}.

We define
\begin{equation}
  O_{\vec k,\vec q}=\langle a_{-\vec k+\vec q\downarrow} a_{\vec k\uparrow}\rangle
\end{equation}

 The equilibrium conditions are given by

 \begin{equation}
 \frac{\partial O^{eq}_{\vec k,\vec q}}{\partial t}=0 \hspace{15 pt} \mbox{and} \hspace{15 pt}  \frac{\partial b^{eq}_{\vec q}}{\partial t}=0
\end{equation}

Now, we consider the intrinsic and spontaneous quantum fluctuations on top of the equilibrium base states. We write, $O_{\vec k,\vec q}=O^{eq}_{\vec k,\vec q}+\Tilde O_{\vec k,\vec q}$ and
$b_{\vec q}=b^{eq}_{\vec q}+\Tilde b_{\vec q}$. Here
$\Tilde O_{\vec k,\vec q}$ is the fluctuation in $O_{\vec k,\vec q}$ and $\Tilde b_{\vec q}$ is the fluctuation in $b_{\vec q}$. The scheme resembles the one used in \cite{continuous} for studying population-imbalanced uniform systems as well as trapped systems under local density approximation. However, we have extended that to incorporate the effects of optical lattices that are more relevant for present-day experiments.

By calculating the commutation relations of $a_{\vec k\uparrow}$, $a_{-\vec k+\vec q\downarrow}$,  $b_{\vec q}$ and $O_{\vec k,\vec q}$ respectively with the Hamiltonian, and imposing the  equilibrium conditions which mentioned as above,  we arrive at two dynamical equations :
\begin{equation}
\begin{split}
\label{commutation of O}
i\hbar\frac{\partial\Tilde O_{\vec k,\vec q}}{\partial t}= & \sum_{\vec k}(\epsilon_{\vec k\uparrow}+\epsilon_{-\vec k+\vec q \downarrow})\Tilde O_{\vec k,\vec q}
  -(g_1\sum_{\vec k}\Tilde O_{\vec k,\vec q}+g_2 \Tilde b_{\vec q})
\end{split}
\end{equation}
\begin{equation}
\label{boson}
i\hbar\frac{\partial \Tilde b_{\vec q}}{\partial t}=g_2\sum_{\vec k}\Tilde O_{\vec k}+(\epsilon_b-(\mu_\uparrow+\mu_\downarrow))\Tilde b_{\vec q}
\end{equation}
By taking the fourier transform of Eqs. (\ref{commutation of O}, \ref{boson}), we obtain
\begin{equation}
\label{momentum space representation of O}
  \sum_{\Vec{k}} \Tilde O_{\vec k,\vec q}(\omega)=\frac{g_2 \Tilde b_{\vec q}(\omega) f_1(\omega)}{1-g_1 f_1(\omega)}
\end{equation}
Here
\begin{equation}
\label{f1_omega}
 f_1(\omega)=\sum_{\Vec{k}}\frac{1}{\epsilon_{\vec k\uparrow}+\epsilon_{-\vec k+\vec q\downarrow}+\hbar \omega}
\end{equation}

By putting   Eq. (\ref{momentum space representation of O}) in Eq. (\ref{boson}) we  get
\begin{equation}
   b_{\vec q}(\omega)\Big( \epsilon_b+\hbar \omega+\frac{g^2_2  f_1(\omega)}{1-g_1 f_1(\omega)}\Big)=0
\end{equation}

Thus, in the Fourier expansion of $b_{\vec q}(t)$, only those $b_{\vec q}(\omega)$ would survive for which
\begin{equation}
\label{freq}
\epsilon_b+\hbar \omega+\frac{g^2_2  f_1(\omega)}{1-g_1 f_1(\omega)}=0
\end{equation}

Therefore, if there is a single solution  of $\omega$ for which Eq. (\ref{freq}) is satisfied, $b_{\vec q}(t)=b_0+b_1 e^{i\omega t}$, and the condensate fraction $|b_{\vec q}(t)|^2$ would show an oscillatory dynamics with frequency $\pm \omega$. If there are two such solutions $\omega_1$ and $\omega_2$ :  $b_{\vec q}(t)=b_0+b_1 e^{i\omega_1 t}+ b_2 e^{i\omega_2 t}$,  the condensate fraction would have three periodic components :   $\pm \omega_1$, $\pm \omega_2$, $\pm (\omega_1-\omega_2)$, and similarly for higher number of allowed $\omega$ values. If no real solution is found for Eq. (\ref{freq}), that  the condensate fraction would have no oscillatory dynamics at all.

As evident from Eq. (\ref{f1_omega}), the quantity $f_1(\omega)$ involves a sum over all paired region in the momentum space and is dependent on (i) the  dimension of the system, (ii) the amount of population imbalance present in the system, and  (iii) the nature of pairing. Thus, the solution for $\omega$ obtained from Eq. (\ref{freq}), that translates to the frequency of oscillation of the condensate fraction, would contain information about these factors as well.

\section {Different variants of the population imbalanced system }\label{3rd section}
\subsection{Dimension of the system}
Most of the ultracold-atomic experiments offer a tremendous flexibility when it comes to the dimension of the system. For example, 1-dimensional (1D) optical lattice can be formed when the potential  is  shallower in one direction but  it is much  deeper  in the other two directions \cite{1D1,1D2,1D3,1D4}.
Converting the summation in Eq. (\ref{f1_omega}) to an integral,
\begin{equation}
\label{f1 1D}
  f_1(\omega)=\frac{L}{2\pi\hbar} \int \frac{dk}{\epsilon_{\vec k\uparrow}+\epsilon_{-\vec k+\vec q\downarrow}+\hbar \omega}
\end{equation}
Here  $\epsilon_{\vec k\uparrow}=2t(\cos{k a})-\mu_{\uparrow}$,
$\epsilon_{-\vec k\downarrow}=2t(\cos{k a})-\mu_{\downarrow}$. $L$ and  $a$ are the quantization length and lattice constant of the optical lattice respectively.
Similarly, 2-dimensional (2D) optical lattice can be formed when the potential is deeper in one direction but  shallower  in the other two directions \cite{2D}.
\begin{equation}
\label{f1 2D}
    f_1(\omega)=\frac{A}{(2\pi\hbar)^2} \int \frac{dk_x dk_y}{\epsilon_{\vec k\uparrow}+\epsilon_{-\vec k+\vec q\downarrow}+\hbar \omega}
\end{equation}

In this case,
 $\epsilon_{\vec k\uparrow}=2t(\cos{k_x a}+\cos{k_y a})-\mu_{\uparrow}$,
 $\epsilon_{-\vec k\downarrow}=2t(\cos{k_x a}+\cos{k_y a})-\mu_{\downarrow}$ and $A$ is the quantization area of the optical lattice.

In 3-dimensional (3D) optical lattices, the potential depths are comparable in three possible directions \cite{3D}.
Here,
\begin{equation}
\label{f1 3D}
    f_1(\omega)=\frac{V}{(2\pi\hbar)^3} \int \frac{dk_x dk_y dk_z}{\epsilon_{\vec k\uparrow}+\epsilon_{-\vec k+\vec q\downarrow}+\hbar \omega}
\end{equation}
In this case,
$\epsilon_{\vec k\uparrow}=2t(\cos k_x a+\cos k_y a+\cos k_z a)-\mu_\uparrow$,
$\epsilon_{\vec k\downarrow}=2t(\cos k_x a+\cos k_y a+\cos k_z a)-\mu_\downarrow$ and $V$ is the quantization volume of the optical lattice. In Eqs. (\ref{f1 2D}) and (\ref{f1 3D}), it is assumed that the lattice constant is of the same value in all the directions.
\subsection{Possible pairing phases}
\label{different phases}
\begin{figure}[t]
    \includegraphics[width=0.6\linewidth]{Fig_1_a.pdf}
    \includegraphics[width=0.7\linewidth]{Fig_1_b.pdf}
   \includegraphics[width=0.6\linewidth]{Fig_1_c.pdf}
    \includegraphics[width=0.6\linewidth]{Fig_1_d.pdf}
    \caption{The structure of different pairing phases in population imbalanced fermionic systems in 1D :   (a) BP1 (b) BP2 (c) PS (d) FFLO }
    \label{all phases}
\end{figure}
If a system has population imbalance between up-spin and down-spin then different exotic phases  arise, because the paring cannot entirely take place in the BCS-way. The phases that we consider here are
(i) Breached-Pair-1 state (BP1) or Sarma phase  (ii) Breached pair-2 state (BP2) (iii) phase separation (PS) , and (iv) FFLO.
\begin{figure}[b]
    \centering
    \includegraphics[width=0.5\linewidth]{BP1_2D.pdf}
    \caption{The momentum-space structure of BP1 phase in 2D, the blue shell contains superfluid and the red region contains free fermi gas}
    \label{BP1 2D}
\end{figure}
\begin{figure}[t]
    \centering
    \includegraphics[width=0.5\linewidth]{BP2_2D.pdf}
    \caption{The momentum-space structure of BP2 phase in 2D, the blue regions contains superfluid and the red region contains  free fermi gas}
    \label{BP2 2D}
\end{figure}

In Fig. \ref{all phases}, the  blue region and the red dashed regions indicate the paired and unpaired region respectively for a one-dimensional Fermi gas, and different types of pairing structures are depicted. Breached pair or Sarma phase  appears at $T=0$, $\Delta\neq 0$, $P\neq 0$ (here, $\Delta$ is the energy gap and $P=\frac{n_{\uparrow}-n_{\downarrow}}{n_{\uparrow}+n_{\downarrow}}$ is the population imbalance). In the BP1 (Fig. \ref{all phases} (a)) structure for a 1-dimensional system,   we get the paired region near the ends of the momentum-space line, and the unpaired free fermi gas remains in the middle. In BP2 (Fig. \ref{all phases} (b)) structure for 1-dimension, we get the paired region  at both ends of the line and also in the middle, and the ``breached" regions with a single species exist in-between.  Phase separation (Fig. \ref{all phases} (c))  is defined as the coexistence of superfluid ($\Delta\neq 0$, $P= 0$) and normal states ($\Delta=0$, $P\neq 0$). In Fig. \ref{all phases} (c), $-k_1$ to $k_1$ is the range of paired region and  $-k_2$ to $k_2$ is the range of the unpaired region in a momentum space.  FFLO (Fig. \ref{all phases} (d))  state is a BCS-like paired state but  with a non-zero total momentum. FFLO state is not translationally invariant like BCS and BP phases. Here, the gap parameter is $\Delta(r)=\Delta e^{2i\vec q .\vec r}$ ($q>0$) where $\Vec{r}$ is the distance from the origin.
\begin{figure}
    \centering
    \includegraphics[width=0.5\linewidth]{Fig_4_a.pdf}
    \includegraphics[width=0.5\linewidth]{Fig_4_b.pdf}
    \caption{(a) Usual BCS type of pairing (b) FFLO type of pairing}
    \label{FFLO 2D}
\end{figure}
Out of these phases, BP1, BP2 and FFLO involve specific structures in momentum-space. Although these phases have been predicted long back, direct experimental signatures of these phases are  difficult to obtain. In contrast, phase separation appears in the real space, and is easier to identify experimentally \cite{ps_reference3}.
\begin{figure}
  \includegraphics[width=0.5\linewidth]{Fig_5_a.pdf}
    \includegraphics[width=0.5\linewidth]{Fig_5_b.pdf}
    \caption{(a) low attraction strength between red and blue atoms (b) moderate attraction strength between red and blue atoms which generate the phase separation state}
    \label{Normal}
\end{figure}
In 2D or 3D, the BP1 (Fig. \ref{BP1 2D}) is a two-shell structure of the free fermi gas and the superfluid in momentum space. The BP2 state (Fig. \ref{BP2 2D}) in 2D and 3D is a three-shell structure of the free fermi gas and the superfluid in a momentum space.  In  Fig. \ref{FFLO 2D}, the nature of FFLO pairing is depicted where there is the non-zero momentum  pairing between red (up-spin) and blue atoms (down-spin) \cite{BCS-FFLO} in 2D. The phase-separated state, on the other hand, originates when the up-spin (red) and down-spin (blue) atoms pair with each other and repel the redundant atoms \cite{ps_diagram}, leading to a phase separation between the paired superfluid fermions and the excess unpaired fermions (as shown in Fig. \ref{Normal}).
The integration limits for solving Eqs. (\ref{f1 1D}), (\ref{f1 2D}) and (\ref{f1 3D}) later have been
chosen in accordance with the pairing structures described here, and is elaborated further in \ref{4th section}.
\subsection{System parameters}
In this subsection, we set up the parameter values for different types of paring structures to numerically solve Eq. (\ref{freq}). We choose the following values for the  optical lattice depth : $V_x=6 E^F_R$ and $V_y=V_z=40E^F_R$ for 1D; $V_x=V_y=6E^F_R$ and $V_z=40E^F_R$ for 2D; and  $V_x=V_y=V_z=6E^F_R$ for 3D. Here $E^F_R$ is the recoil Fermi energy. The wavelength of the laser beam that creates the optical lattice is chosen to be  $\lambda=825$nm as in \cite{hopping9_and_lambda}.
In  $x$ direction the hopping amplitude is given by  \cite{hopping1,hopping2,hopping3,hopping4, hopping5,hopping6,hopping7_and_g1,hopping8,hopping9_and_lambda}
\begin{equation}
  \frac{t}{E^F_R}=\frac{4}{\sqrt{\pi}}\Big(\frac{V_x}{E^F_R}\Big)^{\frac{3}{4}} e^{-2\sqrt{V_x/E^F_R}}
\end{equation}
 Hopping in $y$ and $z$ directions can be calculated in the similar fashion. The interatomic coupling $g_1$ is given by \cite{hopping7_and_g1}
\begin{equation}
 g_1=\sqrt{\frac{8}{\pi}} k a_{bg} E^F_R \Big(\frac{V_x}{E^F_R}\Big)^\frac{1}{4} \Big(\frac{V_y}{E^F_R}\Big)^\frac{1}{4} \Big(\frac{V_z}{E^F_R}\Big)^\frac{1}{4}
 \end{equation}

The background scattering length is $a_{bg}$ for different fermionic species ($^6\mbox{Li}$ and $^{40}\mbox{K}$) are given in Table \ref{building elements of g1}.
If the optical lattice potential is quite deep ($\geq 5E^F_R$), only  the lowest two energy bands remain relevant. In this case, the additional interaction in the form of Feshbach coupling is given by,
\begin{equation}
\begin{split}
g_2=&\sqrt{\frac{4\pi\hbar^2 a_{bg}\Delta B\mu_{co}}{m_F}}
({\int dx \mathcal{W}^B_x[\mathcal{W}^F_x]^2} {\int dy \mathcal{W}^B_y[\mathcal{W}^F_y]^2}\\ & {\int dz \mathcal{W}^B_z[\mathcal{W}^F_z]^2})
\end{split}
\end{equation}
 Here $\mathcal{W}^B_x$, $\mathcal{W}^B_y$, $\mathcal{W}^B_z$ are Wannier functions for bosons along x, y, z direction respectively, and $\mathcal{W}^F_x$, $\mathcal{W}^F_y$, $\mathcal{W}^F_z$ are their fermionic counterparts \cite{wannier1}. Again, for sufficiently deep lattice potentials, the Wannier functions can be approximated by harmonic oscillator wavefunctions. For  $x$-direction, we thus take the bosonic Wannier function to be
 \begin{equation}
 \mathcal{W}^B_x=\Big(\frac{m_B\omega^B_x}{\pi\hbar}\Big)^{\frac{1}{4}}e^{-(m_B\omega^B_x x^2)/2\hbar}
 \end{equation}
 and fermionic Wannier function to be
 \begin{equation}
 \mathcal{W}^F_x=\Big(\frac{m_F\omega^F_x }{\pi\hbar}\Big)^{\frac{1}{4}}e^{-(m_F\omega^F_x x^2)/2\hbar}
 \end{equation}
 as in \cite{wannier2}.
Here, $m_B$ and $m_F$ are the masses of the bosonic molecule and fermionic atoms respectively. In addition, we consider $\mu_{co}=2\mu_B$ \cite{source_of_table}, $\omega^B_x=\frac{2}{\hbar}\sqrt{V_x E^B_R}$ \cite{omega}, $\omega^F_x=\frac{2}{\hbar}\sqrt{V_x E^F_R}$ \cite{omega}. Here, $E^B_R(=\frac{\hbar^2 k^2}{2 m_B})$ is the bosonic recoil energy and $\Delta B$ is the resonance width. Similarly, $\mathcal{W}^F_y$ and $\mathcal{W}^B_y$ can be computed from $V_y$; while $\mathcal{W}^F_z$ and $\mathcal{W}^B_z$ can be computed from $V_z$.

Combining all these, the coupling $g_2$ can be expressed as :
\begin{equation}
\begin{split}
\label{g2}
   g_2=\sqrt{\frac{4\pi a_{bg}\Delta  B\mu_{co}}{m_F\hbar}} \Big(\frac{2m_B}{\pi}\Big)^{\frac{3}{4}}\Big( E^B_R\Big)^{\frac{3}{8}} \Big(V^B_x V^B_y V^B_z\Big)^\frac{1}{8}\\
\Big(\frac{2}{2+\sqrt{\frac{m_B}{m_F}}}\Big)^{\frac{3}{2}}
\end{split}
\end{equation}
 The complete list of system parameters for 1D, 2D and 3D that we compute, and use for solving the dynamical equations are given in Tables \ref{system parameters for 1D}, \ref{system parameters for 2D} and \ref{system parameters for 3D} in Appendix \ref{appendix A}.

\section{Variation of frequency with detuning}\label{4th section}
In this section, we present the numerical solutions of Eq. (\ref{freq}) in different dimensions, for different types of pairing structures in a population-imbalanced system. The result is presented in the form of oscillation frequency ($\omega$) of the condensate fraction vs. detuning ($\epsilon_b$) of the Feshbach resonance for two distinct fermionic systems : $^6\mbox{Li}$ and $^{40}\mbox{K}$.

\subsection{Dynamics in 1-dimensional systems}
\label{dynamics numerical}

Here the value of $f_1(\omega)$ in Eq. (\ref{freq}) has to be evaluated by using Eq. (\ref{f1 1D}). The limit of the integration is chosen in accordance with the structure of the pairing phases as shown in Fig. \ref{all phases}. As an example, for fractional imbalance, $P=0.5$, we put the limits of the $k$ value of the  paired region  as $-\frac{\pi}{2}$ to $-\frac{\pi}{3}$ and $\frac{\pi}{3}$ to $\frac{\pi}{2}$ for BP1. The paired region in BP2 is taken as $k=-\frac{\pi}{2}$ to $k=-\frac{5\pi}{12}$, $k=-\frac{\pi}{12}$ to $k=\frac{\pi}{12}$ and $k=\frac{5\pi}{12}$ to $k=\frac{\pi}{2}$. The paired region in PS spans from $k=-\frac{\pi}{6}$ to $k=\frac{\pi}{6}$. In FFLO, the paired region  is $k=\frac{\pi}{6}$ to $k=\frac{\pi}{2}$.

The $\omega$ vs. $\epsilon_b$ plots for $^6\mbox{Li}$ near the narrow resonance (resonance width $\Delta B =0.1$ G , resonance position $B_0=543.25$ G) is shown in Fig. \ref{1d,0.5,0.1g}. Here we find that for all values of $\epsilon_b$, there exists a single value of $\omega$. It implies that the dynamics of $b_{\vec q}$ is periodic :
\begin{equation}
\label{periodic}
b_{\vec q}(t)=b_0+b_1 e^{i\omega t }
\end{equation}
Therefore, the condensate fraction goes as
\begin{equation}
\label{periodic_cond}
|b_{\vec q}(t)|^2=|b_0|^2+|b_1|^2 +b_0^\dagger b_1 e^{i\omega t }+b_1^\dagger b_0 e^{-i\omega t}
\end{equation}
So, essentially it would contain a $\cos{\omega t}$
component. Thus, the dynamics of the condensate fraction would be periodic as well.

\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{1d_0_point_5_bp1_0_point_1g_new.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_5_bp2_0_point_1g_new.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_5_ps_0_point_1g_new.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_5_fflo_0_point_1g_new.pdf}
\caption{$\omega$ vs. $\epsilon_b$ plot at fractional imbalance $P=0.5$ for various phases in 1D : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^6\mbox{Li}$ near narrow resonance. Inset : the enlarged view of the no oscillation region indicated by the gap near the resonance point.}
\label{1d,0.5,0.1g}
\end{center}
\end{figure}

The $\omega$ vs. $\epsilon_b$ curve is almost linear in nature. There is a small region with no real solution for $\omega$ in the central region ( Inset of Fig. \ref{1d,0.5,0.1g}).  The plot for $^{40}\mbox{K}$ (resonance width $\Delta B= 7.8$ G, $B_0=202.10$ G) shows a similar pattern in Fig. \ref{1d,0.5,7.8g}, but the region with no solution for $\omega$ is wider. This suggests that the condensate fraction shows no oscillation if the magnitude of the detuning is lower than a threshold value, and there is a periodic oscillation  beyond it.
\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{1d_0_point_5_bp1_7_point_8g.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_5_bp2_7_point_8g.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_5_ps_7_point_8g.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_5_fflo_7_point_8g.pdf}
\end{center}
\caption{$\omega$ vs. $\epsilon_b$ plot at fractional imbalance $P=0.5$ for various phases : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^{40}\mbox{K}$ in narrow resonance width at 1D}
\label{1d,0.5,7.8g}
\end{figure}
We would like to mention that the $\omega$ vs. $\epsilon_b$ curve for the broad resonance in $^6\mbox{Li}$ ($\Delta B= 300$ G, $B_0=834.15$ G) shows a similar behavior. However, the range over which no real solution of $\omega$ is obtained is much wider in this case, and one obtains the straight-line like curve only when the magnitude of the detuning is very large. For the relatively broader resonance of $^{40}\mbox{K}$ ($\Delta B = 9.7 $ G, $B_0=224.21$ G), the $\omega$ vs. $\epsilon_b$ curve is almost similar to the narrower resonance case. These sets of plots (broad resonances of $^6\mbox{Li}$ and $^{40}\mbox{K}$) are not included in the present manuscript.

One point to note is that, in all these plots, no remarkable difference is spotted among the $\omega$ vs. $\epsilon_b$ curves belonging to the different pairing structures. Thus, these cannot directly be used to differentiate among possible exotic phases.

\subsection{Dynamics in 2-dimensional systems }

Here, the limits of the integration  in Eq. (\ref{f1 2D}) is chosen keeping in mind that the momentum space phase diagram in 2D is a circle. We set the total radius of the circle as unity, i.e., all other momenta are scaled by the Fermi momentum $k_F$ of the majority species. For  fractional imbalance $P=0.5$, we put the limits of the  paired region  as $k=\sqrt{2/3}$ to $1$ for BP1. The paired region in BP2 spans from $k=0$ to $k=\sqrt{1/6}$, and $k=\sqrt{5/6}$ to $1$. The paired region in PS extends from $k=0$ to $k=\sqrt{1/3}$. The paired region in FFLO is $k=1-\sqrt{1/3}$ to $1$.

The $\omega$ vs. $\epsilon_b$ plots for $^6\mbox{Li}$ near the narrow resonance  is shown in Fig. \ref{2d,0.5,0.1g}. Here, too,  we find that for all values of $\epsilon_b$, there exists a single value of $\omega$, meaning the dynamics of the condensate fraction is periodic as in Eq. (\ref{periodic}).

\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{2d_0_point_5_bp1_0_point_1g_new.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_5_bp2_0_point_1g_new.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_5_ps_0_point_1g_new.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_5_fflo_0_point_1g_new.pdf}
\caption{$\omega$ vs. $\epsilon_b$ plot at fractional imbalance $P=0.5$ for various phases in 2D : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^6\mbox{Li}$ near narrow resonance.  Inset : the enlarged view of the no oscillation region, indicated by the gap near the resonance point of negative detuning side.}
\label{2d,0.5,0.1g}
\end{center}
\end{figure}

The $\omega$ vs. $\epsilon_b$ curve is almost linear in nature, and again, there is a region with no real solution near the resonance point (Inset of Fig. \ref{2d,0.5,0.1g}) of negative detuning side.  The plot for $^{40}\mbox{K}$  shows a similar pattern in Fig. \ref{2d,0.5,7.8g}. This no-solution region appears mostly in the negative detuning side in case of $^6\mbox{Li}$, and along  both the detuning sides in case of $^{40}\mbox{K}$. However, the region of no solution near the $\epsilon_b=0$ point is less wider here than the analogous region in the 1D case. Though the curve is nearly linear, there is a very small curvature near the resonance.
\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{2d_0_point_5_bp1_7_point_8g.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_5_bp2_7_point_8g.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_5_ps_7_point_8g.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_5_fflo_7_point_8g.pdf}
\caption{$\omega$ vs. $\epsilon_b$ plot at fractional imbalance P=0.5 for various phases : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^{40}\mbox{K}$ in narrow resonance width at 2D}
\label{2d,0.5,7.8g}
\end{center}
\end{figure}

In all these plots, too, the different pairing structures cannot be distinguished solely from the $\omega$ vs. $\epsilon_b$ curves.

\subsection{Dynamics in 3-dimensional systems.}
The $\omega$ vs. $\epsilon_b$ plots for $^6\mbox{Li}$ (narrow resonance)  and $^{40}\mbox{K}$ (narrow resonance) are shown in Fig. \ref{3d,0.5,0.1g} and Fig. \ref{3d,0.5,7.8g} respectively. The momentum space phase diagram in 3D is  a sphere, and we set the total radius of the sphere as unity. As for the limits of the integration in Eq. (\ref{f1 3D}) for fractional imbalance $P=0.5$, the paired region spans from $k=(2/3)^\frac{1}{3}$ to 1 for BP1, from  $k=0$ to $k=(1/6)^\frac{1}{3}$, and $k=(5/6)^\frac{1}{3}$ to $1$ for BP2. The one for  PS is $k=0$ to $k=(1/3)^\frac{1}{3}$ and for FFLO is $k=1-(1/3)^\frac{1}{3}$ to $1$.

Here, too,  we find that for all values of $\epsilon_b$, there exists a single value of $\omega$, meaning the oscillation of the condensate fraction is periodic as in Eq. (\ref{periodic}).
\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{3d_0_point_5_bp1_0_point_1g_new.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_5_bp2_0_point_1g_new.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_5_ps_0_point_1g_new.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_5_fflo_0_point_1g_new.pdf}
\caption{$\omega$ vs. $\epsilon_b$ plot at fractional imbalance P=0.5 for various phases in 3D : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^6\mbox{Li}$ in narrow resonance. Inset : the enlarged view of the no oscillation region, indicated by the gap near the resonance point of negative detuning side.}
\label{3d,0.5,0.1g}
\end{center}
\end{figure}

The same trend (i.e., a periodic oscillation except near the resonance point) is observed in this case as well. The span of the region of no solution  is less than the corresponding regions in 1D or 2D.  The $\omega$ vs. $\epsilon_b$ curve is almost linear.
\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{3d_0_point_5_bp1_7_point_8g.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_5_bp2_7_point_8g.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_5_ps_7_point_8g.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_5_fflo_7_point_8g.pdf}
\caption{$\omega$ vs. $\epsilon_b$ plot at fractional imbalance P=0.5 for various phases : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^{40}\mbox{k}$ in narrow resonance width at 3D}
\label{3d,0.5,7.8g}
\end{center}
\end{figure}

Thus, by comparing  Fig. \ref{1d,0.5,0.1g} , Fig. \ref{2d,0.5,0.1g} and Fig. \ref{3d,0.5,0.1g}, and also by comparing  Fig. \ref{1d,0.5,7.8g}, Fig. \ref{2d,0.5,7.8g} and \ref{3d,0.5,7.8g},  we see that the phase plots do not differ much if the dimension is altered. Also, no major difference vis \`a vis the choice of the atoms ($^6\mbox{Li}$ vs. $^{40}\mbox{k}$) is found.

We would like to highlight here that the frequencies of oscillation of the condensate fraction falls well within experimentally detectable range. If the range of detuning is $5$ G to $50$ G \cite{magnetic_field_for_detuning} away from the resonance point ($B_0=543.25$ G \cite{source_of_table} for narrow resonance width of $^6\mbox{Li}$. ),  we get the frequency in GHz range for all possible dimensions.  As an example, in Fig. \ref{1d,0.5,0.1g}, \ref{1d,0.5,7.8g}, \ref{2d,0.5,0.1g}, \ref{2d,0.5,7.8g}, \ref{3d,0.5,0.1g}, \ref{3d,0.5,7.8g} the $\epsilon_b=1000$ point in our scale implies a detuning of $0.5$ G apart from the resonance point. The corresponding oscillation frequency is $\omega=1000$ in our convention, that translates to a frequency $0.01$ GHz.

Additionally, we observe that  if the  amount of fractional imbalance is changed, then the qualitative nature of the curves remains the same, but some quantitative changes are involved. To quantify these changes, we cast each  $\omega$ versus $\epsilon_b$ plot in the form
\begin{equation}
\label{linearize}
\omega=a\epsilon^2_b+b\epsilon_b+c
\end{equation}

where $a$ is found to be very small ($a\rightarrow 0$). Here $b$ and $c$ are the slope and intercept of the line respectively.

In the next section we  study the effects of population imbalance on the  slope ($b$) of these  $\omega$ vs. $\epsilon_b$ plots.

\section{variation of slope with population imbalance}\label{5th section}

In this section, we study how the $\omega$ vs. $\epsilon_b$ plots in Sec. \ref{4th section} change with a varying amount of population imbalance. For this, we extract the slope $b$ from the linearized form of the curves as in Eq. (\ref{linearize}). The result is presented in the form of slope ($b$)  versus fractional imbalance ($P$)  plot for $^6\mbox{Li}$. We observe that these curves are of parabolic in nature, either concave or convex : depending on the dimensionality and the nature of the pairing structure.

\subsection{Characteristics of 1-dimensional systems}

 Fig. \ref{1d,0.1g,+ve,slope} shows that if $P$  increases then the absolute value of the slope decreases for all four phases in $^6\mbox{Li}$ in 1D for frequency range $16 \leq \omega \leq 1000$. Moreover, all these plots look parabolic, and are concave in nature.
\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{1d_0_point_1g_plus_ve_slope_bp1.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_1g_plus_ve_slope_bp2.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_1g_plus_ve_slope_ps.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_1g_plus_ve_slope_fflo.pdf}
\caption{$b$ vs. $P$ plots for various phases in positive frequency and negative detuning domain : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^6\mbox{Li}$ at narrow resonance width in 1D}
\label{1d,0.1g,+ve,slope}
\end{center}
\end{figure}

\subsection{Characteristics of 2-dimensional systems}
Fig. \ref{2d,0.1g,+ve,slope} corresponds to $b$ vs $P$ values for different phases in population-imbalanced $^6\mbox{Li}$ in 2D for frequency range $9 \leq \omega \leq 1000$, and it shows that if $P$ increases, the magnitude of the slope decreases. The plots of BP1, BP2 and PS phases are parabolic and concave in nature. The plot for FFLO phase, however, offers a striking difference : it is convex.
\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{2d_0_point_1g_plus_ve_slope_bp1.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_1g_plus_ve_slope_bp2.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_1g_plus_ve_slope_ps.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_1g_plus_ve_slope_fflo.pdf}
\caption{$b$ vs. $P$ plot for various phases in positive frequency and negative detuning domain : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^6\mbox{Li}$ at narrow resonance width in 2D}
\label{2d,0.1g,+ve,slope}
\end{center}
\end{figure}
 \subsection{Characteristics of 3-dimensional systems}
Fig. \ref{3d,0.1g,+ve,slope} shows that, like the 2D case, the $b$ vs. $P$ plots for 3D, too, are concave in nature for BP1, BP2 and PS for frequency range $13 \leq \omega \leq 1000$. The corresponding plot for FFLO phase is convex.
\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{3d_0_point_1g_plus_ve_slope_bp1.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_1g_plus_ve_slope_bp2.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_1g_plus_ve_slope_ps.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_1g_plus_ve_slope_fflo.pdf}
\caption{$b$ vs. $P$ plot for various phases in positive frequency and negative detuning domain : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^6\mbox{Li}$ at narrow resonance width in 3D }
\label{3d,0.1g,+ve,slope}
\end{center}
\end{figure}

We would like to point out that Fig. \ref{1d,0.1g,+ve,slope}, Fig. \ref{2d,0.1g,+ve,slope} and Fig. \ref{3d,0.1g,+ve,slope}  correspond to the positive frequency (and negative detuning) sides of Fig. \ref{1d,0.5,0.1g}, Fig. \ref{2d,0.5,0.1g} and Fig. \ref{3d,0.5,0.1g} respectively (and their equivalent plots for other values of $P$). The plots for the negative frequency (and positive detuning) sides are qualitatively similar, and therefore, are not presented separately. The plots for $^{40}\mbox{K}$ near the narrow resonance show the same tendencies, and are not included either.

 Thus, we find that the $b$ vs. $P$ plots for FFLO provide a different signature from the other novel phases in 2D and 3D. This can be used as a tool to detect FFLO, the only phase with non-zero momentum pairing. Moreover, the slope of the $\omega$ versus $\epsilon_b$ curve in any dimension can be used to extract the population-imbalance present in the system by means of Fig. \ref{1d,0.1g,+ve,slope}, Fig. \ref{2d,0.1g,+ve,slope} and Fig. \ref{3d,0.1g,+ve,slope}.

In the next section we try to study the effect of $P$ on the intercept of the linearized form of the $\omega$ vs. $\epsilon_b$ curve.

\section{Variation of intercept with fractional imbalance} \label{6th section}
In this section, the linearized form of the $\omega$ vs. $\epsilon_b$ plot is studied and the intercept is extracted. Here we investigate the effect of the population imbalance on the value of this intercept. The result is presented in the form of intercept ($c$) versus fractional imbalance ($P$) plot for $^6\mbox{Li}$. We find that these plots are of parabolic in nature, and are  either concave or convex.

\subsection{Characteristics in 1-dimensional system}
 Fig. \ref{1d,0.1g,+ve,intercept} and Fig. \ref{1d,0.1g,-ve,intercept} show that if the fractional imbalance increases then the absolute value of the intercept decreases for all the phases in $^6\mbox{Li}$ for frequency range $16 \leq \omega \leq 1000$. The plots are concave in nature for both  positive frequency (negative detuning), and negative frequency  (positive detuning) domains.
\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{1d_0_point_1g_plus_ve_intercept_bp1.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_1g_plus_ve_intercept_bp2.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_1g_plus_ve_intercept_ps.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_1g_plus_ve_intercept_fflo.pdf}
\caption{$c$ vs. $P$ plot for various phases in positive frequency and negative detuning domain : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^6\mbox{Li}$ at narrow resonance width in 1D }
\label{1d,0.1g,+ve,intercept}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{1d_0_point_1g_minus_ve_intercept_bp1.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_1g_minus_ve_intercept_bp2.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_1g_minus_ve_intercept_ps.pdf}
\includegraphics[width=0.45\linewidth]{1d_0_point_1g_minus_ve_intercept_fflo.pdf}
\caption{$c$ vs. $P$ plot for various phases in negative frequency and positive detuning domain : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^6\mbox{Li}$ at narrow resonance width in 1D }
\label{1d,0.1g,-ve,intercept}
\end{center}
\end{figure}

\subsection{Characteristics in 2-dimensional system}
Fig. \ref{2d,0.1g,+ve,intercept} and Fig. \ref{2d,0.1g,-ve,intercept} show that if the fractional imbalance increases then the absolute value of the intercept  decreases for different phases in $^6\mbox{Li}$ for frequency range $9 \leq \omega \leq 1000$. The plots for BP1, BP2 and PS are concave in nature for both positive frequency (negative detuning) and negative frequency (positive detuning) domain. However, the plot corresponding to  the FFLO phase is convex
\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{2d_0_point_1g_plus_ve_intercept_bp1.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_1g_plus_ve_intercept_bp2.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_1g_plus_ve_intercept_ps.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_1g_plus_ve_intercept_fflo.pdf}
\caption{$c$ vs. $P$ plot for various phases in positive frequency and negative detuning domain : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^6\mbox{Li}$ at narrow resonance width in 2D }
\label{2d,0.1g,+ve,intercept}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{2d_0_point_1g_minus_ve_intercept_bp1.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_1g_minus_ve_intercept_bp2.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_1g_minus_ve_intercept_ps.pdf}
\includegraphics[width=0.45\linewidth]{2d_0_point_1g_minus_ve_intercept_fflo.pdf}
\caption{$c$ vs. $P$ plot for various phases in negative frequency and positive detuning domain : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^6\mbox{Li}$ at narrow resonance width in 2D }
\label{2d,0.1g,-ve,intercept}
\end{center}
\end{figure}
\subsection{Characteristics in 3-dimensional system}
 Here we find that (Fig. \ref{3d,0.1g,+ve,intercept} and Fig. \ref{3d,0.1g,-ve,intercept}) the $c$ vs. $P$ plots for BP1, BP2 and PS are concave for both positive frequency (negative detuning) and negative frequency (positive detuning) side for frequency range $13 \leq \omega \leq 1000$.  The plot for FFLO is, on the other hand, convex for  both positive frequency (negative detuning)  and negative frequency (positive detuning) side. \\
\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{3d_0_point_1g_plus_ve_intercept_bp1.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_1g_plus_ve_intercept_bp2.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_1g_plus_ve_intercept_ps.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_1g_plus_ve_intercept_fflo.pdf}
\caption{$c$ vs. $P$ plot for various phases in positive frequency and negative detuning domain : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^6\mbox{Li}$ at narrow resonance width in 3D}
\label{3d,0.1g,+ve,intercept}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=0.45\linewidth]{3d_0_point_1g_minus_ve_intercept_bp1.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_1g_minus_ve_intercept_bp2.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_1g_minus_ve_intercept_ps.pdf}
\includegraphics[width=0.45\linewidth]{3d_0_point_1g_minus_ve_intercept_fflo.pdf}
\caption{$c$ vs. $P$ plot for various phases in negative frequency and positive detuning domain : (a) BP1 (b) BP2 (c) PS (d) FFLO of $^6\mbox{Li}$ at narrow resonance width in 3D}
\label{3d,0.1g,-ve,intercept}
\end{center}
\end{figure}

 Thus, the slope ($b$ vs. $P$) and the intercept  ($c$ vs. $P$)  of the frequency vs. detuning graph lead to a distinct signature for the  FFLO phase in both 2D and 3D.

\section{Analyzing the nature of the slopes and the intercepts } \label{7th section}
In this section, we try to explain the results presented in \ref{5th section} and \ref{6th section}, by means of providing analytical arguments and physical justifications.  We first establish that there would be only one non-trivial solution for $\omega$ from Eq. (\ref{freq}). Then we argue that the absolute value of the slope and the intercept of the $\omega$ vs. $\epsilon_b$ curve increases with a decreasing population imbalance. Finally, we demonstrate how (i) the slope($b$)/intercept($c$) of the $\omega$ vs. $\epsilon_b$ curve changes from one phase to another, and (ii) the $b$ vs. $P$ and $c$ vs. $P$ curves give a distinct signature for the FFLO phase in 2D and 3D.

\subsection{A single branch of frequency }
The function $f_1(\omega)$ can be approximated by a simpler function
\begin{equation}
f_1(\omega)\approx \frac{\alpha}{\omega}.
\label{approx}
\end{equation}

Here the value of $\alpha$ depends on (i) dimension of the system, (ii) amount of population imbalance present, and (iii) nature of the pairing. In Fig. \ref{last diagram} we plot $f_1(\omega)$ and $\alpha$/$\omega$ for 1D to show how close the two functions are. In this plot, the blue solid line and red dashed line represent  $f_1(\omega)$ and $\frac{\alpha}{\omega}$ respectively.  So, Eq. (\ref{freq}) will now be
\begin{equation}
\epsilon_b=-\omega-\frac{g^2_2 \alpha}{\omega-g_1 \alpha}
\end{equation}
Simplifying above equation we  get
\begin{equation}
\label{omega}
\begin{split}
 \omega_{\pm}=\frac{g_1 \alpha-\epsilon_b}{2} \pm\frac{\sqrt{(\epsilon_b-g_1 \alpha)^2-4\alpha(g^2_2-g_1\epsilon_b)}}{2}
 \end{split}
\end{equation}
From the above equation we obtain two roots. One root is close to zero and the other root is prominent. For positive $\epsilon_b$, $\omega_{-}>>\omega_{+}$ and $\omega_{+}\approx 0$. . Similarly, For negative $\epsilon_b$, $\omega_{+}>>\omega_{-}$ and $\omega_{-}\approx 0$.
Thus, only one non-trivial solution for $\omega$ exists, and that explains the nature of the plot in  Fig. \ref{1d,0.5,0.1g}, \ref{1d,0.5,7.8g}, \ref{2d,0.5,0.1g},  \ref{2d,0.5,7.8g}, \ref{3d,0.5,0.1g}, \ref{3d,0.5,7.8g}. No oscillation region of these plots implies imaginary $\omega$. The fluctuation of the condensate fraction is, therefore,  purely periodic as in Eq. (\ref{periodic_cond}).

\subsection{Effect of the population imbalance on the slope and intercept of  frequency vs. detuning curves }\label{slope and intercept}
The limits of the integration in Eqs. (\ref{f1 1D}), (\ref{f1 2D}) and (\ref{f1 3D}) are controlled by the fractional imbalance $P$. So, $f_1(\omega)$ depends on $P$. If $f_1(\omega)$ is approximated as in Eq. (\ref{approx}), we find that if $P$ decreases then $\alpha$ increases and vice versa. A set of  $\alpha$ vs. $P$ values is presented in Table [\ref{variation of f1 with fractional imbalance}] for the PS state, and the overall trend remains the same for the other phases as well.

\subsubsection{Variation of slope  with fractional imbalance}

The slope of the  $\omega$ vs. $\epsilon_b$ curve, as per Eq. (\ref{omega}) is
\begin{equation}
\label{omega vs. epsilon_b}
\begin{split}
 b=\frac{d\omega_{\pm}}{d \epsilon_b}=-\frac{1}{2}  \pm \frac{1}{2}\frac{\epsilon_b+g_1 \alpha}{\sqrt{(\epsilon_b-\alpha g_1)^2-4\alpha (g^2_2-\epsilon_b g_1)}}
 \end{split}
\end{equation}
We consider the non-trivial branch always, i.e., for $\epsilon_b>0$, $b$ is calculated using $\omega_-$; while for $\epsilon_b<0$, $b$ is calculated using $\omega_+$. Accordingly, $b$ vs. $\alpha$ is plotted in Fig. \ref{slope}. It shows that $b$ increases  linearly with $\alpha$ when other parameters are held constant. Thus, if $P$ decreases then $b$  increases. This can explain the traits shown in  Fig. \ref{1d,0.1g,+ve,slope}, \ref{2d,0.1g,+ve,slope},  and \ref{3d,0.1g,+ve,slope}.

  \begin{figure}
      \centering
      \includegraphics[width=0.4\textwidth]{fig_21.pdf}
      \caption{ Plot of $f_1(\omega)$ denoted by blue solid line and its nearest function denoted by red dashed line vs. $\omega$}
      \label{last diagram}
  \end{figure}
  \begin{figure}
      \centering
      \includegraphics[width=0.4\textwidth]{modified_fig_22.pdf}
      \caption{ $b$ vs. $\alpha$ plot for $\epsilon_b<0$}
      \label{slope}
  \end{figure}
 \subsubsection{Variation of intercept  with fractional imbalance }\label{intercept}
 \begin{figure}
     \centering
     \includegraphics[width=0.4\textwidth]{describing_equation3.pdf}
     \caption{$\omega$ vs. $\epsilon_b$ plot for $^6\mbox{Li}$ (narrow resonance width) at $P=0.5$ for $\epsilon_b>0$ of PS state in 1D.}
     \label{+ve intercept}
 \end{figure}
 The $\omega_+$ vs. $\epsilon_b$ curve for $^6\mbox{Li}$ (narrow resonance width) at $P=0.5$ , as obtained from Eq. (\ref{omega}) is presented in Fig. \ref{+ve intercept}, corresponding to  $\epsilon_b>0$ and PS state in 1D.  It clearly indicates a positive intercept, i.e., $c > 0$. It can be argued that if $\omega_{-}\approx 0$ in Eq. (\ref{omega}), then
 \begin{equation}
 \label{eqn.26}
    \omega_{+}\approx -\epsilon_b+g_1 \alpha
\end{equation}
 Thus, $c=g_1\alpha$.
Since  if $P$ decreases then $\alpha$ increases, it implies that  if $P$ decreases then the intercept $c$ increases, and vice versa.
 \subsection{The nature of pairing as reflected in Eq. \ref{f1_omega}}
If we consider different types of pairing (e.g.,  BP1/BP2/PS/FFLO), then  the limits of the integrations in Eqs. (\ref{f1 1D}), (\ref{f1 2D}), (\ref{f1 3D})  differ from one phase to another, affecting the value of $f_1(\omega)$. For example, if we choose a particular fractional imbalance ($P$=0.5, say) then the order  of $f_1(\omega)$ is  PS$>$BP2$>$BP1 for both positive and negative frequency regions and all possible dimensions, as shown in  Table \ref{f1 for positive frequency in 1D}, and \ref{f1 for negative frequency in 1D}. The magnitude of $f_1(\omega)$ for FFLO lies between the $f_1(\omega)$ values of  BP2 and BP1 phase  in 1D. On the other hand, for 2D and 3D, the magnitude of $f_1(\omega)$  for FFLO phase is greater than the magnitudes of $f_1(\omega)$ for all other phases for positive frequency and less than all other phases for negative frequency respectively.  The corresponding numerical values  for a few specific frequency values are presented in Tables \ref{order of f1 for positive frequency in 2D}, \ref{order of f1 for positive frequency in 3D}, \ref{order of f1 for negative frequency in 2D} and \ref{order of f1 for negative frequency in 3D}.
\subsection{Concavity or convexity of slope vs. fractional imbalance plot}
One striking feature of the slope ($b$) vs. population imbalance ($P$) and intercept ($c$) vs. population imbalance ($P$) plots is that the FFLO phase gives a completely distinct signature. In  both the positive and negative frequency domains, the curves are concave for all other phases, viz.,  BP1, BP2 and PS, while for FFLO it is convex.

We try to elaborate on this concavity and convexity. The slope of the $b$ vs $P$ curve is given by
\begin{equation}
  \frac{\partial b}{\partial P}= \frac{\partial b}{\partial \alpha} \frac{\partial \alpha}{\partial P}
\end{equation}
$\frac{\partial b}{\partial \alpha}$ is almost constant for all phases. So, whether the curve will be a concave or a convex one depends solely on $\frac{\partial \alpha}{\partial P}$. For BP1, BP2 and PS phases, the $\alpha$ vs. $P$ curve in the positive frequency domain shows a concave nature (the curve corresponding to the PS phase is shown in Fig. \ref{concave-convex} (a)). In contrast, the curve for FFLO state shows a convex nature (Fig. \ref{concave-convex} (b)).
  \begin{figure}
\centering
\includegraphics[width=0.45\linewidth]{fig_24_a.pdf}
\includegraphics[width=0.45\linewidth]{fig_24_b.pdf}
\centering
\caption{$\alpha$ vs. $P$ plot for various phases in positive frequency and negative detuning domain : (a) BP1/BP2/PS (b) FFLO of $^6\mbox{Li}$ (narrow resonance width) }
\label{concave-convex}
\end{figure}
 Similar arguments can be put forward for the intercept vs. population imbalance curves as well. All these explain the concavity and convexity corresponding to the exotic phases as reported in Sec.\ref{5th section} and Sec.\ref{6th section}.

\section{Discussion} \label{8th section}
In this work, we  theoretically studied the dynamics of the population-imbalanced ultracold fermions in an optical lattice. Considering a Feshbach coupling that enables the fermionic atoms to form bosonic molecules, the fluctuation of the condensate fraction is probed. We showed that beyond a small threshold value of the Feshbach detuning, the dynamics of the condensate fraction is always periodic in nature. This implies  that the oscillation is not just a small ``fluctuation", rather, it is embedded in the mean-field description of the system itself : the periodic dynamics is sustained throughout the course of time. This  result is independent of the amount of imbalance present in the system, as well as of the pairing structure.

In particular, we considered different pairing structures that arise in population-imbalanced fermionic systems, viz., Breached Pair (BP1 and BP2),  Phase Separation (PS), Fulde-Ferrel-Larkin-Ovchinikov(FFLO). Using the exact values of the Feshbach-resonance parameters of $^6\mbox{Li}$ and $^{40}\mbox{K}$, we calculated the oscillation frequencies of the condensate fraction for each of these phases, and in different dimensions. The result is presented in the form of frequency ($\omega$) versus Feshbach detuning ($\epsilon_b$) plots. Again, we observed that neither the nature of the pairing (i.e., whether it is BP1/BP2/PS/FFLO), nor the choice of the samples ( e.g., $^6\mbox{Li}$ or $^{40}\mbox{K}$) / dimensionality (1, 2 or 3 dimension) affects the overall qualitative description : The $\omega$ vs. $\epsilon_b$ curve always shows a straight-line like behaviour. This also implies  that the dynamics of the population-imbalanced system in an optical lattice is markedly different from that in a harmonic trap/ homogeneous system.

However, the equation of this straight line is dependent on the dimension of the system, the nature of the pairing, and most importantly, the amount of imbalance present. Thus, the slope and intercept of this line can be used to characterize the population-imbalanced system : like the amount of imbalance can be deduced for a certain structure from the slope of this line, if the slope/ intercept is known for the same system at some other fixed imbalance values.
Next the slope and intercept of this straight line are plotted with respect to the fractional imbalance.  In this case,  the  FFLO phase gives a distinct signature. In 2D or 3D, the curve for FFLO shows a convex nature (while the curve for all other phases are concave ones).

The FFLO phase is fundamentally different from the other phases in the sense that it hosts a non-zero momentum pairing (and a broken translational symmetry) and that is reflected in this pattern. Thus, if one can construct the  The $\omega$ vs. $\epsilon_b$ lines for a fixed population imbalance by experimentally varying the Feshbach magnetic field, and then repeat the same with different values of imbalance, one can have a wonderful insight in the nature of the pairing. As the detection of exotic phases like FFLO is quite a difficult and elusive process, such an indirect probe might turn out to be of great use in the long run.