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



