\section{Introduction}
\label{sec:intro}
Over the course of the last two decades, there has been a groundswell of theoretical and experimental interest in ultracold gases of alkali metals confined in optical and/or magnetic traps. These systems of ultracold gases have proved to be extremely robust and tunable systems for studying condensed matter physics in regimes that are inaccessible in solid state systems~\cite{lewenstein:review}. During this time, physicists also became specifically interested in condensates of Fermionic alkali atoms (such as $^{6}Li$), and obtaining a \textit{BCS superfluid} of Cooper pairs similar to those seen in solid state superconductors (BCS theory). If the effective attraction between Cooper pairs can be rendered sufficiently strong, Cooper pairs of Fermions are no longer merely correlated and far apart as in traditional BCS systems, but have much smaller correlation lengths approaching the interparticle spacing. Thus, they can be treated as \textit{composite Bosons}, causing the system to undergo a \textit{BCS-BEC crossover} to a BEC superfluid.

Theoretical work on the single channel model of the BCS-BEC crossover over the course of the late $1980$s and $1990$s~\cite{springerlink:10.1007/BF00683774,randeria:fermigas,randeria:bcsbec3,randeria:bcsbec2d,dreschler:bcsbec} motivated numerous experiments with laser cooling and trapping of Fermions, culminating in the observation of this crossover in the early $2000$s~\cite{zwierlein:bcsbecexpt,thomas:superfluid,fermionic:cond:K, collective1}. In the late $1990$s, the more informative \textit{dual-channel} model of \textit{resonant superfluidity} was proposed, by Holland \textit{et al}~\cite{holland:bcsbec}, and by Timmermans \textit{et al}~\cite{timmermans}. This model, built from the \textit{Timmermans' Hamiltonian}, is a generalization of the Dicke (Tavis - Cummings) model in quantum optics~\cite{yuzbashyan}. The attractive interaction between the Fermions can be controlled by tuning a homogeneous magnetic field to their \textit{Feshbach resonances}~\cite{timmermans}, which are caused in a two-particle system  by coupling bound states in a close channel with states in the scattering continuum. The effective scattering length can be tuned simply by varying the external magnetic field that controls the net magnetic moments of the different channels. In this dual channel model, sufficiently strong resonances cause bound states to form in the closed channel. Thus, Cooper pairs of Fermions can \textit{physically combine} into Bosonic molecules~\cite{bcsbec:ohashi}, entities referred to variously as \textit{composite Bosons}, \textit{diatomic molecules}, or \textit{quasimolecular states} in the literature. These composite Bosons have extremely long lifetimes near a Feshbach resonance~\cite{cubizolles:feshbachlife1,strecker:feshbachlife2}, and repel each other~\cite{petrov:repulsivebosons}, facilitating their condensation to a \textit{BEC superfluid}~\cite{donley:becquench,bourdel}. This area of research has gained enormous interest over the last two decades. An overview of the phase portrait of a population balanced Fermi gas can be found in~\cite{randeria:bcsbecnature}, and reviews of the current status of research can be found in~\cite{stringari:review,bloch:review}.

The possibility of observing the nonequilibrium dynamics of coherent quantum states via \textit{quenching} (diabatic variations of the system parameters) is one that is unique to systems of ultracold gases, and unavailable in other similar condensed matter systems. Post-quenched dynamics in cold atom systems have been studied and reported. In particular, the \textit{collapse and revival} of a post-quenched coherent BEC state in an optical lattice~\cite{greiner:colrev} generated considerable interest. Quenched dynamics have also been investigated for Fermi-Bose mixtures, experimentally~\cite{olsen:rabifermi}, as well as theoretically~\cite{andreev:noneqmbcsbec,barankov:bcsbecbloch,uys:hyperbolicbcsbecdyn,jackpu:bcsbecdyn,yuzbashyan,yuzbashyan2}. The possibility of observing a \textit{collapse and revival} phenomenon as the system is quenched past the BCS-BEC crossover has been raised~\cite{huang:becbcs2}. More recently, nonequilibrium dynamics in condensed matter physics have also been described by the complex Ginzburg-Landau equation~\cite{ginzburglandau:review}. The nonlinear damping in the dynamics produces very rich and diverse behavior in Fermi-Bose mixtures~\cite{machida:dynamics}. The mean field dynamics of BCS superconductivity have been obtained from microscopic models via Ginzburg-Landau-Abrikosov-Gor'kov (GLAG) theory. A similar approach has proven successful for BCS-BEC systems as well, thus leading to their description by the TDGL equation for the single channel case~\cite{randeria:bcsbec3}. More recently, the applicability of TDGL dynamics have been demonstrated in the dual channel case~\cite{machida:dynamics}. This motivates the  use of this treatment to study the dynamics of coherent matter waves in BCS-BEC systems in this report.

This paper focuses on the dynamics of the TDGL equation in coupled BCS-BEC systems, and the dynamics of quantum quenching therein. The most general state of this system across the phase diagram is that of a \textit{Fermi-Bose mixture}, where the composite Bosons coexist with correlated Fermions, and their dynamics are linked in the mean field by the two-channel scattering process through coupled Ginzburg-Landau and Gross-Pitaevski equations. The mixture is characterized by two distinct phases, the Fermi (BCS) superfluid phase given by the order parameter $\Psi_1$, and the Bose (BEC) superfluid phase given by the order parameter $\Psi_2$. The Fermi superfluid consists of distinct Cooper pairs, and the Bose superfluid consists of quasimolecular Bosons. This 2-channel model is more faithful to the microscopic nature of the system than the single-channel model, especially when Feshbach resonances are involved. The dynamics of the coupled BCS-BEC phases is richer than that described by a single-phase model in the single channel case, since the Ginzburg-Landau dynamics of the latter is only the damped one except in the BEC regime. This is especially true if the Feshbach resonance is narrow, which is the case being discussed in this work.  {Though many current experiments correspond to broad resonances, the consideration of a narrow Feshbach resonance here is not unrealistic at all}~\cite{rembert:stoof,stoof:other}. The paper reconstructs the time dependent Ginzburg Landau equations for the dynamics of a Fermi-Bose mixture in the two channel case starting from the many-body functional field integral for the Timmermans' Hamiltonian. This Hamiltonian describes the Fermions by the BCS Hamiltonian and the composite Bosons by the Bose-Hubbard Hamiltonian~\cite{timmermans}, and includes a resonant coupling between the two species. The path integral is written as a functional of the superfluid gap parameter $\Delta$ and the boson coherent state (all of them are taken to lie in the zero momentum state) $b_0$, taken to be a c-number. The mean-field dynamics is subsequently derived, and the equilibrium densities evaluated at the stationary solutions. The paper then details investigations of the dynamics of quantum quenching, and reports the possibility of collapse and revival in the full matter wave of the BEC at large times.  Section~\ref{sec:dyneqns} begins by outlying the formalism that obtains the TDGL dynamics from the Timmermans' Hamiltonian. Section~\ref{sec:fpchempot} reports the study of the stationary solutions of the TDGL dynamical system in 3 dimensions. The chemical potential and condensate fractions are evaluated as the Feshbach detuning is varied adiabatically. Section~\ref{sec:quench} looks at the dynamical evolution of the BEC (which can be seen directly in the lab via time-of-flight absorption) as the Feshbach detuning is varied \textit{diabatically} (ie 'quenched') to the shallow BEC regime. Concluding remarks are made in Section~\ref{sec:conclusion}.
\section{Dynamical Equations of Motion}
\label{sec:dyneqns}
The treatment that obtains the coupled TDGL dynamics of this system closely follows that which obtains the conventional TDGL dynamics of the single channel model, as applied by Huang, Yu and Yin~\cite{huang:bcsbecgp}. The formalism is applied to the dual channel case in a manner similar to the treatments by Machida and Koyama~\cite{machida:dynamics}, except for the nature of the Bosonic states, which is approximated by a single mode. It is assumed that, at $T=0$, the time-scales of the dynamics are sufficiently weak so as to not induce Boson excitations above the ground state, an assumption justified in the context of dynamics in greater detail in the literature~\cite{andreev:noneqmbcsbec}. This approximation greatly simplifies the dynamics by removing any spatial information in the composite Boson amplitudes from the beginning. This approximation also allows for rapid transitions near unitarity.

The $D-$dimensional Fermion and zero momentum composite Boson fields in this system are represented by the operators  $\phi_\sigma(x)$ and $b_0$  respectively, with the index $\sigma = \uparrow,\downarrow$ representing the Fermion pseudospin. The dual-channel \textit{Timmermans' Hamiltonian}~\cite{timmermans} $H_{tm}$ for a Fermi-Bose mixture at $T=0$ for a unit volume is
\begin{equation}
 H_{tm} = \int{{\mathrm d}^Dx} \times {\mathcal H_{tm}} (x),
\end{equation}
where
\begin{multline}
\label{timmermans}
{\mathcal H_{tm}}(x) =  \sum_{\sigma} \bigg\{ \phi^\dagger_\sigma(x) \left[h({\bf r})-\mu_F\right]\phi_\sigma(x) \bigg \} - |u_F| \phi^\dagger_\uparrow(x)\phi^\dagger_\downarrow(x)\phi^{}_\downarrow(x)\phi^{}_\uparrow(x)\\
+\left[ 2\nu-\mu_B \right]b^\dagger_0(t')b^{ }_0(t')+ u_B b^\dagger_0b^{ }_0\left(b^\dagger_0b^{ }_0-1 \right)\\
+g_r\left[b^\dagger_0(t') \phi_\uparrow(x) \phi_\downarrow(x) + h.c\right].
\end{multline}
Here, $x=({\bf r},t')$. Equation~\ref{timmermans} describes a system of ultracold electrically neutral two-component Fermions interacting attractively. The first line in equation~\ref{timmermans} represents the \textit{Fermi-BCS} part of the Hamiltonian. Here, $h({\bf r})=\left[ - \frac{\nabla^2}{2m}+V({\bf r}) \right]$ is the single particle Hamiltonian, and $m$ is the Fermion mass (the mass of the composite Bosons is thus $2m$). The second line represents the Hamiltonian of the composite Bosons~\cite{timmermans, huang:becbcs2}. Here, the Feshbach threshold energy (also called the Feshbach 'detuning' from the molecular channel to the continuum~\cite{timmermans}) is represented by $2\nu$, and $u_B$ represents the amplitude of the repulsion between the composite Bosons. The final line describes the Feshbach resonance that leads to the Fermions binding to (or dissociating from) the composite Bosons, with  $g_r$ representing the atom-molecule coupling. Note that the chemical potentials satisfy $\mu_B=2\mu_F$. The path integral grand partition function $Z$ is defined by
\begin{equation}
\label{pathintegral}
 Z=\int{{\mathrm D}[\bar{\phi},\phi]{\mathrm D}[b^*,b]}e^{-S_{\phi b^{ }_0}},
\end{equation}
where ${\mathrm D}[\bar{\phi},\phi]$ and ${\mathrm D}[b^*,b]$ are the path integral measures of the Fermion and Boson fields respectively, and $t'$ is the imaginary time. The action $S_{\phi b^{ }_0}$ is given by
\begin{multline}
\label{pathintegral:action}
 S_{\phi b^{ }_0} = \sum_{\sigma} \int{{\mathrm d}^Dx} \int^\infty_0{{\mathrm d}t'} \bigg[ \bar{\phi}_{\sigma}(x)\partial_{t'} \phi_{\sigma}(x)+ b^*_0(t')\partial_{t'}b^{ }_0(t') +
{\mathcal H_{tm}}(\bar{\phi},\phi,b^*_0,b^{ }_0)\bigg].
\end{multline}
This integral can be evaluated by introducing the macroscopic gap parameter $\Delta(t')$ into equation~\ref{pathintegral} via the Gaussian identity
\begin{equation}
\int {\mathcal D}[\Delta^*,\Delta] \exp{\left[-\int^\infty_0{{\mathrm d}t'}\frac{\Delta^*(t')\Delta(t')}{|u_F|}\right]}=1,
\end{equation}
where the Fermion field (characterized by the gap parameter) is spatially homogeneous due to it's coupling  to the spatially homogeneous Bose field $b_0$ by momentum conservation. Performing the Hubbard-Stratonovich transformation, $\Delta\rightarrow \Delta + |u_F|\phi_\downarrow\phi_\uparrow $ and $\Delta^*\rightarrow \Delta^* + |u_F|\bar{\phi}_\uparrow\bar{\phi}_\downarrow$, cancels out the four-Fermion term from $H_{tm}$ in equation~\ref{pathintegral:action}. Now, integrating out the Fermion fields $\phi,\bar{\phi}$ remaining in equation~\ref{pathintegral} using formal Grassman calculus yields
\begin{equation}
 Z = \int{{\mathrm D}[b^*,b] {\mathrm D}[\Delta^*,\Delta]} e^{-S_{\Delta b^{ }_0}},
\end{equation}
where
\begin{multline}
\nonumber
 S_{\Delta b^{ }_0} = \int{{\mathrm d}t'} \bigg\{ \left[ 2\nu-\mu_B \right] |b^{ }_0(t')|^2 + u_B |b_0(t')|^2\left[|b_0(t')|^2-1 \right] + b^*_0(t')\partial_{t'}b^{ }_0(t')+ \frac{|\Delta(t')|^2}{|u_F|} \bigg\} + \\
\int{{\mathrm d}^Dx \times {\mathrm d}t'} \ln\det{\bf M}(x).
\end{multline}
Here,
\begin{equation}
  {\bf M}(x) \equiv \left[ \begin{array}{cc}
							    \partial_{t'}-\mu_F+h({\bf r})  & \Delta(t')+g_r b^{ }_0(t') \\
							    \Delta^*(t')+g_r b^*_0(t')  & \partial_{t'}+\mu_F-h({\bf r})
                                                             \end{array} \right].
\end{equation}
The action $S_{\Delta b^{ }_0}$ is split into a field independent part $S_0= \int{{\mathrm d}^Dx {\mathrm d}t'} \ln\det{M_0(x)}$, where
\begin{equation}
{\bf M}_0(x)\equiv\left[ \begin{array}{cc}\partial_{t'}-\mu_F+h({\bf r})&0\\
0&\partial_{t'}+\mu_F-{h}({\bf r})\end{array}\right],
\end{equation}
and a field dependent part $S_{eff}$ which vanishes when $\Delta$ and $b_0$ do. Expanding $S_{eff}$ to the fourth order in
$\Delta + g_rb_0$~\cite{huang:bcsbecgp} and performing gradient expansion results in
\begin{multline}
\label{gdexp}
S_{eff}\approx\int{d}^{\mathcal{D}}x\bigg\{\,d\,\left[\Delta^*(t')+g_rb^*_0(t')\right]\partial_{t'}\left[ \Delta(t') + g_rb^{ }_0(t') \right]+
 \frac{|\Delta(t')|^2}{|u_F|}\\
-a|\Delta(t')+g_rb^{ }_0(t')|^2+ \frac{1}{2}b|\Delta(t')+g_rb^{ }_0(t')|^4+ \\
\left[2\nu-\mu_B\right]|b_0(t')|^2  + u_B|b_0(t')|^2\left[|b_0(t')|^2-1 \right]  +b^*_0(t')\partial_{t'}b^{ }_0(t')\bigg\}.
\end{multline}
Here,
\begin{eqnarray}
\label{abc}
a&=& \int{\mathrm d}^{\mathcal{D}}x'\times Q(x-x'/2,x+x'/2),\nonumber\\
b&=&\int{\prod_{i=1}^D{\mathrm d}^{\mathcal{D}}x_i} \times R(x,x_1,x_2,x_3),
\end{eqnarray}
the coefficient $d$ is obtained from~\cite{machida:dynamics},
\begin{equation}
\label{d}
d=\lim_{\omega\rightarrow 0}\int{\mathrm d}^{\mathcal{D}}x'\times {e^{i\omega{t'}}-1\over i\omega} Q(x-x'/2,x+x'/2).
\end{equation}
Here, ${\mathrm d}^{\mathcal{D}}x={\mathrm d}^Dx \times dt'$, where ${\mathrm d}^{D}x$ is the measure of integration over all $D$ spatial degrees of freedom, and
\begin{eqnarray}
Q(x_1,x_2) &=& G_+(x_1,x_2)G_-(x_2,x_1) \nonumber \\
R(x_1,...,x_4)&=&G_+(x_1,x_2)G_-(x_2,x_3)G_+(x_3,x_4)G_-(x_4,x_1),\nonumber \\
\end{eqnarray}
with $x'=({\bf r}',t')$, and the Gor'kov Green's function, $ {\bf G}(x)= {\bf M}^{-1}_0(x)=\bigg[\begin{array}{cc}G_+(x)&0\\
0&G_-(x)\end{array}\bigg]$, which is obtained from the Green's function of a noninteracting Fermi gas i.e.
\begin{equation}
\label{greens:function}
 \left[ \partial_{t'} \mp \mu \pm h({\bf r}) \right] G_\pm(x-x_1) = \delta^D({\bf r}-{\bf r}_1) \delta(t'-t'_1).
\end{equation}
Note that the product in the expression for $b$ above is over the spatial dimensions $D$ only, whereas the measure is over the space-time dimensions $\mathcal{D}$. The mean field equations of motion of the order parameters can be obtained by equating the functional derivatives $\delta S_{eff}/\delta\Delta^*(t')$ and   $\delta S_{eff}/{\delta b^*_0(t)}$ to $0$ after analytically continuing the time to the real axis by substituting $it \rightarrow t$. This yields the final dynamical equations for this system
\begin{eqnarray}
\label{dynamical_system}
\dot \Psi_1 + i\gamma \left( \Psi_1-\Psi_2 \right) - i\alpha \Psi_1 +i\beta |\Psi_1|^2 \Psi_1 &=& 0 \nonumber \\
\dot \Psi_2 +2 i \lambda \Psi_2 +2 i\chi|\Psi_2|^2\Psi_2 - i\kappa\gamma \left(\Psi_1-\Psi_2 \right) &=& 0.
\end{eqnarray}
Here,
\begin{eqnarray}
\label{vartrans}
&&{\Psi}_1 \equiv \frac{\Delta + g_rb_0}{|\mu_F|\sqrt{\mathcal N}}, \nonumber \\
&&{\Psi}_2 \equiv \frac{g_rb_0}{|\mu_F|\sqrt{\mathcal N}},
\end{eqnarray}
where $\mathcal{N}$ is the total (Fermion) particle number. The constants expressed as Greek letters in equation~\ref{dynamical_system} are given by
\begin{equation}
\label{greek_consts}
\begin{array}{ll}
 \alpha\equiv  a|\mu_F|,		&\lambda \equiv \bigg[\nu+|\mu_F|-\frac{u_B}{2} \bigg]|\mu_F|d,  \\
					&					  			 \\
\beta\equiv  b|\mu_F|^3{\mathcal N},    &\sigma \equiv   \frac{|\mu_F|}{g_r},			  	 \\
					&					  			 \\
\gamma\equiv  |\frac{\mu_F}{u_F}|, 	&\chi\equiv  d u_B|\mu_F|\sigma^2{\mathcal N}.			 \\
					&					  			 \\
\kappa  \equiv  g^2_rd,  		& 						  		 \\
\end{array}
\end{equation}
Finally, note that time has been rendered dimensionless via the transformation $t\rightarrow{1\over |\mu_F|d}\times t$, and the chemical potential is presumed to be negative. Thus, the dynamics of the Fermi-Bose mixture is that of a system where the Fermi and Bose fields evolve according to \textit{coupled} Ginzburg-Landau Gross-Pitaevski-Bogoliubov dynamics. The coupling is caused in the Ginzburg-Landau case by the order parameter $\Delta$ getting \textit{nonlinearly dressed} by the Bose field $g_rb_0$, and in the Gross-Pitaevski case by a harmonic coupling to $\Delta$.

It is also noted that the dynamical system bears a resemblance to that of a molecular BEC of atomic Bosons in a resonance effective field theory as proposed by Kokkelmans and Holland in 2002~\cite{kokkelmans:molbecdyn}. The role of the pairing field of noncondensed atoms (represented by the 'anomalous density' of noncondensed pairs) in that system is assumed by the Fermion gap parameter $\Delta$ (related to the anomalous Cooper pair density by $\Delta^* \sim \sum_{\bf k}\langle a_{{\bf k}\uparrow} a_{-{\bf k}\downarrow}\rangle$) in this one. This strengthens the analogy between Cooper pairs in Fermi-Bose systems and noncondensed atoms in Bose-Bose systems.

The dynamics is investigated for a system confined in a 3 dimensional box where the confinement $\mathcal{V}=\mathcal{L}^3$ is much larger than the inter-particle spacing, effectively treating the trap as homogeneous. The constants $a,b$ and $d$ can be evaluated from equations~\ref{abc} and~\ref{d} by using the Green's function for a free particle. This yields~\cite{machida:dynamics,huang:bcsbecgp}
\begin{equation}
\label{Dl0_equation}
\begin{array}{lll}
a=\sum_{|{\bf k}|<k_R}\frac{1}{2\epsilon_{k}}, & b=\sum_{|{\bf k}|<k_R}\frac{1}{4\epsilon_{k}^3},& d=\sum_{|{\bf k}|<k_R}\frac{1}{4\epsilon^2_{\bf k}}.
\end{array}
\end{equation}
Here, $\epsilon_k=\frac{k^2}{2m}-\mu_F$ is the Fermion energy and $k_R={2\pi\over R}$ is the renormalization cutoff in momentum, placed to counter ultraviolet divergences. An important caveat here is that the integrals in equation~\ref{Dl0_equation} contain singularities if $\mu_F\geq\frac{k^2_R}{2m}$. Thus, this formalism breaks down in that regime, which corresponds to regions where the BCS state dominates over the BEC state, well beyond the current region of interest at a negative $\mu_F$. Going to the continuum limit by substituting for the formal sum $\sum \rightarrow \int\mathrm{d}^3k$ in equations~\ref{Dl0_equation} and performing the integrals results in
\begin{eqnarray}
\label{greek_consts_chempot}
\alpha\left(\epsilon_F\right)  &=& \frac{mk_R\epsilon_F}{2\pi^2} - \frac{\left(2m\epsilon_F\right)^{3/2}}{4\pi^2} \arctan{\frac{k_R}{\sqrt{2m\epsilon_F}}},\nonumber \\
\beta\left(\epsilon_F\right)   &=& \frac{\left( {2m\epsilon_F} \right)^{3/2}}{128\pi}\times\mathcal{N}, \nonumber \\
\gamma\left(\epsilon_F\right)  &=&  |\frac{\epsilon_F}{u_F}| =\frac{mk_R\epsilon_F}{2\pi^2}+\frac{\epsilon_F}{|u^0_F|} , \nonumber \\
\kappa\left(\epsilon_F\right)  &=&  \frac{g^2_r}{32\pi}\times \frac{\left(2m\epsilon_F\right)^{3/2}}{\epsilon_F^2},\nonumber \\
\lambda_\nu\left(\epsilon_F\right) &=&  \frac{\nu+\epsilon_F-\frac{1}{2}u_B}{32\pi}\times \frac{\left(2m\epsilon_F\right)^{3/2}}{\epsilon_F}, \nonumber \\
\sigma\left(\epsilon_F\right)  &=&   \frac{\epsilon_F}{g_r},\nonumber \\
\chi\left(\epsilon_F\right)    &=&  \frac{u_B\sigma^2\left(\epsilon_F\right)}{32\pi}\times \frac{\left(2m\epsilon_F\right)^{3/2}}{\epsilon_F}\times\mathcal{N},
\end{eqnarray}
where the dimensionless constants in Greek letters are now functions of the chemical potential $\epsilon_F\equiv|\mu_F|$\footnote{Note that, in general, $\epsilon_F$ as it is defined here does not equal to the Fermi energy.}, and the $\nu$-dependence on $\lambda_\nu$ has been emphasized by subscript. In the equations above, the relation
\begin{equation}
\frac{1}{|u_F|} = \frac{1}{|u^0_F|}+\sum_{|{\bf k}|<k_R}\frac{1}{2\epsilon^0_k},
\end{equation}
has been used to obtain the expression for $\gamma$. Here, $\epsilon^0_k = \frac{k^2}{2m}$, $u^0_F$ is the bare interaction $4\pi a_s/m\mathcal{V}$, and $a_s$ is the s-wave scattering length controlled by Feshbach resonance. This relation comes about as a consequence of renormalizing the BCS gap equation so as to counter ultraviolet divergences in the single channel model~\cite{randeria:fermigas,randeria:bcsbec2d,huang:bcsbecgp}.  { Finally, note that the inability to renormalize this theory analytically disallows taking the theoretical limit of $k_R\rightarrow\infty$ for some of the constants above. This problem will be discussed further in the subsequent section.}
\section{Stationary solutions and Chemical Potential}
\label{sec:fpchempot}
For a complete phenomenological description of the dynamics, number conservation has to be satisfied and used as a constraint at $t=0$ to obtain the chemical potential $\epsilon_F=|\mu_F|$ where $\mu_F = {\mu_B \over 2}$ . In order to do this, the action $S_{eff}$ from equation~\ref{gdexp} is investigated at small temperature $\beta=1/k_BT$. Since we only consider regimes with negative chemical potential, all the Fermions in the gas are correlated and no 'free Fermions' remain~\footnote{The number density of uncorrelated Fermions, obtained from the 'free Fermion' term $S_0= \int{{\mathrm d}^Dx {\mathrm d}t'} \ln\det{M_0(x)}$ in section~\ref{sec:dyneqns}, is $n_{\bf k}=\Theta(\mu_F-\epsilon_{\bf k})$, which is $0$ if $\mu_F<0$.}. Thus, in the stationary case,
\begin{equation}
\label{gdexp:hom}
S_{eff}(\beta)\approx\int^\beta_0{\mathrm d}t'\bigg\{\frac{|\Delta|^2}{|u_F|}-a|\Delta+g_rb^{ }_0|^2 + \left[2\nu-\mu_B\right]|b_0|^2  + u_B|b_0|^2\bigg[|b_0|^2-1 \bigg]\bigg\},
\end{equation}
where the quartic contribution has been neglected. The temperature dependence of all constants can also be neglected and the expression above simplified to get
\begin{equation}
S_{eff}(\beta)\approx \beta \bigg\{\frac{|\Delta|^2}{|u_F|}-a|\Delta+g_rb^{ }_0|^2 + \left[2\nu-\mu_B\right]|b_0|^2  + u_B|b_0|^2\bigg[|b_0|^2-1 \bigg]\bigg\},
\end{equation}
The Helmholtz free energy at $T=0$ is calculated from
\begin{equation}
 \Omega = \lim_{\beta\rightarrow\infty} -\frac{1}{\beta}\ln{Z(\beta)},
\end{equation}
where
\begin{equation}
\label{chempot_partfn}
Z(\beta) = \int D[b^*, b] D[\Delta^*, \Delta] e^{-S_{eff}(\beta)}.
\end{equation}
Simplifying equation~\ref{chempot_partfn} by taking the mean field (ignoring fluctuations) yields
\begin{equation}
\label{chempot_partfn_mf}
Z(\beta) \approx e^{-S_{eff}(\beta)}.
\end{equation}
Thus, the Helmholtz free energy $\Omega$ at $T=0$ is given by
\begin{equation}
\Omega \approx \frac{|\Delta|^2}{|u_F|}-a|\Delta+g_rb^{ }_0|^2 +\left[2\nu-\mu_B\right]|b_0|^2  + u_B|b_0|^2\bigg[|b_0|^2-1 \bigg].
\end{equation}
Imposing number conservation by using $\mathcal{N} = - {\partial\Omega}/{\partial\mu_F}$,
\begin{equation}
\label{eq:numfb}
\mathcal{N} = \frac{\partial a}{\partial\mu} |\Delta + g_rb_0|^2 + 2 |b_0|^2.
\end{equation}
where $\mu=\mu_F=-\epsilon_F$, and the dependence of $a$ on $\mu$ is obtained from equations~\ref{Dl0_equation}. Writing the equation above  in terms of dimensionless variables,
\begin{equation}
\label{chempot}
 2 \sigma^2|\Psi_2|^2 + \xi^2 |\Psi_1|^2 = 1,
\end{equation}
where
\begin{equation}
\label{xisq}
\xi^2 \equiv \mu^2\frac{\partial a}{\partial \mu}.
\end{equation}
Equation~\ref{chempot} is the constraint that fixes the number of particles $\mathcal{N}$ via the chemical potential $\mu_F$. Furthermore, in the stationary case, equations~\ref{dynamical_system} reduce to
\begin{eqnarray}
\label{fpcomplex}
\gamma \left(  \bar{\Psi}_1- \bar{\Psi}_2 \right) - \alpha  \bar{\Psi}_1 +\beta | \bar{\Psi}_1|^2  \bar{\Psi}_1 &=& 0, \nonumber \\
2\lambda_\nu  \bar{\Psi}_2 +2 \chi| \bar{\Psi}_2|^2 \bar{\Psi}_2 - \kappa\gamma \left( \bar{\Psi}_1- \bar{\Psi}_2 \right) &=& 0,
\end{eqnarray}
where $ \bar{\Psi}_1, \bar{\Psi}_2$ are the stationary solutions of the dynamical system. The trivial stationary solutions of the dynamics in the $\Psi_1,\Psi_2$ phase space, $ \bar{\Psi}_1, \bar{\Psi}_2=0$, are realized at temperatures above the critical temperature $T_c$~\cite{huang:bcsbecgp}. The nontrivial ones are the locus of points satisfying the equations above plus the chemical potential equation~\ref{chempot}. Thus, the 3 unknowns $ \bar{\Psi}_1, \bar{\Psi}_2,\epsilon_F$ are solved from the 3 simultaneous equations
\begin{eqnarray}
\label{fixedpoints}
\gamma\left(\epsilon_F\right) \left(  \bar{\Psi}_1- \bar{\Psi}_2 \right) - \alpha\left(\epsilon_F\right)  \bar{\Psi}_1 +\beta\left(\epsilon_F\right) | \bar{\Psi}_1|^2  \bar{\Psi}_1 &=& 0, \nonumber \\
2\lambda_\nu\left(\epsilon_F\right)  \bar{\Psi}_2 +2 \chi\left(\epsilon_F\right)| \bar{\Psi}_2|^2 \bar{\Psi}_2 - \kappa\left(\epsilon_F\right)\gamma\left(\epsilon_F\right) \left( \bar{\Psi}_1- \bar{\Psi}_2 \right) &=& 0, \nonumber \\
2 \sigma^2(\epsilon_F) | \bar{\Psi}_2|^2 + \xi^2(\epsilon_F) | \bar{\Psi}_1|^2 -1 &=& 0.
\end{eqnarray}
As explained in section~\ref{sec:dyneqns}, the BCS-dominant regime is not entirely accessible in this formalism. However, the predominantly BEC gapless regime can be obtained from equations~\ref{fixedpoints}. In that regime, $\Delta=0$ ie $ \bar{\Psi}_1= \bar{\Psi}_2$. In the case of noninteracting Bosons, ie $u_B=0$ ($\chi=0$), the second equation of~\ref{fixedpoints} necessitates that $\lambda_\nu=0$ ie $\epsilon_F=-\nu$, in agreement with mean field results~\cite{huang:becbcs2}. This is also consistent with the physics of the system, since the chemical potential is the energy required to remove one Fermion from the system. In the gapless and noninteracting BEC dominant case, the majority of the Fermions are dimerised with binding energy $-\nu$ , and a dimerised Fermion requires $\nu$ energy to dissociate from the dimer and free itself. For the more general case,
\begin{equation}
\label{bec_fp}
|\bar{\Psi}_{1,2}|^2 = -\frac{\lambda_\nu}{\chi}=\frac{\alpha}{\beta}.
\end{equation}
Clearly, $\lambda_\nu$ and so $\nu$ has to be negative for this to be true ($\chi$ is always positive). Also
\begin{equation}
\label{bec_fp2}
 \frac{\alpha}{\beta}+\frac{\lambda_\nu}{\chi}=0.
\end{equation}
Since $k_R \gg 1$, so is $\alpha$. Therefore, for the equation above to hold, $\lambda_\nu$ and therefore $\nu$ must be large and negative for a gapless BEC dominant state to exist. This criterion is in agreement with the physics of the system. For large negative $\nu$, the binding energy and therefore the molecular affinity of the Boson dimers will be large and negative, facilitating the dimerisation of the majority of the Fermions~\cite{bcsbec:ohashi}.
 In order to calculate the chemical potential in this regime, the relation $\alpha=\beta| \bar{\Psi}_2|^2$ is substituted in equation~\ref{bec_fp2}. Then,  equations~\ref{greek_consts_chempot} are applied, yielding
\begin{equation}
 -\frac{\nu+\epsilon_F-\frac{1}{2}u_B}{u_B} = \frac{g_r |b_0|^2}{\epsilon_F}.
\end{equation}
Solving this quadratic equation in $\epsilon_F$ and rejecting the unphysical root,
\begin{equation}
\label{chempot_purebec}
\epsilon_F = -\frac{1}{2}\left(\nu-\frac{1}{2}u_B\right)\bigg\{1+ \bigg[1-\frac{4g_ru_B|b_0|^2}{\left(\nu-\frac{1}{2}u_B\right)^2} \bigg]^{1/2}\bigg\} \approx - \left(\nu-\frac{1}{2}u_B\right) \bigg[1-\frac{g_ru_B|b_0|^2}{\left(\nu-\frac{1}{2}u_B\right)^2} \bigg].
\end{equation}
In the case of noninteracting Bosons, equation~\ref{chempot_purebec} gives $\epsilon_F=-\nu$ as above. Furthermore, in this regime $|b_0|^2=\mathcal{N}/2$ (no BCS state, all Fermions are dimerised), or $|\bar{\Psi}_{1,2}|^2 = 1/2\sigma^2(\epsilon_F)$. This simplifies equation~\ref{chempot} to $\xi = 0$ or  $\partial a/\partial\mu = 0$. Figure~\ref{fig:xi} shows plots of $\partial a/\partial\mu $ as functions of $\nu$ (where $\nu= \mu_F$) for $m=1$, $\mathcal{V}=1$ and several large values of $k_R$. Note that  $\partial a/\partial\mu $ and therefore $\xi$ indeed vanishes in the limit $\nu \rightarrow -\infty$, which is where the BEC regime is expected.

Continuing with the case of noninteracting Bosons ($u_B=0$ ie $\chi=0$), equations~\ref{fixedpoints} reduce to
\begin{eqnarray}
\label{fixedpoints_noint}
\gamma\left(\epsilon_F\right) \left(  \bar{\Psi}_1- \bar{\Psi}_2 \right) - \alpha\left(\epsilon_F\right)  \bar{\Psi}_1 +\beta\left(\epsilon_F\right) | \bar{\Psi}_1|^2  \bar{\Psi}_1 &=& 0, \nonumber \\
\lambda_\nu\left(\epsilon_F\right)  \bar{\Psi}_2 - \kappa\left(\epsilon_F\right)\gamma\left(\epsilon_F\right) \left( \bar{\Psi}_1- \bar{\Psi}_2 \right) &=& 0, \nonumber \\
2 \sigma^2(\epsilon_F) | \bar{\Psi}_2|^2 + \xi^2(\epsilon_F) | \bar{\Psi}_1|^2 -1 &=& 0.
\end{eqnarray}
with $\lambda_\nu(\epsilon_F)=\frac{(2m\epsilon_F)^{3/2}}{32\pi}[1+\frac{\nu}{\epsilon_F}]$. Solving the first two equations of~\ref{fixedpoints_noint} yields
\begin{eqnarray}
\label{fixedpoints_noint_final}
 \bar{\Psi}_1 &=& \bigg[\frac{1}{\beta} \left(\alpha - \gamma\eta_\nu \right)\bigg]^{1/2} e^{i\theta} \nonumber \\
 \bar{\Psi}_2 &=& \left( 1-\eta_\nu\right)  \bar{\Psi}_1,
\end{eqnarray}
where $\eta_\nu=\frac{\lambda_\nu/\kappa\gamma}{1+\left(\lambda_\nu/\kappa\gamma\right)}$.  {The problem of renormalizability of this dual-channel theory mentioned in the previous section manifests here. The inability to renormalize happens because of the quantity $\alpha-\gamma\eta_\nu$ in the equations above. In the single channel case, when $g_r, \kappa\rightarrow0$ and $\eta_\nu\rightarrow 1$, this expression converges as $k_r\rightarrow \infty$, as can be seen using equations}~\ref{greek_consts_chempot},  {yielding}
\begin{equation}
\lim_{k_R\rightarrow\infty}\left[\lim_{g_r\rightarrow0}\left(\alpha-\eta\gamma\right)\right] = -\left[\frac{\epsilon_F}{|u^0_F|}+\frac{1}{8\pi}(2m\epsilon_F)^{3/2}\right].
\end{equation}
 {However, finite nonzero $g_r$ leads to $\eta$ vanishing as $k_R\rightarrow\infty$, yielding $\lim_{k_R\rightarrow\infty}\left(\alpha-\eta\gamma\right) \sim \alpha$ which diverges for arbitrarily large $k_R$} (see equations~\ref{greek_consts_chempot})\footnote{ {It should be noted, however, that this divergence goes away at $T\neq 0$ due to the use of thermal Green's functions (instead of the ones from equation}~\ref{greens:function})  {in the evaluation of the constants in equations}~\ref{abc} ~\cite{machida:dynamics}~\cite{huang:bcsbecgp}.}.

 {The difficulty described above can be worked around by noting that the allowed momenta need to have a finite cutoff $k_R$, which is governed by the physics of the system. For instance, in the case of solid state BCS systems, the cutoff is given by the Debye frequency on account of the physical origins of the electron-electron attraction in the slowly relaxing lattice vibrations. Ordinarily, the choice of cutoff is decided by the range of the interatomic potential, with $k_R\sim 8\pi^2\nu/g^2_r$ if the bare repulsion between the Fermions is neglected~\cite{timmermans,pethick:bec}. However, for small confinements and small detuning, the maximum momentum allowed in these cold atom systems cannot be greater than $2\pi/R$, where $R$ is the size of the atoms. Experimentally, these atoms are very cold, and never move so fast as to cause inelastic scattering beyond the Feshbach resonance, since such excitations will dynamically alter the internal degrees of freedom of the atoms themselves. Thus, this choice of $k_R$ keeps a reasonable upper bound in all of the integrations. For Fermionic atoms like $Li$, $R$ is ${\mathcal O}(10^{-10}m)$, which needs to be rescaled with respect to the unit length viz. the confinement dimensions. The confinement, however, is highly tunable in cold atom systems, with a fairly wide range of permitted values. Thus, using a wide variety of $k_R$'s to get quantitatively correct results should be permissible, once the trap sizes are adjusted accordingly. In this report,  $k_R$ are chosen to be ${\mathcal O}(10)$ in units of inverse trap size. These should be attainable in small volume systems of ${\mathcal O}(10^2)$ atoms (obtained experimentally via laser culling by Chuu et al}~\cite{raizen:culling})  {confined to length scales of ${\mathcal O}(10^{-9}m)$  (obtained experimentally by confining ultracold gases using atom chips}~\cite{doublewell:chip}).

Equations~\ref{fixedpoints_noint_final} are nothing more than the gap equation for the BCS-BEC system in position space. In all these equations, $\alpha,\beta,\gamma,\kappa,\lambda_\nu,\eta_\nu$ are functions of $\epsilon_F$. Note from the above that if $\kappa\rightarrow0$ ie $g_r \rightarrow0$, then $\eta_\nu\rightarrow 1$, and  $ \bar{\Psi}_2\rightarrow0$, which will give  a pure BCS until $\epsilon_F$ goes to a regime where $\alpha$ can no longer be evaluated without running into singularities (as per section~\ref{sec:dyneqns}), causing this formalism to break down. Plugging the values of $\bar{\Psi}_{1,2}$ above to the final equation of~\ref{fixedpoints_noint} yields
\begin{equation}
\label{chempot_final_noint}
1=\frac{1}{\beta(\epsilon_F)}\bigg[\alpha(\epsilon_F) - \gamma(\epsilon_F)\eta_\nu(\epsilon_F) \bigg] \bigg\{2\sigma^2(\epsilon_F)\big[1-\eta_\nu(\epsilon_F)\big]^2 +\xi^2(\epsilon_F)\bigg\}.
\end{equation}
\pagebreak
\begin{figure}[h!bt]
\ \epsfig{file=fig1.eps,height=3.0in,width=4.5in}
\caption{(Color online) Plots of $\partial a/\partial\mu$ as a function of $|\mu_F|$ for $k_R=100$ (red), $500$ (green), and $10^9$ (blue). Here, $m=1$ and $\mathcal{V}=1$. Notice the rapid rate of convergence as $k_R\rightarrow\infty$.}
\label{fig:xi}
\end{figure}
Equation~\ref{chempot_final_noint} is a transcendental equation and needs to be solved numerically for $\mu_F$ ($-\epsilon_F$).
Figure~\ref{fig:chempot} contains the numerical results of evaluating the condensate fractions and chemical potentials. These have been obtained by setting $\mathcal{N}$ to $100$ and using representative values of the interaction parameters. The chemical potential $\mu_F$ was obtained by numerically solving equation~\ref{chempot_final_noint} using Newton-Raphson methods. The default working precision was kept to the eighth place of decimal. Figures~\ref{fig:chempot} (a), (b) and (c) contain plots for $g_r = 25$, and figures~\ref{fig:chempot} (d), (e) and (f) contain plots for $g_r = 40$. Figures~\ref{fig:chempot}(a) and (d) are plots for $\mu_F/\nu$ as a function of $\nu$. Note from these figures that, for sufficiently large $-\nu$, $\mu_F \approx \nu$ in accordance with the analytical results in the BEC regime. Figures~\ref{fig:chempot}(b) and (e) are plots for the condensate fractions $n_F$ (Fermions in BCS) and $n_B$ (molecular composite Bosons in BEC). The fractions are computed after solving equation~\ref{chempot_final_noint}, and substituting the corresponding values of $|\bar{\Psi}_{1,2}|^2$ (obtained from
equations~\ref{fixedpoints_noint_final}) into the relations (see equation~\ref{chempot})
\begin{eqnarray}
n_F &=& \xi^2 | \bar{\Psi}_1|^2, \nonumber \\
n_B &=& \sigma^2| \bar{\Psi}_2|^2.
\end{eqnarray}
The rise in the BCS superfluid density can be clearly seen in figure~\ref{fig:chempot}(b) around $\nu \approx -50$, and in figure~\ref{fig:chempot}(e) of $\nu \approx -150$. These results agree qualitatively with more sophisticated theories of BCS-BEC systems~\cite{leggett:book}, considering that numbers are highly sensitive to the choice of the renormalization cutoff $k_R$.
In the case of interacting bosons, equations~\ref{fixedpoints} are simplified by operating in a regime where the chemical potential is weak enough so that $\beta$ can be ignored. This yields
\begin{eqnarray}
\label{fixedpoints_nobeta}
\gamma\left(\epsilon_F\right) \left(  \bar{\Psi}_1- \bar{\Psi}_2 \right) - \alpha\left(\epsilon_F\right)  \bar{\Psi}_1 &=& 0, \nonumber \\
2\lambda_\nu\left(\epsilon_F\right)  \bar{\Psi}_2 +2 \chi\left(\epsilon_F\right)| \bar{\Psi}_2|^2 \bar{\Psi}_2 - \kappa\left(\epsilon_F\right)\gamma\left(\epsilon_F\right) \left( \bar{\Psi}_1- \bar{\Psi}_2 \right) &=& 0, \nonumber \\
2 \sigma^2(\epsilon_F) | \bar{\Psi}_2|^2 + \xi^2(\epsilon_F) | \bar{\Psi}_1|^2 -1 &=& 0.
\end{eqnarray}
Solving the first two equations of~\ref{fixedpoints_nobeta} yields the gap equations in this regime. The solutions are
\begin{eqnarray}
\label{fixedpoints_nobeta_final}
\bar{\Psi}_1 &=& \frac{\alpha}{\gamma-\alpha} \bar{\Psi}_2, \nonumber \\
\bar{\Psi}_2 &=& \bigg \{-\frac{\lambda_\nu}{\chi} \bigg[1+\frac{\kappa\gamma}{2\lambda_\nu\left(\gamma-\alpha \right)} \bigg]  \bigg \}^{1/2} e^{i\theta},
\end{eqnarray}
where a necessary condition is that $\lambda_\nu$ be negative. Substituting these results into the final equation of~\ref{fixedpoints_nobeta} yields the transcendental equation for $\epsilon_F$ in this case,
\begin{equation}
\label{chempot_final_nobeta}
1+\frac{2\sigma^2\left(\epsilon_F\right)\lambda_\nu\left(\epsilon_F\right)}{\chi\left(\epsilon_F\right)} \bigg\{1+\frac{\xi^2\left(\epsilon_F\right)\gamma^2\left(\epsilon_F\right)}{2\sigma^2\left(\epsilon_F\right)\left[\gamma\left(\epsilon_F\right)-\alpha\left(\epsilon_F\right) \right]^2} \bigg\} \bigg\{1+\frac{\kappa\left(\epsilon_F\right)\gamma\left(\epsilon_F\right)}{2\lambda_\nu\left(\epsilon_F\right)\left[\gamma\left(\epsilon_F\right)-\alpha\left(\epsilon_F\right) \right]} \bigg\} = 0.
\end{equation}
\pagebreak
\begin{figure}[h!bt]
\hspace*{-0.3in}
\ \epsfig{file=fig2.eps,height=6in,width=6in}
\caption{(Color online) Plots of the  Fermion/Boson number and  chemical potential $\mu_F$  as a function of the Feshbach detuning $\nu$ in the deep-BEC and shallow-BEC regimes. Here, $\hbar = m =\mathcal{V} = 1$. The renormalization cutoff $k_R$ is set to $18$. The curves are evaluated for two values of the atom molecule coupling $g_r$. Figures (a), (b) and (c) contain plots for the ratio $\mu_F/\nu$, the condensate fractions and $\mu_F$ respectively for $g_r=25$, $u_F = -0.3$, $u_B = 0$  and $\mathcal{N} = 100$. Figures (d), (e) and (f) contain plots for the ratio $\mu_F/\nu$, the condensate fractions and $\mu_F$ respectively for $g_r=40$, $u_F = -0.3$, $u_B = 0$  and $\mathcal{N} = 100$. In figures (b) and (e), the Fermion condensate fraction $n_F$ is shown in red, and the Boson molecular condensate fraction $n_B$  is shown in blue. Note that the pure molecular BEC is achieved when $n_B = 0.5$, half the total number $\mathcal{N}$. Also, note the onset of the shallow BEC regime at $\nu \approx -50$ for (b) and at $\nu \approx -150$ for (e). Also, note from figures (a) and (d) that $\mu_F$ approaches $\nu$ for sufficiently large values of $-\nu$, indicating the onset of a pure molecular BEC in accordance with the results in Equation~\ref{chempot_purebec}. Finally, note that figs (c) and (f) show that $\mu_F$ is approximately equal to $\nu$ for large $|\nu|$, but deviates from this at smaller $|\nu|$, as can be seen more clearly in figs (a) and (d).}
\label{fig:chempot}
\end{figure}
\pagebreak
\begin{figure}[h!bt]
\ \epsfig{file=fig3.eps,height=6in,width=4.5in}
\caption{(Color online) Plots of the  Fermion/Boson number and  chemical potential $\mu_F$  as a function of the Feshbach detuning $\nu$ in the deep-BEC and shallow-BEC regimes. Here, $\hbar = m =\mathcal{V} = 1$, $u_F= -0.3$, $u_B = 0.6$ and $\mathcal{N}=100$. The renormalization cutoff $k_R$ is set to $14.5$. The curves are evaluated for $g_r=25.0$. Figure (a) contains plots for the condensate fractions as a function of the Feshbach detuning $\nu$. The Fermion condensate fraction $n_F$ is shown in red, and the Boson molecular condensate fraction $n_B$  is shown in blue. Note that the pure molecular BEC is achieved when $n_B = 0.5$, half the total number $\mathcal{N}$. Also, note the onset of the shallow-BEC regime at $\nu \approx -100$. Figure (b) contains plots of the chemical potential $\mu_F$ as a function of $\nu$. Note here that $\mu_F$ approaches $\nu$ for sufficiently large values of $-\nu$, indicating the onset of a pure molecular BEC in accordance with the results in Equation~\ref{chempot_purebec}.}
\label{fig:chempot2}
\end{figure}
\pagebreak
\begin{figure}[h!bt]
\ \epsfig{file=fig4.eps,height=3.0in,width=4.5in}
\caption{(Color online) Plots of $\Omega_{b_0}$ as a function of chemical potential $\mu_F$ for $g_r=40$, $u_F=-0.3$, $u_B=0$, and $\mathcal{N}=100$. Compare order of magnitudes with that of $\Omega_m$ in Table~\ref{tab:colrev}.}
\label{fig:omega}
\end{figure}
Equation~\ref{chempot_final_nobeta} is solved numerically using the same algorithms and tolerances as equation~\ref{chempot_final_noint} for representative values of the parameters. The results are shown in figures~\ref{fig:chempot2}. Figure~\ref{fig:chempot2} (a) contains plots of the condensate fractions $n_F$ (red) and $n_B$ (blue) for $\mathcal{N}=100$, $g_r = 25$, $|u_F| = 0.3$, and $u_B = 0.6$. Here, $k_R$ is chosen to be $14.5$. Note the shallow BEC point at $\nu \approx -100$ in Figure~\ref{fig:chempot2} (a). In this case, the fermion population appears to increase much faster than in the case of noninteracting Bosons as $\nu$ is varied adiabatically past this regime. The system should populate to a full BCS state shortly after unitarity.
\section{Quenched Time Evolution}
\label{sec:quench}
Equations~\ref{dynamical_system} provide a complete description of mean field dynamics of the Timmermans' Hamiltonian in equation~\ref{timmermans}. The $\Psi_i$s will not, in general, represent densities when the system departs from equilibrium. The definition of $\Psi_2$ in equation~\ref{vartrans} indicates that  while $|\Psi_2|^2(t)$ will be proportional to the Boson density, $|\Psi_1|^2(t)$ will no longer represent the Fermion density out of equilibrium (which will only depend on $1-2\sigma^2|\Psi_2|^2$), since equation~\ref{eq:numfb} breaks down away from equilibrium. The general quasiharmonic solutions of $\Psi_{1,2}$ have been discussed by Machida and Koyama~\cite{machida:dynamics}, and the global existence of weak solutions have been established recently by Chen \textit{et al}~\cite{chen:tdglsoln}. This paper focuses on the dynamics in the shallow BEC regime close to the BCS-BEC crossover as the system is quenched from the deep-BEC regime.

The continued use of the single mode approximation in this case can be justified by evaluating the maximum characteristic time scale $t_q$ for the quench below which sizable excitations of the closed channel Bosons and free Fermions will occur. In order to do so, the Timmermans' Hamiltonian in equation~\ref{timmermans} is written in momentum space with the closed channel excitations included, yielding
\begin{multline}
\label{eq:timmermans:kspace}
H_{tm} = \sum_{{\bf k},\sigma} \left(\epsilon_{\bf k} - \mu_F \right)a^\dagger_{{\bf k}\sigma}a^{ }_{{\bf k}\sigma} - |u_F|\sum_{{\bf k}{\bf k'}}a^\dagger_{{\bf k}\uparrow}a^\dagger_{-{\bf k}\downarrow}a^{ }_{{\bf k'}\downarrow}a^{ }_{{\bf k'}\uparrow}\\
+\sum_{\bf q} \left(E^0_{\bf q}+2\nu-\mu_B \right)b^\dagger_{\bf q}b^{ }_{\bf q} + u_B\sum_{{\bf q_1}{\bf q'_1} {\bf q_2} {\bf q'_2}}b^\dagger_{\bf q'_1}b^\dagger_{\bf q'_2} b_{\bf q_2} b_{\bf q_1}\delta_{{\bf q_1}+{\bf q_2},{\bf q'_1}+{\bf q'_2}}\\
+g_r\sum_{{\bf k}{\bf q}}\left(b^\dagger_{\bf q}a_{{\bf p+\frac{q}{2}}\uparrow}a_{{\bf -p+\frac{q}{2}}\downarrow} + h.c. \right).
\end{multline}
Here, the first line represents the Fermions in momentum space, with the creation (annihilation) operators $a^\dagger_{{\bf p}\sigma}$ ($a^{ }_{{\bf p}\sigma}$) obtained by Fourier transforms of $\phi^\dagger_\sigma$ ($\phi^{ }_\sigma$) from equation~\ref{timmermans}, and the free Fermion energies are represented by $\epsilon_{\bf k}$. The second line represents the Hamiltonian for the closed channel quasimolecular bosons $b^\dagger_{\bf q}$ ($b^{ }_{\bf q}$) with energies $E^0_{\bf q}+2\nu-\mu_B$. Finally, the last line represents the atom-molecule coupling discussed in section~\ref{sec:dyneqns}. The Fermi and Bose condensates exchange momentum ${\bf q}$ through this coupling. The Hamiltonian above can be written in terms of the generators of the $SU(2)\otimes SU(1,1)$ Lie algebra~\cite{huang:becbcs2}. Assuming that most of the Bosons are in the ground state given by ${\bf q}=0$, treating $b_0$ as a c-number, retaining only terms up to second order in $b^\dagger_{\bf q}$ ($b^{ }_{\bf q}$) and neglecting fluctuations about the mean field for all generators, the Hamiltonian can be diagonalized in a manner similar to that done for a pure Boson system in Bogoliubov theory~\cite{huang:becbcs2}, with the eigenstates given by the generalized $SU(2)\otimes SU(1,1)$ coherent state. If the system is in a pure quantum state, the energy eigenvalues are given by~\cite{huang:becbcs2}
\begin{equation}
E_{sn} = \left(2\nu-\mu_B \right)|b^2_0| + \sum_{\bf k} \left(2sE_{\bf k} +\epsilon_{\bf k}-\mu_F \right) +\sum_{\bf k\neq 0} \left[ \left( n+\frac{1}{2}\right)E^b_{\bf k} -\frac{1}{2}\left(E^b_{\bf k} + 2\nu-\mu_B+2u_B|b_0|^2\right)\right].
\end{equation}
Here, $E_{\bf k} = \sqrt{(\epsilon_{\bf k}-\mu_F)^2+|\Delta+g_rb_0|^2}$ are the energies of the Bogoliubov quasiparticles from conventional BCS theory, $E^b_{\bf k} = \sqrt{E^0_{\bf k}(E^0_{\bf k} + 2u_B |b_0|^2)}$ are the energies of the quasiparticle excitations of  the pure Boson system from Bogoliubov theory~\cite{huang:becbcs2}, and $s$/$n$ are indices that correspond to the ladder operators of the $SU(1,1)$/$SU(2)$ Lie algebras respectively, with $s=\pm\frac{1}{2}$ and $n\in\mathbb{N}^0$. If the energies are referenced from the ground state $E_{-\frac{1}{2}0}$, then the excitation energies are given by
\begin{eqnarray}
\label{eq:excitations}
 \delta E_{\frac{1}{2}n}&=& 2\sqrt{\mu^2_F + |\Delta + g_rb_0|^2} + \sum_{{\bf p}\neq 0}
\left(2E_{\bf p}+nE^b_{\bf p} \right),  \nonumber \\
\delta E_{-\frac{1}{2}n} &=& \sum_{{\bf p}\neq0} nE^b_{\bf p},
\end{eqnarray}
where $\epsilon_{{\bf k}=0}$ is taken to be $0$. Evaluating equation~\ref{eq:excitations} in $3$D with $\epsilon_{\bf k} = k^2/2m$, $E^0_{\bf p}=p^2/4m$ and a momentum cutoff at $p=p_c\equiv\frac{2\pi}{R}$ to counter ultraviolet divergences shows that the sums in these equations go as positive powers of the particle number density (noting that $\mathcal{V}=1$ here) $\rho\sim u_B |b_0|^2\sim\frac{1}{R^3}$, in particular,  $\delta E_{-\frac{1}{2}n}\sim \rho^{5/3}$. Thus, a lower bound for the excitation energies is given by
\begin{equation}
E_c = 2\sqrt{\mu^2_F + |\Delta + g_rb_0|^2}.
\end{equation}
$E_c\sim \rho^{1/2}$ (see equation~\ref{eq:numfb}), which is fairly macroscopic due to the presence of the effective gap $\Delta + g_rb_0$. The corresponding relaxation time $t_c\equiv 1/E_c$ defines the time scale below which excitations of energy greater than $E_c$ are imparted, and  excitations away from the ground state (such as closed channel bosons and unpaired Fermions dissociated from them) are possible. In all subsequent analyses, it is assumed that the quench time $t_q \leq t_c$ and both are neglected in comparison to the time scales of the post quench dynamics. This effectively makes the quench an instantaneous impulse, and an analysis of the dynamical system based only on a Boson field and Fermion pair field completely captures the non-equilibrium state.

If there is no BCS state in this system for finite values of $\epsilon_F$, then $\Psi_1 = \Psi_2 \equiv \Psi$, and the dynamics is equivalent to one with no atom-molecule coupling $g_r$. Thus, equations~\ref{dynamical_system} simplify to
\begin{eqnarray}
\label{dynamical_system_nogr}
\frac{\partial\Psi}{\partial t} - i\alpha \Psi +i\beta |\Psi|^2 \Psi &=& 0, \nonumber \\
\frac{\partial\Psi}{\partial t} +2 i \lambda_\nu \Psi +2 i\chi|\Psi|^2\Psi &=& 0.
\end{eqnarray}
Substituting quasiharmonic trial solutions $\Psi = \tilde{\Psi} e^{i\Omega t}$ yields
\begin{eqnarray}
\label{dynamical_system_nogr_ft}
\Omega \tilde{\Psi} - \alpha \tilde{\Psi} +\beta |\tilde{\Psi}|^2 \tilde{\Psi} &=& 0, \nonumber \\
\Omega\tilde{\Psi} +2  \lambda_\nu \tilde{\Psi} +2 \chi|\tilde{\Psi}|^2\tilde{\Psi} &=& 0.
\end{eqnarray}
Both of these equations need to hold for two unknowns $\tilde{\Psi}$, $\Omega$. Thus, in the absence of Boson interactions ($\chi = 0$), $\Omega=- 2 \lambda_\nu$, and the phase of the Boson condensate fraction oscillates without any population transfer out of the BEC in all regimes.

The shallow-BEC regime is one where the condensate fractions $n_F$ and $n_B$ are comparable (see figure~\ref{fig:chempot}). If no \emph{a-priori} assumptions are made about any of the other constants, equations~\ref{greek_consts_chempot} indicate that $\beta$ can be ignored. In the case of noninteracting Bosons, a quench from a pure BEC, where $\psi_{1,2}(0) = \left(1/\sqrt{2} \right)g_r/\epsilon_F$, to this regime (where the dynamics is linear) produces \begin{eqnarray}
\label{dynamical_system_nobeta}
\frac{\partial\Psi_1}{\partial t}+ i\gamma \left( \Psi_1-\Psi_2 \right) - i\alpha \Psi_1 &=& 0, \nonumber \\
\frac{\partial\Psi_2}{\partial t} +2i\lambda_\nu\Psi_2 - i\kappa\gamma \left(\Psi_1-\Psi_2 \right) &=& 0.
\end{eqnarray}
Substituting trial solutions $\Psi_{1,2}(t) = \tilde{\Psi}_{1,2} \times e^{i\Omega t}$, and solving the  resulting characteristic equation for $\Omega$,
\begin{equation}
\label{omegapm}
\Omega_\pm\left(\epsilon_F \right) = -\frac{1}{2} {\mathrm A}\left(\epsilon_F\right) \bigg[ 1 \pm \sqrt{1-4\frac{{\mathrm B}\left(\epsilon_F\right)}{{\mathrm A}^2\left(\epsilon_F\right)}}\bigg].
\end{equation}
where
\begin{eqnarray}
\label{AB}
{\mathrm A}\left(\epsilon_F\right) &\equiv& \gamma -\alpha +2\lambda -\kappa\gamma , \nonumber \\
{\mathrm B}\left(\epsilon_F\right) &\equiv& \left(\gamma-\alpha \right)\left(2\lambda+\kappa\gamma \right)-\kappa\gamma^2.
\end{eqnarray}
Here, as always the constants in Greek letters are functions of $\epsilon_F$. It is noted that, for sufficiently large $k_R$, $\mathrm{B}(\epsilon_F)$ is negative, and so $\Omega_\pm$ are both real (no damping or exponential growth). The BEC condensate fraction,  proportional to $|\Psi_2|^2(t)$,  oscillates with frequency $\Omega_{b_0} = |\Omega_+ - \Omega_-|$. Substituting the values above,
\begin{equation}
\Omega_{b_0} = {\mathrm A}\left(\epsilon_F\right)\sqrt{1-4\frac{{\mathrm B}\left(\epsilon_F\right)}{{\mathrm A}^2\left(\epsilon_F\right)}}.
\end{equation}
For sufficiently large cutoff frequency $k_R$, $\Omega_{b_0} \sim k_R$ and thus it is believed to represent Rabi oscillations of solitons that are similar to the ones reported in $2004$ by Andreev, Gurarie, and Radzihovsky~\cite{andreev:noneqmbcsbec}, as well as others since $2004$~\cite{barankov:bcsbecbloch,yuzbashyan2,miyakawa} for a quantum quenched noninteracting Fermi gas with a Feshbach resonance. Numerical results for $\Omega_{b_0}$ are plotted for $g_r=40$ in figure~\ref{fig:omega}.

A linear stability analysis of small displacements away from equilibrium should provide some qualitative insight into the long term dynamics of the system even for large displacements from equilibrium, such as impulse quenches. Absent the interference term $\Psi_1-\Psi_2$, eqns~\ref{dynamical_system} resemble the normal form of the Hopf bifurcation~\cite{strogatz:book}. However, this does not actually occur due to the complex nature of the coefficients. Writing the dynamical variables in polar coordinates as $\Psi_j = a_j e^{i\phi_j}$ and dropping $\Psi_1-\Psi_2$ simplifies the dynamics to yield $\dot{a}_j=0$ and $\phi_j\propto t$, trivial trajectories that do not admit to bifurcations of any kind.

Investigating the more general cases, however, do yield interesting results. Substituting polar coordinates into the dynamics in eqns~\ref{dynamical_system}, and splitting up the real and imaginary parts, results in a nonlinear dynamical system evolving in a $4$-dimensional phase space spanned by $a_{1,2}$ and $\phi_{1,2}$ viz.
\begin{eqnarray}
\label{dynamics:radial}
\dot{a}_1 &=& -\gamma a_2 \sin{\left(\phi_2-\phi_1\right)},\nonumber \\
\dot{a}_2 &=& \kappa \gamma a_2 \sin{\left(\phi_2-\phi_1\right)},\nonumber \\
\dot{\phi}_1 &=& \alpha - \beta a^2_1 - \gamma\left[1-\frac{a_2}{a_1}\cos{\left(\phi_2-\phi_1\right)} \right],\nonumber \\
\dot{\phi}_2&=& \kappa\gamma\left[\frac{a_1}{a_2}-\cos{\left(\phi_2-\phi_1\right)} \right]-2\lambda_\nu.
\end{eqnarray}
Here, only the first of two cases discussed in the previous section are considered, where the boson repulsion term $\chi$ has been ignored. The radial fixed points of $\dot{a}_j=0$ can be realized for nonzero $a_j$ if $\phi_2(t) = \phi_1(t)+n\pi$.  Note that, in this limit, $\cos{\left(\phi_2-\phi_1\right)}=\pm 1$ and the rhs of the last two of eqns~\ref{dynamics:radial} must be equal to each other. This is satisfied by eqn~\ref{chempot_final_noint}, once the fixed point values of $a_j$ from eqns~\ref{fixedpoints_noint_final} are taken as $a_1=r_\nu$, $a_2=(1-\eta_\nu)r_\nu$, with $r_\nu \equiv\left[\frac{1}{\beta}\left(\alpha -\gamma\eta_\nu \right)\right]^{1/2}$. Defining $\omega_\nu\equiv\dot{\phi}_1=\dot{\phi}_2$, and applying these expressions to the fixed points above, the system evolves at the radial fixed points as
\begin{eqnarray}
 \Psi_1(t) &=& r_\nu e^{-i\omega_\nu t},\nonumber \\
 \Psi_2(t) &=& (1-\eta_\nu)r_\nu e^{-i\omega_\nu t},\nonumber \\
\omega_\nu &=& \bigg\{\begin{array}{l}
                         \lambda_\nu\\
			 \lambda_\nu - 2\kappa\gamma
                       \end{array}.
\end{eqnarray}
Small displacements away from this orbit can be provided with complex quantities $\delta\Psi_j$ as follows.
\begin{eqnarray}
 \Psi_1(t) &=& \left[r_\nu + \delta\Psi_1(t)\right]e^{-i\omega_\nu t},\nonumber \\
 \Psi_2(t) &=& \left[\left(1-\eta_\nu\right)r_\nu + \delta\Psi_1(t)\right]e^{-i\omega_\nu t}.
\end{eqnarray}
Substituting this into the dynamics in eqns~\ref{dynamical_system} yields
\begin{eqnarray}
 \delta\dot{\Psi}_1 &=& -i\left(\omega_\nu+\gamma-\alpha+2\beta r^2_\nu\right)\delta\Psi_1-i\beta r^2_\nu \delta\Psi^\ast_1+i\gamma\delta\Psi_2 + C_1,\nonumber \\
\delta\dot{\Psi}_2 &=& i\kappa\gamma\delta\Psi_1 -i\left(\omega_\nu+2\lambda+\kappa\gamma\right)\delta\Psi_2 + C_2,
\end{eqnarray}
where $C_j$ are constants, and terms $\mathcal{O}(|\delta\Psi_j|^2)$ and higher have been dropped. Separating the system into real and imaginary parts yields $\delta \dot{x}_j \sim \sum_k \mathcal{J}^{jk}_\nu\delta x_k$, where $\delta\Psi_{1,2}=\delta x_{1,3} + i \delta x_{2,4}$, and ${j,k}=1-4$. The Jacobian matrix elements that describe this dynamics are given by $\mathcal{J}^{jk}_\nu\equiv\partial\dot{x}_j/\partial x_k$. The matrix evaluates to
\begin{equation}
\label{jacobian:nobosonint}
\mathcal{J}_\nu \approx \begin{pmatrix}
 0 & -[\gamma -\alpha]-3 \beta r^2_\nu  & 0 & \kappa \gamma  \\
 [\gamma -\alpha ]+\beta r^2_\nu  & 0 & \kappa \gamma  & 0 \\
 0 & \gamma  & 0 & -\kappa\gamma \left[\frac{\eta_\nu+1}{\eta_\nu-1}\right] \\
 -\gamma  & 0 & -\kappa\gamma\left[\frac{\eta_\nu-3}{\eta_\nu-1}\right] & 0
\end{pmatrix},
\end{equation}
at the fixed points in eqn~\ref{fixedpoints_noint_final}. The eigenvalues $\Omega_\nu$ are given by the characteristic equation $|\mathcal{J}_\nu-\Omega_\nu I|=0$. In this case, the equation is biquadratic and yields $\Omega_\nu^4 + \mathcal{B}_\nu\Omega_\nu^2-\mathcal{C}_\nu =0$, where
\begin{eqnarray}
\label{roots:jacobian:nobosonint}
\Omega_\nu &=& \pm\sqrt{-\frac{\mathcal{B}_\nu}{2}}\left(1\pm\sqrt{\mathcal{D}_\nu}\right)^{1/2},\nonumber \\
 \mathcal{D}_\nu&\equiv& 1+\frac{4\mathcal{C}_\nu}{\mathcal{B}_\nu^2},\\
\label{nobosonint:coeff}
\mathcal{B}_\nu &\equiv& \beta r_\nu^2 \left[2\left(\alpha-\gamma\right)+3\beta r_\nu^2 \right]-\left(\alpha-\gamma\right)^2-
\frac{\left(\eta_\nu-3\right)\left(\eta_\nu+1\right)}{\left(\eta_\nu-1\right)^2}\times\kappa^2\gamma^2,\nonumber \\
\mathcal{C}_\nu &\equiv& \kappa^2\gamma^4+\frac{\left(\eta_\nu-3\right)\left(\eta_\nu+1\right)}{\left(\eta_\nu-1\right)^2}
\left[\left(\alpha-\gamma-2r^2_\nu\beta\right)^2-r^4_\nu\beta^2 \right]\kappa^2\gamma^2+2\kappa^2\gamma^3\left[\alpha-\gamma-\left(\frac{2\eta_\nu-1}{\eta_\nu-1}\right)r^2_\nu\beta\right].\nonumber \\.
\end{eqnarray}
The quantity $\mathcal{D}_\nu$ is the \textit{discriminant} of the quartic equation. The linear dynamics is given by~\cite{strogatz:book} $|\delta x (t) \rangle \sim \sum^{4}_{j=1} c_j |\Omega^j_\nu\rangle e^{\Omega^j_\nu t}$, where $c_j$ are constants, $|\Omega^j_\nu\rangle$ are the eigenvectors of $\mathcal{J}_\nu$  that correspond to eigenvalues $\Omega^j_\nu$, and $|\delta x (t) \rangle$ is the vector given by $\delta x_{1-4}$. The sign of the quartic discriminant $\mathcal{D}_\nu$ plays a key role in determining the nature of these trajectories. The regions of interest and the nature of the roots $\Omega_\nu$ are tabulated in table~\ref{tab:roots}. If $\Omega_\nu$ are pure imaginary, the trajectories correspond to orbital motion about the radial fixed points. The orbits constitute two oscillations of frequencies $\Omega^{\pm}_\nu=|\sqrt{-{\mathcal{B}_\nu}/{2}}(1\pm\sqrt{\mathcal{D}_\nu})^{1/2}|$. Their combination results in periodicity iff $\Omega^\pm_\nu$ have rational winding numbers $n^+/n^- = \Omega^+_\nu/\Omega^-_\nu$, where $n^\pm$ are integers with no common factors. In that case, the period of the resultant oscillations are $T_\nu = \pi (n^+/\Omega^+_\nu + n^-/\Omega^-_\nu)$. The changes in the number density of the BEC, denoted by $|\delta \Psi_2|^2=\delta x^2_3 + \delta x^2_4$, will have $4$ oscillations viz $2\Omega^\pm_\nu$ and $\Omega^+_\nu\pm\Omega^-_\nu$. The net period is much smaller than that of the $\Psi_j$s. When the $\Omega_\nu$ start to pick up real parts, then the trajectories corresponding to negative real parts are stable and decay in spirals to the orbits of the radial fixed points, thus creating a stable limit cycle.
\begin{table}
 \begin{center}
\begin{tabular}{|l|l|l|l|}
\hline
		& $\mathcal{D}_\nu\leq0$    		             & $0<\mathcal{D}_\nu<1$		    &$\mathcal{D}_\nu											       		     \geq1$\\
\hline
		& $\Omega_\nu = \pm\sqrt{\frac{|\mathcal{B}_\nu|}{2}}
		  \left(1\pm i\sqrt{|\mathcal{D}_\nu|}\right)^{1/2}$ & $\Omega_\nu=\pm\sqrt{\frac{|\mathcal{B}_\nu|}{2}}
								\left(1\pm \sqrt{|\mathcal{D}_\nu|}\right)
												     ^{1/2}$
													    & $\Omega_\nu=\pm 														    \mbox{ }
													    \sqrt{\frac{|														    \mathcal{B}_\nu|}{2}}
												            \left(1+ 															    \sqrt{|			  												    \mathcal{D}_\nu|}														     \right)^{1/2}$\\
$\mathcal{B}_\nu<0$&						&					    & $\Omega_\nu=\pm i
													    \sqrt{\frac{|														    \mathcal{B}_\nu|}{2}}
												            \left(\sqrt{
													    |\mathcal 															    {D}|}-1\right)
													    ^{1/2}$\\
		& All roots are complex 			& All roots are real with		    & One pair of real
													      roots with \\			& in pairs of conjugates			& opposing signs			    &opposing signs
													     and one pair of\\
		& with opposing real parts			&		      	                    &imaginary roots
													     with opposing 												     		     signs\\
\hline
		&$\Omega_\nu = \pm i \sqrt{\frac{|\mathcal{B}_\nu|}{2}}
		  \left(1\pm i\sqrt{|\mathcal{D}_\nu|}\right)^{1/2}$&$\Omega_\nu=\pm i\sqrt{\frac{|\mathcal{B}_\nu|}{2}}
								\left(1\pm \sqrt{|\mathcal{D}_\nu|}\right)
												     ^{1/2}$&$\Omega_\nu=\pm i
													      \sqrt{\frac{|														      \mathcal{B}_\nu|}														       {2}}\left(1+ 															\sqrt{|			  												    	\mathcal{D}_\nu|}														      \right)
													      ^{1/2}$\\
$\mathcal{B}_\nu>0$ & 						&					    &$\Omega_\nu=\pm
													      \mbox{ }
													      \sqrt{\frac{|														      \mathcal{B}_\nu|}														       {2}}\left( 															\sqrt{|			  												    	\mathcal{D}_\nu|}														      -1 \right)
													      ^{1/2}$\\
		& All roots are complex 			& All roots are imaginary with		    & One pair of real
													      roots with \\			& in pairs of conjugates			& opposing signs			    &opposing signs
													     and one pair of\\
		& with opposing real parts			&		      	                    &imaginary roots
													     with opposing 												     		     signs\\
\hline \end{tabular}
\end{center}
\caption{Table of the eigenvalues $\Omega_\nu$ of the Jacobian $\mathcal{J}_\nu$. The table breaks the parameter space into $2$ regions (depicted as rows) demarcated by the sign of $\mathcal{B}_\nu$. In each zone, $\mathcal{D}_\nu$ further demarcates to $3$ sub-zones (depicted as columns). The eigenvalues, together with their natures, constitute the elements of the table.}
\label{tab:roots}
\end{table}
The eigenvalues with positive real parts depart further away from equilibrium in an unstable limit cycle. However, these trajectories are not expected to actually occur in these systems since they correspond to displacements that violate number conservation. In general, however, the radial fixed points have a saddle-like nature in this regime. In the case when $\Omega_\nu$ are purely real, then the radial fixed points remain saddles, but the decay to the radial orbit is non-spiral and very fast.

 Note that the frequency $\Omega_\nu$ depends on the renormalization cutoff $k_R$ through $\alpha$, and so depends on the range of the Feshbach resonance induced interaction. For small ranges, $\Omega_\nu\sim k_R$ and so this differs from the oscillations about equilibrium seen in BEC systems with no Feshbach resonance. In general, a profile of the discriminant $\mathcal{D}_\nu$ and the coefficient $\mathcal{B}_\nu$ as functions of the detuning $\nu$, obtained using eqns~\ref{nobosonint:coeff} and~\ref{roots:jacobian:nobosonint}, will give us the bifurcation diagram of this system for any choice of parameters. Figure~\ref{fig:bif:nobosonint} (a) and (b) show plots of $\mathcal{B}_\nu$ and $\mathcal{D}_\nu$ as functions of detuning for the choice of parameters shown in the caption. Here, $\mathcal{B}_\nu$ is positive for small detuning and drops below zero at a critical value as $|\nu|$ is increased. Before that happens, however, $\mathcal{D}_\nu$ starts from a positive value that is less than unity. As table~\ref{tab:roots} indicates, this means that the eigenvalues are pure imaginary. While $\mathcal{B}_\nu$ is still positive, the discriminant drops to negative values at $\nu\approx--88.6(0)$, and the roots  are now complex conjugate pairs. The value of $\mathcal{B}_\nu$ drops down to zero at around $\nu=-943.77(0)$, with the roots now vanishing at the origin. As $\mathcal{B}_\nu$ becomes negative, $\mathcal{D}_\nu$ is also negative and the roots retain their ongoing structure until $\mathcal{D}_\nu$ switches back to positive at around $\nu\approx-6715.(0)$, after which the roots are real with opposite signs as per table~\ref{tab:roots}. All these trends can be seen in the numerical evaluation of the eigenvalues plotted in the bifurcation diagram in figure~\ref{fig:bif:nobosonint} (c). Thus, a bifurcation is seen at the locus of $\mathcal{D}_\nu=0$ in the shallow-BEC regime. Here the limit cycles coalesce into orbital trajectories around the fixed points as the detuning $\nu$ is varied. However, this is not a Hopf bifurcation. The real parts of the eigenvalues do not cross over from the negative to the positive side as would be the case with a Hopf bifurcation~\cite{strogatz:book}.

The paragraph above demonstrates that the 'deep-BEC ' and 'shallow-BEC' regimes can be differentiated by the sign of the discriminant $\mathcal{D}_\nu$ inside these regimes. Finally, a rapid quench from the former to the latter can be expected to yield orbital phase space dynamics of $\Psi_j$ around the equilibrium fixed points after a sufficiently long time. Variations in the BEC number density, given by $|\Psi_2|^2 = x^2_3 + x^2_4$, is thus expected to show rapid nontrivial oscillations that can qualitatively resemble some sort of 'collapse and revival' of the matter wave packet.
\begin{figure}[h!bt]
\ \epsfig{file=fig7.eps,height=7.5in,width=4.0in}
\caption{(Color online) Bifurcation diagram in the complex plane of $\Omega_\nu$ obtained by numerically solving the BCS gap and number equations, applying the results to the Jacobian, and diagonalizing it. The numerics are performed for $100$ particles with $g_r=25$,$u_F = -0.3$,$u_B=0$, and renormalization cutoff $k_r=18.0$ in a unit volume for unit particle mass, where $\hbar$ has also been set to unity. Figure (a) shows the coefficient $\mathcal{B}$ defined in eq~\ref{nobosonint:coeff} as a function of detuning $\nu$. Note that $\mathcal{B} > 0$ until $\nu \approx -943.77(0)$ and then switches sign. Figure (b) shows the magnitude of the discriminant $\mathcal{D}$ as a function of detuning $\nu$ in a semi-log plot. The regions of positive (negative) $\mathcal{D}$ are indicated in blue (red). Figure (c) shows the actual bifurcation diagram, with the roots $\Omega$ shown in the complex plane. The origin, where $\mathcal{B}=0$ (see table~\ref{tab:roots}), is indicated, as well as the points where the discriminant $\mathcal{D}$ vanishes. The latter occurs at $\nu\approx-88.6(0)$ and $\nu\approx-6715.(0)$. The arrows indicate the direction where the magnitude of the detuning $\nu$ (given by $|\nu|$) increases in the range shown in figs (a) and (b).}
\label{fig:bif:nobosonint}
\end{figure}

Now, the dynamics is analyzed numerically for larger values of $\epsilon_F$ in the shallow-BEC regime. In this region, numerical methods are utilized to solve equations~\ref{dynamical_system}. The system is evolved numerically using Runge-Kutta-Fehlberg methods using the initial conditions $\Psi_{1,2} = g^2_r / 2 \epsilon_F$, corresponding to $|b_0|^2 = \mathcal{N}/2$ with $\mathcal{N}$ taken to be $100$. The quenching is assumed to have taken place by a sharp variation of $\epsilon_F$ (via the Feshbach detuning $\nu$) from a large negative value (in the deep-BEC regime) to the shallow-BEC region. Figures~\ref{fig:colrev}(a)-(h) show the time variation of the BEC condensate fraction for the two representative parameter sets from figure~\ref{fig:chempot} in section~\ref{sec:fpchempot}. Figures~\ref{fig:colrev}(a) and (e) show the time evolution for very small times ($ \sim 10^{-2}$ units). These signals are the master signals of the dynamics. The master signals appear to be Rabi oscillations about the equilibrium value $|\bar{\Psi}_2|^2$ with a frequency of $\Omega_m$. In figures~\ref{fig:colrev}(b) and (f), discrete time samples of the actual signals are plotted in steps of $2\pi/\Omega_m$. The presence of another oscillation(s) of frequency $\Omega_e$ indicates that the actual signals contain a superposition of several comparably fast Rabi oscillations interfering with each other, and that $\Omega_e$ is their \textit{beat frequency}. For larger times, note from figures~\ref{fig:colrev}(c) and (g) that these beats damp out over a long period of time $t_d$, until a sudden onset of partial \textit{collapse and revival} of the matter wave begins after a short initial relaxation time $t_r$ (indicated in the respective figures). The \textit{revival time} $t_R$ is also indicated in the respective figures. Approximate values of $t_d$, $t_r$, $t_R$ and $\Omega_{m,e}$ have been obtained from figures~\ref{fig:colrev} and can be seen in table~\ref{tab:colrev}.
The numerics plotted in fig~\ref{fig:colrev} indicate that the onset of this collapse and revival takes place at around $10^3$ numerical units. The values of $\epsilon_F$ , which are of the order of $10^2$ in numerical units, can be obtained from table~\ref{tab:colrev}. Note that these times are expressed in numerical units. In order to obtain them in actual units, it is necessary to start with the system trap size $L_u$, and note that the Fermion mass ($M$ in actual units, $m$ in numerical units) is taken to be unity in these numerical units. Thus, the relationship between numerical time $t$ and actual time $T$ can be obtained using eqns~\ref{greek_consts}, ~\ref{greek_consts_chempot}, the dimensionless scaling of time by $(\epsilon_F d)^{-1}$ in equations~\ref{dynamical_system} , and the fact that the base unit of time is $T_u = 2 M L^2_u/\hbar$. The relationship is
\begin{equation}
 L_u=2\left(\sqrt{\frac{2}{\epsilon_F}}\frac{\pi \hbar }{M}\frac{T}{t}\right)^{1/2},
\end{equation}
where the chemical potential $\epsilon_F$ is expressed in numerical units. Thus, the trap size $L_u$ can be adjusted for a particular $t$ until an experimentally accessible $T$ is attained. As an example, if the atoms involved are Fermionic isotopes of Lithium (with mass $M\sim 10^{-26}$ kg), and the onset of collapse and revival (at around $t\sim 10^3$) is desired to be of the order of milliseconds, then the equation above necessitates that the trap size be of the order of $300$ nanometers, with the larger confinements requiring longer times $T$ growing as  $\sqrt{T}$. Ultracold BCS systems with such dimensions have already been obtained using superconducting atom chips~\cite{atom:chips}, where trap lifetimes can range from seconds to minutes~\cite{atomchips:traps}. Therefore, such time scales are experimentally accessible. The similarities in order of magnitude with $\Omega_m$ from table~\ref{tab:colrev} and $\Omega_{b_0}$ from figure~\ref{fig:omega} suggest that the collapse and revival effect may be soliton-like Rabi oscillations reported in the literature. Also, note from these figures that all three time scales, $t_d$, $t_r$ and $t_R$, are larger for smaller $g_r$. Thus, it is expected that $t_d \rightarrow \infty$ as $g_r \rightarrow 0$, removing the collapse and revival effect. The master signal frequency is also expected to reduce to $2\lambda_\nu$ in this limit. It is also noted from figure~\ref{fig:chempot}(d) that the revival process characterized by $t_R$ seems to be interrupted by a signal that has the same temporal characteristics as the relaxation signal characterized by $t_r$, thereby inhibiting a total revival of the matter wave.

In the case of interacting Bosons, the parameters from figure~\ref{fig:chempot2} are used and  the coefficient $\beta$ neglected to simplify the dynamics. Thus, equations~\ref{dynamical_system} simplify to
\begin{eqnarray}
\label{dynamical_system_nobeta_int}
\dot \Psi_1 + i\gamma \left( \Psi_1-\Psi_2 \right) - i\alpha \Psi_1  &=& 0 \nonumber \\
\dot \Psi_2 +2 i \lambda \Psi_2 +2 i\chi|\Psi_2|^2\Psi_2 - i\kappa\gamma \left(\Psi_1-\Psi_2 \right) &=& 0,
\end{eqnarray}
This system is evolved numerically using the same initial conditions and numerical methods as detailed in the previous paragraph (however, the tolerance required for convergence was an order of magnitude lower). The Boson interaction amplitude was taken to be fairly small so as to see the effects of departing slightly from the $u_B=0$ case for equations~\ref{dynamical_system_nobeta}. For small values of $u_B$, only a damped oscillation was seen for all the time intervals that were allowed by the numerical tolerances. For $u_B = 0.6$, the collapse and revival effect was seen for large times, just as it was with the case of noninteracting bosons.

In both of these cases, the nonlinear contribution to the dynamics go as $|x|^2 x$, absent which the dynamics is oscillatory. Therefore, it is believed that this effect is caused by the nonlinear dynamics from this type of expression alone. The results of the numerical simulation for $u_B = 0.6$ are shown in figures~\ref{fig:colrev2}(a) through (d). The frequencies $\Omega_m$  and $\Omega_e$, as well as the time scales $t_d$, $t_r$ and $t_R$, are tabulated in table~\ref{tab:colrev}. Note that the decay time $t_d$ is much higher for the case of interacting bosons and no $\beta$, suggesting that the contribution of the Boson interaction to the collapse and revival is much smaller than that of the term that is nonlinear in the gap parameter and weighed by $\beta$. These striking phenomena should manifest in the time-of-flight absorption images of the gas.
\begin{table}
\begin{center}
\begin{tabular}{ | l  l  l  l|| l | l | l | l | l |}
      \hline
  $g_r$ 	& $-\nu$  & $\epsilon_F$ &$u_B$	& $\Omega_m$    & $\Omega_e$     & $t_d$ 	& $t_r$ & $t_R$ \\  \hline  \hline
  $25$		& $55$	  & $125.78$ 	 &$0$ 	& $1400$	& $10$		 & $400$ 	& $9$	& $6.0$   \\  \hline
  $40$		& $140$	  & $192.40$ 	 &$0$	& $3500$	& $6$		 & $125$	& $4$	& $2.5$ \\  \hline
  $25$		& $110$	  & $85.5$	 &$0.6$ & $1050$	& $1$		 & $600$	& $12$  & $8.0$ \\  \hline
     \end{tabular}
\caption{Table of master and envelope frequencies $\Omega_{m,e}$, as well as decay, relaxation and revival times ($t_{d,r,R}$ respectively) of the quenched dynamics shown in figures~\ref{fig:colrev} and~\ref{fig:colrev2}. The system is evolved for two values of $g_r$. The corresponding $\nu$s and $\epsilon_F$s are also shown.  {The values of the parameters are chosen as a representative sample, based on our choice of length and atom number scales as explained in section}~\ref{sec:fpchempot},  {and the conclusions are not affected qualitatively if the ranges are varied within these orders of magnitude and kept in the regions of interest.}}
\label{tab:colrev}
\end{center}
\end{table}

\section{Conclusion}
\label{sec:conclusion}
In this paper, the dynamics of an ultracold gas of Fermions with a narrow Feshbach resonance has been formulated in the dual-channel case after making the single mode approximation for Bosons. The mean field nonlinear and complex TDGL dynamics of the Fermions, coupled to the Gross Pitaevski dynamics for the composite Bosons, have been obtained. The Fermion dynamics is encapsulated in the temporal variation of the superfluid gap parameter $\Delta$, and that of the Bosons in the ground state amplitude $b_0$.

The equilibrium states for the system in the BEC side have been evaluated by looking at the adiabatic evolution of the stationary solutions of the dynamics, which leads to a population transfer from the BEC superfluid to the BCS superfluid for nondivergent Fermion interactions. The dynamics of the system as it is quenched from a pure BEC state to a state in the shallow-BEC  regime (accomplished by a rapid variation in the Feshbach detuning $\nu$) has been analyzed numerically. Nonlinearities in the dynamics cause the Rabi oscillations to relax back to it's equilibrium state. However, at large times, the relaxation gets interrupted and a
\textit{collapse and revival} type phenomenon of the Bosonic matter wave field ensues. This effect seems to be caused by the interference between the multiple modes of the dynamics, and shows a striking phenomenon that is analogous to the collapse and revival effect that has been reported experimentally for Bosons.
\section{Acknowledgements}
This paper was supported by a postdoctoral fellowship from the Department of Science and Technology, Government of India,  at the SN Bose National Centre for Basic Sciences. The author thanks both institutions for financing this work. The author also thanks Prof J.K. Bhattacharjee for his discussions and insights into the problem.
\begin{figure}[h!bt]
\hspace*{-0.3in}
\ \epsfig{file=fig5.eps,height=6.2in,width=6.2in}
\caption{(Color online) Plots of the signal $|\Psi_2|^2(t)$  from the solutions to equations~\ref{dynamical_system} for $u_B = 0$ after the system is quenched from a pure BEC to the shallow-BEC regime (see figure~\ref{fig:chempot}). Here, $\hbar = m =\mathcal{V} = 1$, and $u_F = -0.3$. Figures (a) - (d) show plots of the time evolution for $g_r=25.0$ and $\nu = -55.0$. Figures (e) - (h) show plots for $g_r=40.0$ and $\nu = -140.0$. Figures (a) and (e) show plots for very small times ($\sim 10^{-2}$). Figures (b),(c) and (f),(g) show plots of discrete time samples of the signal taken in units of the local time period from figures (a) and (e) respectively. Thus, figures (b),(c) and (f),(g) show the signals that envelope the master signal in (a) and (e) respectively. For figure (a), the master signal has frequency $\Omega_m \simeq 1400$, and for figure (e) $\Omega_m \simeq 3500$. Figures (b),(c) and (f),(g) show the same corresponding envelope plots for two different time intervals. Partial collapse and revival of the matter wave can be seen after a lengthy decay time, where the original envelope damps out as it oscillates. The envelope oscillation frequency $\Omega_e \simeq 10$ for figure (b),(c) and $\simeq 6$ for figures (f),(g). The decay times $t_d$ for figures (c) and (g) are approximately $400$ and $125$ respectively. The region where collapse and revival take place is magnified and shown in figures (d) and (h), where the entire signal (not just the envelope) is now being plotted. The presence of partial collapse and revival of the matter wave can be noted after a relaxation time $t_r$ as indicated in figures (d) and (h). The revival times $t_R$ are also indicated in these figures. In figure (d), $t_r\simeq9$, $t_R\simeq6$. In figure (h), $t_r\simeq4$, $t_R\simeq2.5$}
\label{fig:colrev}
\end{figure}
\pagebreak
\begin{figure}[h!bt]
\hspace*{-0.3in}
\ \epsfig{file=fig6.eps,height=4in,width=6in}
\caption{(Color online) Plots of the signal $|\Psi_2|^2(t)$  from the solutions to equations~\ref{dynamical_system_nobeta} for $u_B = 0.6$ after the system is quenched from a pure BEC to the shallow-BEC regime (see figure~\ref{fig:chempot}). Here, $\hbar = m =\mathcal{V} = 1$, $u_F = -0.3$, and $g_r = 25$.Figures (a) shows plots for very small times ($\sim 10^{-2}$). Figure (b) shows plots of discrete time samples of the signal taken in units of the local time period from figure (a) . Thus, figure (b) shows the signal that envelopes the master signal in (a). For figure (a), the master signal has frequency $\Omega_m \simeq 1050$. Figures (b) and (c) show the same corresponding envelope plots for two different time intervals. Partial collapse and revival of the matter wave can be seen after a lengthy decay time, where the original envelope damps out as it oscillates. The envelope oscillation frequency $\Omega_e \simeq 1$. The region where collapse and revival take place is magnified and shown in figure (d). The presence of partial collapse and revival of the matter wave can be noted after a relaxation time $t_r$ as indicated in figure (d). The revival times $t_R$ are also indicated in these figures. In figure (d), $t_r\simeq 12$, $t_R\simeq 8$.}
\label{fig:colrev2}
\end{figure}
\pagebreak


\section{First Light and the Rise of Primeval Galaxies}

State of the art numerical simulations reveal that the first stars in the universe
formed in isolation in small cosmological halos (spheroidal bound clumps of dark
matter and H/He gas) at redshifts $z \sim$ 20 - 30, or 200 - 400 Myr after the Big
Bang \citep{abn02,bcl02,nu01}. Soon thereafter, gravitational mergers congregated
these star forming halos into the first primitive galaxies. Preliminary results from
minihalo collapse models suggest that primordial (or Pop III) stars were very massive,
25 - 500 M$_{\odot}$ \citep{on07}.  They created H II regions that were enormous in
comparison to those in the Galaxy today, 2.5 - 5 kpc in radius \citep{wan04,ket04,abs06,
awb07}.  Primeval ionization fronts, or I-fronts, were key to the rise of the first
protogalaxies for several reasons. First, Pop III stars lacked strong winds, so their
H II regions determined their circumstellar environments.  They regulated the flow of
the first heavy elements from primordial supernovae (SNe) into the early IGM \citep{
get07,get10} and may even have resulted in the prompt formation of a second generation
of low-mass stars in SNe remnants \citep{wet08b}.  Second, the minihalos forming the
first stars themselves coalesced in small swarms due to cluster bias. Radiation fronts
from one star in the cluster could engulf nearby halos, photoevaporating them to varying
degrees and either suppressing or allowing new star formation in them \citep{sir04,i05,
su06,wet08a,su09,wet10}. Thus, cosmological I-fronts in part governed the number and
character of stars that were swept up into the first primitive galaxies.

Once primeval galaxies formed and sustained cyclical star formation, they too
propagated I-fronts into the IGM, beginning the process of cosmological reionization
by which the cold, neutral, mostly featureless universe of the recombination era was
gradually transformed into the hot, transparent cosmic web of galaxies and clusters
of galaxies we observe today \citep{rgs02,rgs08,wa08,wc09}.  Recent numerical models
have demonstrated that the I-fronts of both individual Pop III stars and protogalaxies
formed H$_2$ in their outer layers, a key coolant at high redshifts that may have
seeded the formation of additional structure in the early universe \citep{rgs01,wn08b}.
High-redshift I-fronts thus influenced the rise of structure on a variety of spatial
scales in the early cosmos.

In the Galaxy today, I-fronts from massive stars in the molecular cores of giant
molecular clouds (GMCs) are thought to trigger new star formation by perturbing
nearby cores on the verge of gravitational collapse \citep{da07,gr09}.  In this
manner, successive fronts can set off a chain of star formation throughout the
cloud.  On the other hand, dynamical instabilities in the fronts may degenerate
into turbulent flows that support gas against collapse, offsetting triggered star
formation. Yet, these same instabilities can form clumps that may later collapse
into stars.  How ionized flows regulate star formation in the Milky Way today
remains to be fully understood.

\section{Ionizing Radiation Hydrodynamics in ZEUS-MP}

We perform our ionization front calculations with ZEUS-MP \citep{wn06}, a massively
parallel Eulerian astrophysical hydrodynamics code that solves the equations of ideal
fluid dynamics\footnote{http://lca.ucsd.edu/portal/codes/zeusmp2}:
\vspace{0.1in}
\begin{eqnarray}
\frac{\partial \rho}{\partial t}  & = & - \nabla \: \cdotp \; (\rho {\bf v})  \\
\frac{\partial \rho v_{i}}{\partial t}  & = & - \nabla \: \cdotp \; (\rho v_{i}
{\bf v}) \: - \: \nabla p \: - \: \rho \nabla \Phi \: - \: \nabla \cdotp {\bf Q}    \\
\frac{\partial e}{\partial t}  & = & - \nabla \: \cdotp \; (e {\bf v}) \: - \: p\nabla \:
\cdotp \: {\bf v} \: - \: \bf{Q} : \nabla  {\bf v}. \\
&  &  \nonumber
\end{eqnarray}
Here, $\rho$, $e$, and the $v_{i}$ are the mass density, internal energy density,
and velocity at each mesh point and $p \;= (\gamma-1) e$ and {\bf{Q}} are the
gas pressure and the von Neumann-Richtmeyer artificial viscosity tensor. ZEUS-MP
evolves these equations with a second-order accurate monotonic advection scheme in
one, two, or three dimensions on Cartesian (XYZ), cylindrical (ZRP), or spherical
polar (RTP) coordinate meshes.  Our augmented version of the publicly-available
code self-consistently couples primordial gas chemistry \citep{wn08a} and
multifrequency photon-conserving UV radiative transfer \citep{wn08b} to fluid
dynamics for the radiation hydroynamical transport of cosmological I-fronts.

\subsection{Primordial H and He Chemistry}

We evolve H, H$^{+}$, He, He$^{+}$, He$^{2+}$, H$^{-}$, H$^{+}_{2}$, H$_{2}$, and
e$^{-}$ with nine additional continuity equations and the nonequilibrium rate
equations of \citet{anet97}: \vspace{0.05in}
\begin{equation}
\frac{\partial \rho_{i}}{\partial t} = - \nabla \: \cdotp \; (\rho_{i} {\bf v})
+ \sum_{j}\sum_{k} {\beta}_{jk}(T){\rho}_{j}{\rho}_{k} - \sum_{j} {\kappa}
_{j}{\rho}_{j}, \vspace{0.05in}
\end{equation}
where ${\beta}_{jk}$ is the rate coefficient of the reaction between species j and k
that creates (+) or destroys (-) species i, and the ${\kappa}_{j}$ are the radiative
rate coefficients.  We assume that the species share a common velocity distribution.
Mass and charge conservation, which are not guaranteed by either chemical or advective
updates, are enforced each time the fluid equations are solved.

Microphysical heating and cooling due to photoionization and gas chemistry is coupled
to the gas energy density by an isochoric update that is operator-split from updates
to the fluid equations:
\vspace{0.05in}
\begin{equation}
{\dot{e}}_{\mathrm{gas}} = \Gamma - \Lambda, \label{eqn: egas}
\vspace{0.05in}
\end{equation}
where $\Gamma$ is the cumulative heating rate due to photons of all frequencies and
$\Lambda$ is the sum of the cooling rates due to collisional ionization and excitation
of H and He, recombinations of H and He, inverse Compton scattering (IC) from the CMB,
bremsstrahlung emission, and H$_2$ cooling.

\subsection{Radiative Transfer}

Our photon-conserving UV transport \citep{mel06,anm99}, which is distinct from the
flux-limited diffusion native to the public release of ZEUS-MP, solves the static
approximation to the equation of transfer in flux form to compute radiative rate
coefficients for the reaction network at every point on the coordinate mesh. As
currently implemented, our code can transport photons from a point source
centered in a spherical grid or in plane waves along the $x$ or $z$-axes of Cartesian
or cylindrical meshes.  To model Pop III spectra, we discretize the blackbody photon
emission rates of a star with 40 uniform bins from 0.755 to 13.6 eV
and 80 logarithmically spaced bins from 13.6 eV to 90 eV, normalizing them by the
total ionizing photon rates by \citet{s02}.  The radiative reactions in our models
are listed in Table 1 of \citet{wn08b}.  We calculuate H$_2$ photodissociation
rates along rays parallel to the direction of radiation flow using self-shielding
functions modified for thermal broadening as prescribed by \citet{db96} to approximate
the effects of gas motion.  They are shown in equations 9 and 10 of \citet{wn08b}.

\subsection{Radiation Forces}

Radiation pressure due to ionizations occurs at the I-front itself and in recombining
gas in the H II region. \citet{wn06} found that the acceleration of fluid elements at
the front was large but momentary, and only alters the velocity of the front by 1 - 2
km s$^{-1}$.  Momentum deposition within the H II region is only prominent where gas
is very dense, like at the center of a minihalo being evaporated by its star.  There,
rapid successive cycles of ionization and recombination can impart
radiation forces to the gas that are hundreds of times the strength of gravity at early
times \citep[lower left panel of Figure 1 in][]{ket04}. As these forces propel gas near
the center of the halo out into the H II region, its densities and recombination rates
fall, so more ionizing photons from the star reach the I-front.  This higher flux
results in I-fronts that are faster than when such forces are not included.  In UV
breakout from the first star-forming halos, radiation forces speed up the D-type I-front
by 10 - 20\%. However, this effect is transient: after the internal rearrangement of gas
deep within the H II region dilutes its interior, radiation forces there sharply fall
 \citep[see also][]{krum09,draine10,lopez10}.  We describe the timescales on which
 momentum updates are performed below.

\subsection{Adaptive Subcycling}

A hierarchy of highly disparate characteristic time scales arises when gas dynamics,
radiative transfer, and primordial chemistry are solved in a given application.  The
three governing times are the Courant time, the chemical time
\begin{equation}
t_{chem} = 0.1 \, \displaystyle\frac{n_{e} + 0.001 n_{H}}{{\dot{n}}_{e}},
\end{equation}
and the photoheating/cooling time \vspace{0.05in}
\begin{equation}
t_{hc} = 0.1 \displaystyle\frac{e_{gas}}{{\dot{e}}_{ht/cool}}.\vspace{0.05in}
\end{equation}
Their relative magnitudes can seamlessly evolve during a single application.
For example, when an I-front propagates through a medium, photoheating times are
often smaller than Courant times, and chemical time scales are usually shorter than
either one. On the other hand, fossil H II regions can cool faster than they
recombine, so cooling times become shorter than chemical times. The key to solving
all three processes self-consistently is to evolve each on its own timescale without
restricting the entire algorithm to the shortest of the times. To successfully deal
with both I-fronts and relic H II regions, an algorithm must adaptively reshuffle
the time scales on which the three processes are solved. Implicit schemes are
sometimes applied to stiff sets of differential equations like those in our
model because they are unconditionally stable over the Courant time.  However,
accurate I-front transport in stratified media often requires restricting updates
to both the gas energy and fluid equations to photoheating times in order to capture
the correct energy deposition into the gas, and linear system solves over such short
time scales would be prohibitive in more than one or two dimensions.  Enforcing photon
conservation in implicit schemes can also be problematic.

We instead subcycle chemical species and gas energy updates over the minimum of the
chemical and heating/cooling times until the larger of the two has been crossed, at
which point we perform full hydrodynamical updates of gas densities, energies, and
velocities.  These times are global minima for the entire grid.  The chemical times
are defined in terms of electron flow to accommodate all chemical processes rather
than just ionizations or recombinations.  Adopting the minimum of the two times for
chemistry and gas energy updates enforces accuracy in the reaction network when $
t_{chem}$ becomes greater than $t_{hc}$, such as in relic H II regions.  Our
adaptive subcycling scheme is described in greater detail in \citep{wn08a}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=19pc]{IF_structure.eps}
\end{center}
\caption{\label{fig1}
The structure of an ionization front. $\chi_e$ denotes ionized fraction and T is gas
temperature.  Top: monoenergetic UV photons just above the ionization edge of H.
Center: a 10$^5$ K Pop III star blackbody spectrum.  Note the significant  broadening
of the front.  Bottom: power-law quasar x-ray source spectrum.  Note that in the first two
cases that temperature tracks the ionization fraction, but not in the quasar case. Here,
hard photons penetrate well ahead of the front, strongly heating but not significantly
ionizing the gas.
}
\end{figure}

\section{Ionization Front Instabilities}

When a point source of monoenergetic ionizing photons appears in a uniform H and
He medium, they are not free to simply stream into the interstellar medium (ISM)
because of the large cross section to photoionization in the gas. Instead, they emerge
as a sharp wall of radiation, behind which gas is almost fully ionized and heated
to a few$\times 10^4$ K and beyond which it is completely neutral. The thickness
of the I-front separating the neutral and ionized gas is a few mean free paths of
the ionizing photons through the neutral gas.  The initial propagation of the
I-front through the gas is supersonic in that the ionization wave advances far
more quickly than any hydrodynamic response by the hot gas in its wake. This is known
as an R-type I-front.  As the front recedes from the star it slows, partly due to
geometric dilution of the ionizing flux and partly because recombinations in the
ionized gas remove photons that would otherwise have gone on to advance the front.
If the medium is static (no gas motion allowed), the I-front comes to a halt when
the H II region it bounds encloses enough ionized gas that the global recombination
rate equals the emission rate of source photons.  The radius at which this occurs
is known as the Str\"omgren radius.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=19pc]{H2_instability.eps}
\end{center}
\caption{\label{fig2}
Thin-shell instability in a primordial I-front mediated by H$_2$ cooling.
}
\end{figure}

In reality, as the I-front slows to the Str\"omgren radius the pressure wave that
has built up in the hot ionized gas behind it overtakes it and passes through it,
steepening into a shock as it does so.  Thereafter, the H II region expands
subsonically with respect to the sound speed in its interior, which is still
supersonic with respect to the ambient neutral gas.  As it expands, it bypasses
the static Str\"omgren radius limit because its interior densities and therefore
recombination rates fall, allowing source photons to continue to stream to the front.
As the H II region grows it pushes aside more and more ambient gas which accumulates
on its surface, causing the shock to detach from and move ahead of the I-front.  This
is known as a D-type I-front.

\subsection{Thin-Shell Instabilities}

If this shocked neutral shell can radiatively cool, it will collapse into a
cold dense layer that is prone to fragmentation via Vishniac, or thin-shell,
instabilities \citep{v83}.  Unlike similar features in the dense shells of
astrophysical blast waves, these erupt into very violent instabilities because
ionizing radiation opportunistically escapes through the cracks in the shell
in jets ahead of the front \citep{gsf96,wn08a}.  These jets rapidly crumple into
complicated flows that form dense clumps that can persist for the life of the
star. In the Galaxy today, metals and dust radiatively cool the shell and trigger
such instabilities.  Cooling by atomic H and He lines in the shocked gas of
cosmological I-fronts by itself is too weak to incite such instabilities. However,
they can instead form by H$_2$ cooling.

\begin{figure*}
\epsscale{2.0}
\plotone{shadow.eps}\vspace{0.15in}
\caption{The shadow instability.}
\vspace{-0.125in}
\end{figure*}

Primordial I-fronts are driven by hard UV sources whose photons have a variety of
mean free paths through neutral gas.  Collectively, they broaden the front, which
leads to ionized fractions of 10\% and temperatures of 2000 - 3000 K in its outer
layers as we illustrate in Figure 1.  These are ideal conditions for gas phase
catalysis of H$_2$ by the H$^-$ and H$_2^+$ channels:
\begin{equation}
H + e^- \rightarrow H^- + \gamma \hspace{0.25in} H^- + H \rightarrow H_2 + e^-
\end{equation}
\begin{equation}
H + H^+ \rightarrow H_2^+  + \gamma \hspace{0.25in} H_2^+ + H \rightarrow H_2 + H^+ \vspace{0.1in}
\end{equation}
Although much less efficient than metals, H$_2$ is a far more effective coolant
than H and He lines at 3000 K; sandwiched between the front and the shocked neutral
shell, it cools the base of the shell sufficiently to trigger thin-shell
overstabilities that erupt into violent radiation-driven instabilities, as we show
in Figure 2 \citep{wn08b}.  In this simulation, we center a 10$^5$ K blackbody UV
source approximating a 120 M$_{\odot}$ Pop III star in a H/He gas profile similar
to that of a cosmological minihalo. When the front becomes D-type, instabilities
due to H$_2$ cooling can be seen to appear at the base of the dense shell.  More
models are required to determine if the turbulent gas clumps survive until the
death of the star and are prone to gravitational collapse.

\subsection{Shadow Instabilities}

In Figure 3 we show the breakout of a second kind of radiation hydrodynamical
instability in primordial I-fronts, the shadow instability \citep{rjw99,wn08b}.
These form when density perturbations are advected through an R-type front. The
perturbation does not trap the front, it merely dimples it.  The walls of the
dimple are exposed to a smaller ionizing flux because they are oblique to it; as
a consequence, the walls of the dimple transition to D-type before the tip does,
causing the tip to advance and elongate the dimple into a jet that abruptly breaks
down into a violent instability.  Such phenomena would have been unavoidable in the
primeval universe, on scales ranging from UV breakout from the first star-forming
halos to 100 kpc protogalactic I-fronts.  The role of H$_2$ formation, thin-shell
instabilities, and shadow instabilities in structure formation in early cosmological
flows remains to be properly investigated.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=19pc]{outside_in.eps}
\end{center}
\caption{\label{fig4}
Angle of incidence (AOI) instabilities in the ionization front of a Pop III star
enveloping a neighbor halo.
}
\end{figure}

\vspace{-0.05in}

\subsection{AOI Instabilities}

We show in Figure 4 the radiation front of a 25 M$_{\odot}$ Pop III star enveloping
a nearby star-forming halo \citep{wet10}. Spikes in the front can be seen all along
the arc surrounding the halo, but are increasingly prominent with altitude above its
axis. This is an example of an I-front instability described by \citet{rjw02}, in
which photons that arrive at the front at any angle other than 90$^{\circ}$ cause it
to be unconditionally unstable: an angle of incidence (AOI) instability.  Detailed linear
perturbation analysis predicts that the greater the angle of deviation of the
photons from pure incidence, the faster the growth of the unstable mode, which true
of the protrusions in Figure 4.  In this example, the unstable modes have little
effect on star formation in the halo because those closest to its core have the
smallest amplitudes and fail to puncture it. Dynamical instabilities have also been
found in I-fronts engulfing Galactic molecular cloud cores in numerical models, but
gas in those simulations could efficiently cool by fine-structure lines \citep{miz05,
miz06} so the unstable modes are likely thin-shell in origin.

\bibliographystyle{spr-mp-nameyear-cnd}