\section{}

\section{Introduction}
\label{SEC-INTRO}

In 1957, Novick and Weiner discovered that {\em Escherichia coli} can
exhibit discontinuous switching in expression of the {\em lac} operon,
with some cells expressing a large amount of $\beta$-galactosidase
($\beta$-gal), other cells expressing a small amount, and an
insignificant number of cells expressing an intermediate
amount~\cite{Novick57}. Recently, this effect was further
characterized using single-cell assays of fluorescence levels in a
population of {\em E. coli} cells carrying a {\em lac::gfp}
reporter~\cite{Ozbudak04}. The population exhibited a bimodal
distribution, with induced cells having over 100 times the
fluorescence level of uninduced cells. These observations have been
attributed to the existence of two steady states, i.e., bistability,
in the induction of {\em lac} in {\em E. coli}.

Recent modeling studies have emphasized the importance of determining
whether bistability in expression of {\em lac} is relevant within a
natural context
\cite{Savageau01,Savageau02,Yildrim03,Yildrim04,vanHoek06}. This
question remains open because experimental studies have focused on the
response of {\em lac} expression to artificial inducers, such as
thio-methylgalactoside (TMG) and isopropyl-$\beta$,
D-thiogalactopyranoside (IPTG), rather than the natural inducer,
allolactose. This difference is critical because artificial inducers
(also known as gratuitous inducers) are not metabolized by the induced
enzyme, whereas the natural inducer is a metabolic intermediate in
lactose degradation, which is catalyzed by the induced enzyme.

Savageau~\cite{Savageau01} found important differences between
induction by IPTG vs. lactose in his theoretical treatment of
bistability in the {\em lac} operon. In Savageau's model, because
production and decay of allolactose are both proportional to the
$\beta$-gal concentration, bistability is forbidden. Expression of
{\em lac} in response to lactose was therefore predicted not to
exhibit bistability. This prediction agreed with the absence of {\em
steady-state} bistability in an experimental study of populations of
{\em E. coli} cells exposed to lactose, described in the Supplementary
Material of Ref.~\cite{Ozbudak04}---in that study, only {\em
transient} bimodal distributions of green fluorescence levels among
cells were observed at some glucose concentrations. It was later noted
that models with operon-independent decay of lactose (e.g., due to
dilution by cell growth) could exhibit
bistability~\cite{vanHoek06}. Several studies using such models found
either a bistable or graded response to lactose, depending on
parameter values or external glucose
levels~\cite{Yildrim03,Yildrim04,vanHoek06,vanHoek07,Santillan07}, and
in agreement with the model of Savageau, a model of van Hoek \&
Hogeweg~\cite{vanHoek06} was explicitly shown to exhibit no
bistability in the absence of operon-independent decay of
allolactose. However, these studies disagree in their assessment of
whether bistability is present~\cite{Yildrim03,Yildrim04,Santillan07}
or absent~\cite{vanHoek06,vanHoek07} in expression of {\em lac} among
{\em E. coli} cells in a natural context.

In addition to predicting whether {\em lac} induction exhibits
bistability, some studies have addressed the question of whether
bistability might enhance or hinder the performance of {\em E. coli}
cells. Both Savageau~\cite{Savageau02} and van Hoek \&
Hogeweg~\cite{vanHoek07} found that bistability increases the time
required to respond to sudden increases in environmental lactose,
which can be a disadvantage in competition for nutrients. These
results argue against the natural relevance of bistability in {\em
lac} expression.

Another important question that has not yet been addressed is whether
the experimental observations of bistability in Ref.~\cite{Ozbudak04}
are consistent with independent biophysical data that characterize
processes relevant to regulation of {\em lac} expression. Although
phenomenological models were developed to reproduce the steady-state
behavior~\cite{Ozbudak04} and the experimentally characterized
dynamics of switching between stable steady states~\cite{Mettetal06},
these models were not constrained by independent biophysical data. For
example, it is unclear whether the phenomenological models are
consistent with independently measured permease transport kinetics. On
the other hand, studies of bistability using more detailed,
biophysical models of {\em lac} induction were either only partially
constrained~\cite{vanHoek06} or did not consider the response to
artificial inducers~\cite{Yildrim03,Yildrim04,Santillan07}.

Here we analyze bistability in an ordinary differential equation (ODE)
model of {\em lac} induction. We use ODEs because we restrict our
analysis to steady-state behaviors, and because the protein
concentrations in fully induced cells are O($10^4$) per cell (see
Parameter Values section) and have negligible fluctuations. Similar
equations describe induction by lactose or artificial inducers. We
first use the model to gain insight into key determinants of
bistability of {\em lac} expression in response to artificial
inducers, and to understand how characteristics of bistability are
controlled by model parameters. We then use the resulting insight to
tune the parameters of the model to match the bistable behavior
observed by Ozbudak et al.~\cite{Ozbudak04}, and to predict mechanisms
by which bistability might be abolished. Finally, like previous
modeling studies, we use the model to address the question of whether
{\em lac} expression might be bistable in a natural context,
contributing to resolution of what is now a long-standing controversy.

\section{Model}
\label{SEC-MOD}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=5.0in]{Fig_1}
\end{center}
\caption{Circuitry for models of {\em lac} induction.  a) Model for
induction by lactose (Eqs.~(\ref{allo})), including the following
processes: (1) proportional production of permease (LacY) and
$\beta$-gal (LacZ); (2) permease-mediated transport of lactose; (3)
dilution of intracellular species by cell growth; (4) $\beta$-gal
catalyzed degradation of lactose, producing both the metabolic
intermediate allolactose, and the ultimate products of degradation,
glucose and galactose; and (5) $\beta$-gal catalyzed degradation of
allolactose, producing glucose and galactose. b) Model for induction
by artificial inducers (Eqs.~(\ref{iptg})), including: (1)
proportional production of permease (LacY) and $\beta$-gal (LacZ); (2)
permease-mediated transport of inducer; (3) dilution of intracellular
species by cell growth and (6) passive transport of inducer.}
\label{myfig}
\end{figure}

In our model of {\em lac} induction (Fig.~\ref{myfig}a), the following
set of coupled ordinary differential equations relate the internal
lactose concentration ($l$), allactose concentration ($a$), and
$\beta$-galactosidase concentration ($z$) to the external lactose
concentration ($l^{*}$)
\begin{subequations}
\label{allo}
\begin{eqnarray}
\dot{l} & = &  \alpha z \frac{\left(l^{*} - \phi l\right)}{K_{i} + l+ l^{*}} - \frac{\beta z l}{\left(1+a/K_{m,a}\right)K_{m,l} + l} - \gamma l, \label{allo1}\\
\dot{a} & = & \frac{\nu\beta z l}{\left(1+a/K_{m,a}\right)K_{m,l} + l} - \frac{\delta z a}{\left(1+l/K_{m,l}\right)K_{m,a} + a} - \gamma a\ \ \mathrm{and}	 \label{allo2}\\
\dot{z} & = & c\gamma +  \frac{\epsilon\gamma a^{n}}{K_{z}^{n} + a^{n}} - \gamma z. \label{allo3}
\end{eqnarray}
\end{subequations}
In Eqs.~(\ref{allo}), $\alpha$ and $\phi\alpha$ are the rate constants for
permease-dependent lactose import and export, respectively, $K_i$ is
the Michaelis constant for permease-dependent lactose transport
(assumed to be the same for import and export), $\beta$ and $K_{m,l}$
are the rate constant and Michaelis constant for lactose degradation,
$\nu$ is the branching fraction of lactose degradation to allolactose,
$\delta$ and $K_{m,a}$ are the rate constant and Michaelis constant
for allolactose degradation, $\gamma$ is the rate of dilution due to
cell growth, $c\gamma$ and $\epsilon\gamma$ are the basal and
inducible rates of $\beta$-galactosidase production, $K_z$ is the
allactose concentration at half-maximal induction of
$\beta$-galactosidase production, and $n$ is the Hill number for
lactose induction of $\beta$-galactosidase production.

The metabolic fluxes in Eqs.~(\ref{allo}) include the effects of
competition between allolactose and lactose for access to
$\beta$-galactosidase ($\beta$ and $\delta$ terms). Because shuttling
of galactosides across membranes occurs through a single permease
channel \cite{Lolkema91}, we also consider the influence of
competition between external and internal lactose for access to
permease ($\alpha$ and $\phi$ terms); however, as a simplification, we
do not consider transitions among distinct internal states of the
permease \cite{Lolkema91}.

To focus on the operating conditions of
the system that are most relevant to lactose utilization by {\em
E. coli}, we only consider regulation in the absence of glucose. This
focus is appropriate because, in the presence of glucose, {\em lac} is
not essential for growth, and induced $\beta$-galactosidase levels are
low \cite{Magasanik87}.

Similarly, the model of artificial induction of {\em lac}
(Fig.~\ref{myfig}b) is given by
\begin{subequations}
\label{iptg}
\begin{eqnarray}
\dot{l} & = & \alpha_{0} \left(l^{*} - l \right) + \alpha z\frac{ \left(l^{*} - \phi  l\right)}{K_{i} + l+ l^{*}}  - \gamma l\ \ \mathrm{and}		\label{iptg1} \\
\dot{z} & = & c\gamma + \frac{\epsilon\gamma l^{n}}{K_{z}^{n} + l^{n}} - \gamma z.	\label{iptg2}
\end{eqnarray}
\end{subequations}
In Eqs.~(\ref{iptg}), variables and parameters have the same meaning as
in Eqs.~(\ref{allo}), except $l$ and $l^{*}$ correspond to the level of
internal and external artificial inducer (e.g., IPTG or TMG),
respectively, and $\alpha_0$ is the rate constant for leakage across
the membrane.

In Eqs.~(\ref{allo}) and Eqs.~(\ref{iptg}), protein expression is
lumped with gene expression, and the dependence of promoter activity
on the level of signal (IPTG, TMG, or allolactose) is modeled using a
simple Hill function, which is significantly simpler than other models
\cite{Wong97,Yildrim03,Yildrim04,vanHoek06,Narang07,vanHoek07,Santillan07}. On
the other hand, Eqs.~(\ref{allo}) considers the effects of competition
among substrates in permease transport and metabolic processes, unlike
other models of {\em lac} induction
\cite{Savageau01,Savageau02,Vilar03,Ozbudak04,Wong97,Yildrim03,Yildrim04,vanHoek06,vanHoek07,Santillan07}. Compared
to the model of Savageau \cite{Savageau01,Savageau02},
Eqs.~(\ref{allo}) considers operon-independent decay of allolactose,
without which bistability in response to lactose is impossible
\cite{Savageau01,Savageau02,vanHoek06}, as discussed above. Overall,
Eqs.~(\ref{allo}) and Eqs.~(\ref{iptg}) are less detailed than the
{\em lac} induction models used in Refs. \cite{Wong97},
\cite{Yildrim03}, \cite{Yildrim04}, \cite{vanHoek06},
\cite{vanHoek07}, and \cite{Santillan07}, and are more detailed than
those used in Refs. \cite{Savageau01}, \cite{Savageau02},
\cite{Vilar03}, and \cite{Ozbudak04}, and they therefore constitute
intermediate complexity equations describing {\em lac}
induction. Compared to the simpler models, the intermediate level of
detail provides increased contact between model parameters and
biophysical measurements, and compared to more detailed models, it
facilitates analysis of the equations and interpretation of the
results.

\section{Parameter values}

We used the parameter values and ranges listed in
Table~\ref{Tableparams} to analyze bistability in Eqs.~(\ref{allo})
and Eqs.~(\ref{iptg}). The values in the table were obtained as follows:

\begin{itemize}
\item $\gamma$. We assume the generation time under the conditions in
Ref.~\cite{Ozbudak04} is 30-60 min. We note, however, that this time
might be very different for {\em E. coli} growing under stress in the
gut; this represents a source of uncertainty concerning the biological
relevance of our predictions.

\item $\alpha_0$. We assume that $\alpha_0=0$ except for the case of
IPTG, where we explore a range consistent with that considered in
Ref.~\cite{vanHoek07}.

\item $\alpha$. An approximate range of 1-100~s$^{-1}$ for sugar
transport turnover numbers was obtained from the review by Wright et
al.~\cite{Wright86}. The range is broader than measured values
\cite{Viitanen84} because measurements were made at 25$^\circ$C rather
than at the physiological temperature of 37$^\circ$C in the host
environment of the gut that we are focusing on here, and at which
measurements in Ref. \cite{Ozbudak04} were performed. The nominal
value of $1000\ \mathrm{min}^{-1}$ was estimated from
Ref.~\cite{Viitanen84} assuming the production rate of permease is the
same as that of functional $\beta$-gal. Because permease is a monomer
while $\beta$-gal is a tetramer, this assumption entails a four-fold
smaller production rate for permease. This seems possible, as (1)
galactoside acetyltransferase (GATase) monomer synthesis is eight-fold
smaller than $\beta$-gal monomer synthesis; (2) due to incomplete
operon transcription and the order of genes in the operon ({\em
lacZYA}), the amount of mRNA transcribed from the GATase gene ({\em
lacA}) and permease gene ({\em lacY}) is smaller than that from the
$\beta$-gal gene ({\em lacZ}); (3) there is some evidence that
permease is made in smaller amounts than $\beta$-gal \cite{Kennedy70}.

\item $\phi$. We assume no export flux through permease in the
artificial induction model, and then examine the consequences of
introducing such a flux on bistability. Guided by
Ref.~\cite{Lolkema91}, for the lactose model, we use a nominal efflux
rate constant ($\phi\alpha$) of half the value of the influx rate
constant $\alpha$, and allow the value to decrease in the search for
bistable conditions.

\item $K_i$. For simplicity, we assume the same Michaelis constant for
permease import and export--a nominal value of $0.5$~mM was obtained
from Ref.~\cite{Viitanen84}. The range was applied as per $\alpha$,
and encompasses measured values \cite{Wright81,Page84,Viitanen84}.

\item $\beta$. A total lactose turnover number for
$\beta$-galactosidase of $2.85\times 10^4$~min$^{-1}$ is estimated
from a measured value of $V_{max}=61.3~\mu$mol~min$^{-1}$~mg$^{-1}$ in
Ref.~\cite{Huber76}. This estimate is an order of magnitude greater
than the value $3.6\times 10^3$~min$^{-1}$ given in
Ref.~\cite{Martinez91}, but the two estimates agree closely when one
considers that $\beta$-gal converts about half of its lactose
substrate to glucose and galactose, rather than allolactose, and that
the enzyme is composed of four monomeric catalytic subunits. The
estimate given in Ref.~\cite{Martinez91} is appropriate for total
turnover of lactose on a per monomer basis. Like for $\alpha$, because
measurements were performed at 30$^\circ$C, we consider a range of
values ten times lower to ten times higher than the nominal value.

\item $K_{m,l}$.  The nominal value was obtained directly from
Ref.~\cite{Huber75}. As for $\beta$, because of temperature
considerations, we use a range from ten times lower to ten times
higher than the nominal value.

\item $\nu$. The value $\nu=0.468$ was calculated from the total rate
of $\beta$-gal degradation of lactose and the partial flux from
lactose to allolactose reported in Ref.~\cite{Huber76}. We take it to
be a constant because the ratio of reaction products was found to be
insensitive to temperature changes between $30^\circ$C and $0^\circ$C.

\item $\delta$. An allolactose turnover number for $\beta$-gal of
$2.3\times 10^4$~min$^{-1}$ is estimated from a measured value of
$V_{max}= 49.6~\mu$~mol~min$^{-1}$~mg$^{-1}$ in
Ref.~\cite{Huber75}. As for $\beta$, because of temperature
considerations, we use a range from ten times lower to ten times
higher than the nominal value.

\item $K_{m,a}$.  The nominal value was obtained directly from
Ref.~\cite{Huber75}. As for $\beta$, because of temperature
considerations, we use a range from ten times lower to ten times
higher than the nominal value.

\item $\epsilon$. Using a production rate of 5 $\beta$-gal tetramers
per cell per second for a 48 min generation time \cite{Kennell77},
14,400 molecules are produced during a generation at full
induction--this is the number of molecules in the cell after doubling
(supporting our choice of a noiseless model). Assuming a 1
$\mu\mathrm{m}^3$ mean cell volume \cite{Kubitschek86} and linear
volume increase in time \cite{Kubitschek90}, the volume after doubling
is approximately 0.7 $\mu\mathrm{m}^3$, leading to a concentration of
34,286 $\mathrm{nM}$.

\item $c$. This value is derived from $\epsilon$, assuming a 1000-fold
increase in $\beta$-galactosidase levels upon induction
\cite{Beckwith87}.

\item $K_z$ and $n$. These values are estimated from IPTG induction
data in permease knockout cells both from Ref.~\cite{Sadler65},
Fig. 15 and from data compiled in Ref.~\cite{Yagil71}, Figs. 1 and
2. The nominal value $n=2$ was estimated from the slopes of the curves
in the figures, and $K_z$ was determined by estimating from the
figures the concentration of IPTG at half-maximal induction. The
nominal value of $10^5\ \mathrm{nM}$ was estimated from data compiled
in Ref.~\cite{Yagil71}. To determine the range, an approximate lower
value of $10^4\ \mathrm{nM}$ was obtained from Ref.~\cite{Sadler65},
and we allowed for an upper value of $10^6\ \mathrm{nM}$ to account
for potential differences between induction by IPTG and TMG or
lactose.

\end{itemize}

\begin{singlespace}
\begin{table}[htd]
\caption{Parameter values.}
\begin{center}
\begin{tabular}{|c|l|c|c|}
\hline
Param & Description & Nominal & Range \\ \hline
	$\gamma$	& growth rate	& -- &	 $0.0116\ \mathrm{min}^{-1}$ -- $0.0231\ \mathrm{min}^{-1}$						\\ \hline
	$\alpha_{0}$ 	& passive transport	&0 &	$0$ -- $1.35\ \mathrm{min}^{-1}$								\\
&rate constant& &\\ \hline
	$\alpha$ 	& permease import	&$600\ \mathrm{min}^{-1}$ & 	$6\times 10^{1}\ \mathrm{min}^{-1}$ -- \\
& turnover number & &$6\times 10^{3}\ \mathrm{min}^{-1}$\\ \hline
	$\phi$		& ratio of permease	&0 (artificial inducers)&	 0 -- $0.5$							\\
&export to import& or 0.5 (lactose)&\\
&turnover numbers& &\\ \hline
	$K_{i}$		& permease Michaelis	& $5\times 10^5\ \mathrm{nM}$&	 $5 \times 10^{4}\ \mathrm{nM}$ --\\
&constant& &$5\times 10^6\ \mathrm{nM}$ \\ \hline
	$\beta$		& $\beta$-gal lactose	&$2.85\times 10^4\mathrm{min}^{-1}$ &	 $2.85\times 10^{3}\ \mathrm{min}^{-1}$ -- \\
&turnover number& &$2.85\times 10^{5}\ \mathrm{min}^{-1}$\\ \hline
	$\nu$		& lactose$\rightarrow$ allolactose 	&0.468 &--	 		\\
&$\beta$-gal branching& &\\
&fraction& &\\ \hline
	$K_{m,l}$	& $\beta$-gal lactose		& 2.53 mM&	 $0.253\ \mathrm{mM}$ -- $25.3\ \mathrm{mM}$				\\
&Michaelis constant & & \\ \hline
	$\delta$	& $\beta$-gal allolactose		& $2.30\times 10^4\mathrm{min}^{-1}$&	 $2.30\times 10^{3}\ \mathrm{min}^{-1}$ -- \\
&turnover number& &$2.30\times 10^{5}\ \mathrm{min}^{-1}$\\ \hline
	$K_{m,a}$	& $\beta$-gal allolactose	& 1.2 mM&	 $0.12\ \mathrm{mM}$ --$12.0\ \mathrm{mM}$ \\
&Michaelis constant& &\\ \hline
	$\epsilon$	& fully induced 	& 	 $34285\ \mathrm{nM}$&-- 									\\
& $\beta$-gal level & & \\ \hline
	$c$		& basal $\beta$-gal
 level		& 	 $34.3\ \mathrm{nM}$ &--									\\ \hline
	$K_{z}$		& signal level at	& $10^5$ nM&	 $10^{4}\ \mathrm{nM}$ -- $10^6\ \mathrm{nM}$									\\
&half-maximal & & \\
&{\em lac} induction& &\\ \hline
	$n$		& Hill number for		& 	 $2$	&--												\\
&signal-dependent& &\\
&{\em lac} induction& &\\ \hline
\end{tabular}
\end{center}
\label{Tableparams}
\end{table}
\end{singlespace}

\section{Results}
\label{SEC-ANAL}

We first used Eqs.~(\ref{iptg}) to determine how parameter values
control bistability in the steady-state response of {\em lac}
expression to artificial inducers. To detect and characterize
bistability for a given set of parameter values, we solved for $z(l)$
and $l^{*}(l)$ as rational functions of $l$. Bistability in {\em lac}
expression exists when the line describing steady-state levels of $z$
vs. $l^{*}$ adopts a characteristic ``S'' shape, as shown in
Fig.~\ref{OzbTPs}. Within the bistable range of $l^{*}$, the highest
and lowest levels of $z$ are stable steady-state solutions and the
intermediate level of $z$ is an unstable steady-state solution of
Eqs.~(\ref{allo}). The bistable range is defined by the lower
($l^{*}=L$) and upper ($l^{*}=U$) turning points, as illustrated in
Fig.~\ref{OzbTPs}. An analogous signature of bistability can be seen
in examining steady-state levels of $l$ vs. $l^{*}$ (not shown). For a
model with given parameter values, $L$ and $U$ can be located by
finding the roots of either $d l^* / d z$ or $d l^{*} / d l$ using an
eigenvalue solver.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=5.0in]{Fig_2}
\end{center}
\caption{An example of a system from Eqs.~(\ref{iptg}) with the upper
($U$) and lower ($L$) turning points consistent with the results in
\cite{Ozbudak04}.  The parameter values are $\gamma = .0231\
\mathrm{min}^{-1}$, $\alpha = 60\ \mathrm{min}^{-1}$, $K_{z} =
123,285\ \mathrm{nM}$ and $K_{i} = 1,077,217\ \mathrm{nM}$.}
\label{OzbTPs}
\end{figure}

We analyzed Eqs.~(\ref{iptg}) for systems with sets of
parameter values drawn from the ranges in Table~\ref{Tableparams},
taking $\alpha_0=0$, $\phi=0$, and $n=2$. Sets of 100 values each for
$K_i$ and $K_z$ were obtained using logarithmically even sampling over
their allowed ranges. Because the steady-state solutions of
Eqs.~(\ref{iptg}) only depend on $\alpha$ and $\gamma$ through the
ratio $\alpha/\gamma$, rather than sampling $\alpha$ and $\gamma$
individually, we obtained 100 values of $\alpha/\gamma$ using
logarithmically even sampling between the upper and lower bound
computed from Table~\ref{Tableparams}. This sampling scheme yielded
$100\times 100\times 100 = 10^6$ systems with different
values of $(\alpha/\gamma,K_i,K_z)$.

We found that all $10^6$ systems exhibited some degree of bistability
in response to induction by artificial inducers. The dependence of the
range of bistability on model parameters was further analyzed using
two measures that we introduce here: the ratio $U/L$, and the product
$UL$. We used these measures to estimate the percentage of systems for
which bistability might be observable in an experiment like that in
Ref.~\cite{Ozbudak04}. By inspecting the measurement errors in
Ref.~\cite{Ozbudak04}, we estimate that systems with $U/L>1.1$ and
$UL>0.01\ \mu\mathrm{M}^2$ exhibit bistability that is favorable for
experimental observation (i.e., difficult to detect), and that systems
with either $U/L<1.1$ or $UL<0.01\ \mu\mathrm{M}^2$ exhibit
bistability that is unfavorable for experimental observation. Among
systems with parameter values sampled as described above, by these
criteria, we predict that experimental observation of bistability
would be favorable for $65\%$ of systems, and unfavorable for $35\%$
of systems.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=5.0in]{Fig_3}
\end{center}
\caption{Effects of variations in the $\alpha_{0} > 0$ parameter on an
artificially induced system with $\phi = 0\ \mathrm{min}^{-1}$ and
$\alpha_{0} = 10^{-k}\ \mathrm{min}^{-1}$, $k = 0,\ldots, 4$. The
other parameters are given by $n = 2,\ \gamma = .0231\
\mathrm{min}^{-1},\ \epsilon = 34286\ \mathrm{nM},\ c = 34.3\
\mathrm{nM},\ K_{i} = 5 \times 10^{6}\ \mathrm{nM},\ K_{z} = 10^4\
\mathrm{nM}$ and $\alpha = 60\ \mathrm{min}^{-1}$.}
\label{alpha}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=5.0in]{Fig_4}
\end{center}
\caption{Effects of variations in the $\phi > 0$ parameter on an
artificially induced system with $\alpha_{0} = 10^{-4}\
\mathrm{min}^{-1}$ and $\phi = 0$ and $10^{-k}\ \mathrm{min}^{-1}$, $k
= 1,\ldots,4$.  All of the scales are in $\mu$M.  The other parameter
values are as in Figure~\ref{alpha}.}
\label{phi}
\end{figure}

Increasing either $\alpha_{0}$ or $\phi$ above zero tends to reduce or
abolish bistability in artificially induced systems.  As
$\alpha_{0}$ is increased (Fig.~\ref{alpha}), first $U$ begins
shifting to lower values of $l^{*}$, then $L$ begins shifting to
higher values of $l^{*}$, leading to an asymptotic behavior in which
bistability is abolished.  Like changes in $\alpha_0$, as $\phi$ is
increased (Fig.~\ref{phi}), $L$ shifts to higher values of $l^{*}$;
however, by contrast, $U$ does not initially show a significant
change.  As $\phi$ is increased further, the entire induction curve
begins to shift to higher levels of $l^{*}$.

\begin{figure}[htd]
\begin{center}
\includegraphics[width=5.0in]{Fig_5}
\end{center}
\caption{Modeling the width as a sigmoid function of $K_{i}/K_{z}$.  Only a sample of data points are shown.}
\label{WidthKiKz}
\end{figure}

To compare Eqs.~(\ref{iptg}) to the data in Ref.~\cite{Ozbudak04}, we
first selected a subset of systems for which the bistable region is in
the same neighborhood as that in Ref.~\cite{Ozbudak04}: from $3\ \mu
\mathrm{M}$ to $30\ \mu \mathrm{M}$ TMG. Considering this range, out
of the $10^6$ systems sampled, we selected 187,108 systems for which
$L>1\ \mu$M and $U<\ 100\ \mu$M for further analysis. Interestingly, we
found that all of these systems collapse to a single curve when
displayed in the space of $\log_{10}(U/L)$ vs. $\log (K_{i}/K_{z})$
(Fig.~\ref{WidthKiKz}), indicating that $U/L$ can be precisely tuned
using the parameter $X = K_{i}/K_{z}$.  As shown in
Fig.~\ref{WidthKiKz}, the dependence was accurately modeled using the
equation
\begin{eqnarray}
\label{widtheq}
\log_{10} (U/L) & \approx &  \frac{\ (K_{i}/K_{z})^{.93}}{(K_{i}/K_{z})^{.93} + (.27)^{.93}} - \frac{1}{10} \geq 0.
\end{eqnarray}
Next, we found that, at a given value of $X=K_{i}/K_{z}$, without
changing the value of $U/L$, $UL$ could be tuned precisely using the
parameter $Y = K_{i}K_{z}\gamma/\alpha$. As shown in
Fig.~\ref{CenterReg}, this dependence was accurately modeled using the
equation
\begin{eqnarray}
\label{centereq}
\log_{10} (UL) & = & C_{0}(X) + C_{1}(X) \log_{10}(Y).
\end{eqnarray}
Figure~\ref{CenterCoeffs} shows the $X$-dependence of the parameters
$C_0(X)$ and $C_1(X)$, obtained numerically using systems with similar
values of $X$.  For the range of systems considered here, we found
that $C_{0}(X)$ could be fit using a third order polynomial in
$\log_{10}(X)$, and that $C_{1}(X)$ could be taken as a constant.

\begin{figure}[htd]
\begin{center}
\includegraphics[width=5.0in]{Fig_6}
\end{center}
\caption{Using the $C_{0}$ and $C_{1}$ in Figure~\ref{CenterCoeffs}
results in a linear relation between the center $\log_{10}(UL)$ and
$\log_{10}(K_{i}K_{z}\gamma/\alpha)$.}
\label{CenterReg}
\end{figure}

\begin{figure}[htd]
\begin{center}
\includegraphics[width=5.0in]{Fig_7}
\end{center}
\caption{Dependence of the coefficients in the regression model of
Eq.~(\ref{centereq}) on $X=K_i/Kz$.}
\label{CenterCoeffs}
\end{figure}

The above phenomenological results provide a prescription for tuning
the range of bistability exhibited by an artificially induced
system. First, the value of $U/L$ can be specified by choosing a value
of the parameter $X=K_{i}/K_{z}$ using Eq.~(\ref{widtheq}). Then,
using this value of $X$, the value of $UL$ can be specified by
choosing a value of the parameter $Y=K_z K_i \alpha/\gamma$ using
Eq.~(\ref{centereq}) and the empirically determined $C_0(X)$ and
$C_1(X)$ (Fig.~\ref{CenterCoeffs}). We used this prescription to
obtain a family of systems that are consistent with the parameter
values in Table~\ref{Tableparams} and that exhibit a range of
bistability consistent with that observed in Ref. \cite{Ozbudak04},
with $\log_{10}(U/L) \approx .86$ and $\log_{10}(UL) \approx 1.92$. An
example of the steady-state behavior of one such system is illustrated
in Figure~\ref{OzbTPs}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=5.0in]{Fig_8}
\end{center}
\caption{Allolactose system with all of the parameters given by their
nominal values in Table~\ref{Tableparams}.}
\label{allosys}
\end{figure}

We used similar methods to analyze Eqs.~(\ref{allo}) which describe
induction by lactose. No bistability was present in the system with
nominal parameter values from Table~\ref{Tableparams}
(Fig.~\ref{allosys}) with $\phi=0.5$, which is consistent with the
theory of Savageau~\cite{Savageau01} and the Supplementary Material of
Ref.~\cite{Ozbudak04}. However, guided by the results for artificial
inducers in Fig.~\ref{phi}, we examined systems with
$\phi=0$. Although the system with nominal parameter values and
$\phi=0$ did not exhibit bistability, other systems that have
parameter values consistent with the ranges in Table~\ref{Tableparams}
did exhibit bistability. We then located the system that exhibits the
largest values of $U/L$ and $UL$; for this case, $\alpha$, $\beta$,
$\delta$ and $K_{z}$ assume their lowest values in
Table~\ref{Tableparams} while $\gamma$, $K_{m,l}$, $K_{m,a}$ and
$K_{i}$ assume their highest values (Fig.~\ref{Biallosys}). The curve
in Fig.~\ref{Biallosys} illustrating bistability characteristics for
this system closely resembles a similar curve shown in van Hoek \&
Hogeweg~\cite{vanHoek06}, Fig. 2B. Thus, although our model is less
detailed than theirs, it can exhibit comparable steady-state behavior.

To estimate the distribution of systems exhibiting the different
qualitative behaviors, as for the case of artificial inducers, we
analyzed $10^4$ systems with randomly sampled parameter values, all
with $\phi = 0$.  We predict $99.87\%$ of these systems to exhibit no
bistability, $0.05\%$ to exhibit bistability favorable for observation
($U/L > 1.1$ and $UL > 0.01\ \mu \mathrm{M}^2$), and $0.08\%$ to
exhibit bistability that is unfavorable for observation ($U/L < 1.1$
or $UL < 0.01\ \mu \mathrm{M}^2$). However, as observed for
Eqs.~(\ref{iptg}), increasing $\phi$ to even a small fraction of its
nominal value rapidly abolishes bistability for all combinations of
other parameter values in Eqs.~(\ref{allo}) (Fig.~\ref{Allophi}).

\begin{figure}[ht]
\begin{center}
\includegraphics[width=5.0in]{Fig_9}
\end{center}
\caption{Bistability in the $\phi = 0$ lactose-induced system with
$\alpha$, $\beta$, $\delta$ and $K_{z}$ at their lowest values in
Table~\ref{Tableparams} and $\gamma$, $K_{m,l}$, $K_{m,a}$ and $K_{i}$
at their highest values. This is the system that exhibits the largest
values of $U/L$ and $UL$ within the allowed ranges of parameter
values.}
\label{Biallosys}
\end{figure}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=5.0in]{Fig_10}
\end{center}
\caption{This is the same system as in Figure~\ref{Biallosys} with $\phi = 0$ and $10^{-k}\ \mathrm{min}^{-1}$, $k = 1,\ldots,4$.}
\label{Allophi}
\end{figure}

\section{Conclusions}
\label{SEC-CONC}

For the equations describing induction by artificial inducers, we
found that the range of external inducer concentrations over which
systems exhibit bistability is precisely controllable by two rational
combinations of model parameters. By adjusting these parameters, we
were able to demonstrate agreement with the bistable range for TMG
induction from Ref.~\cite{Ozbudak04}. However, in achieving this
agreement, we assumed that permease-dependent efflux of artificial
inducers is negligible ($\phi=0$). We have not found independent
biophysical data to constrain this parameter for artificial inducers,
and therefore predict that it has a value much less than the value of
roughly $0.5$ that has been measured for lactose.

To achieve agreement with the bistable range of roughly $3\ \mu
\mathrm{M}$ to $30\ \mu \mathrm{M}$ in Ref.~\cite{Ozbudak04}, $c$ and
$\epsilon$ in Eqs.~(\ref{iptg}) were tuned to exhibit a 1000-fold
induction of protein expression. While this value is reasonable based
on previous studies, it does disagree with the roughly 100-fold
induction of GFP expression reported in Ref.~\cite{Ozbudak04}. We did
analyze systems with alternative values of $c$ and $\epsilon$ that
yield 100-fold induction; however, none of them exhibited bistable
ranges that agree with the range reported in
Ref.~\cite{Ozbudak04}. Further studies will be required to understand
why Eqs.~(\ref{iptg}) does not simultaneously agree with both the
bistable range and maximal induction of the experimental $lac::gfp$
reporter system. In addition to model refinement, it would be fruitful
to seek systematic differences between expression from chromosomal
{\em lac} and the plasmid-based $lac::gfp$ reporter system used in
Ref.~\cite{Ozbudak04}.

The lack of bistability observed for induction by lactose agrees with
modeling studies concluding that bistabity in {\em lac} expression is
irrelevant to {\em E. coli} in a natural
context~\cite{Savageau01,Savageau02,vanHoek06,vanHoek07}. Thus,
although bistable behavior in {\em lac} is now
well-documented~\cite{Novick57,Cohn59c,Ozbudak04}, because it has only
been experimentally observed using artificial inducers, its relevance
within the natural context of {\em E. coli} is doubtful. Indeed, it is
surprising that the {\em lac} operon has been considered to be a
paradigm of bistability in gene regulation, considering the gaps in
understanding that remain after so many careful experimental and
theoretical studies.

The present results predict that bistable behavior can be promoted by
(1) hindering the kinetics of permease transport ($\alpha$, $K_i$) and
$\beta$-gal catalysis ($\beta$, $\delta$, $K_{m,l}$, $K_{m,a}$); (2)
lowering the required level of allolactose for half-maximal {\em lac}
expression ($K_z$); and (3) accelerating cell growth ($\gamma$). These
predictions suggest genetic targets for engineering {\em E. coli}
strains that exhibit a clear signature of bistability. Experiments to
compare the behavior of such strains with wild-type cells would help
to clarify whether bistability in {\em lac} expression is relevant in
a natural context.

\newpage

\begin{acknowledgments}
\label{SEC-ACK}

Supported by the US Department of Energy through contract
DE-AC52-06NA25396, and grant GM 80216 from the National Institutes of
Health. The collaboration was facilitated by the First q-bio Summer
School on Cellular Information Processing, which was sponsored by the
New Mexico Consortium's Institute for Advanced Studies, and the Center
for Nonlinear Studies at Los Alamos National Laboratory.

\end{acknowledgments}

\newpage


\section{INTRODUCTION}

Dark state or coherent population trapping (CPT), which is null
eigenvalue state of the interaction Hamiltonian, has given rise to
growing interests for past few decades, since its first discovery in
1976~\cite{alzetta1976}. If a three-level $\Lambda$ atomic system is
coherently trapped in dark states, there is no absorption of
incident laser fields even in the presence of resonant transitions,
due to destructive quantum interference between two transition
pathways. Well-known dark state or
CPT~\cite{Arimondo-Orriols,gray,Arimondo,EIT-RMP-2005,CPT-RMP-1998,Lukin1999,quantum-optics,dark-state-polariton2000}
involves extensive applications in many fields such as
electromagnetically induced transparency (EIT)~\cite{harris1997},
lasing without inversion~\cite{lasing1,lasing2,lasing3}, the
enhancement of the refractive index~\cite{enhancement of the
refractive index 1,enhancement of the refractive index 2}, adiabatic
population transfer~\cite{oreg1984,bergmann1989}, subrecoil laser
cooling~\cite{cohen1994}, and atom interferometry~\cite{chu}.
Coherent population transfer using stimulated Raman adiabatic
passage(STIRAP)
technique~\cite{CPT-RMP-1998,annu_rev_phys_chem2001,bergmann1989,Gaubatz1990},
make it robust and efficient to prepare atoms and molecules in a
well defined required state in many fields such as spectroscopy,
collision dynamics, and atomic physics. The Stokes and pump laser
fields in a kind of counterintuitive order under two-photon
resonance and adiabatic evolution carry out complete coherent
population transfer, and the STIRAP technique is insensitive to
variations of laser pulse shape, intensity, and laser frequency. It
is also of crucial importance that a required coherent superposition
state is prepared for atoms and molecules. Fractional STIRAP
technique~\cite{Marte1991,J_Phys_B1999}, where the Stokes laser
arrives before the pump laser and both terminate simultaneously
maintaining a constant ratio of relative amplitudes, can create an
arbitrary required coherent superposition state. The technique is
also insensitive to pulse delay, intensity and frequency of laser
fields. Another adiabatic transfer~\cite{Unanyan1998} using a tripod
linkage by three laser fields in a four-level system can also to
obtain good result.

Here we propose an automated method for adiabatic population
transfer and the preparation of an arbitrary quantum superposition
state in a three-level $\Lambda$ configuration atomic system, using
the oscillating dark states (ODS), where probability amplitudes of
the system state oscillate periodically between two ground levels.
The state to be prepared can be retrieved at a predictable time, and
this method is insensitive to the initial state.

The paper is organized as follows. The next section contains the
discussion of the oscillating dark states (ODS) and evolving
adiabatic conditions, under which the system follows the stable ODS
evolution. In Sec. III, we introduce the master equation for a
practical incoherent atomic system, and simulate the ODS evolution,
which then be used for adiabatic population transfer and the
preparation of an arbitrary quantum superposition state. We give
some discussion and summarize our conclusions in Sec. IV.

\section{OSCILLATING DARK STATES}

The three-level $\Lambda$ configuration atomic system is shown in
Fig.~\ref{final-big2}. The ODS is implemented with two pairs of
classical detuning laser fields coupling two lower ground levels
$|1\rangle$ , $|2\rangle$ to a single upper level $|3\rangle$,
respectively. We assume that the $|1\rangle$-$|2\rangle$ transition
is always dipole-forbidden. The time-dependent interaction
Hamiltonian in the interaction picture that describes the atom-laser
coupling within the dipole and rotating wave approximation (RWA) in
the rotating frame reads ($\hbar$=$1$),
\begin{equation}
H_{int}=\left[
\begin{array}{c}~0~~~~~0~~~~P^*\\~0~~~~~\Delta'~~Q^*\\P~~~~Q~~~~\Delta
\end{array} \right]\;, \label{eq:interaction}
\end{equation}
where
$P$=$-ie^{-i\phi_{12}}\Omega_{12}\sin(\frac{\Delta_1-\Delta_2}{2}t)$,
$Q$=$-e^{-i\phi_{34}}\Omega_{34}\cos(\frac{\Delta_3-\Delta_4}{2}t)$,
$\Delta$=$\frac{1}{2}(\Delta_1+\Delta_2)$, and
$\Delta'$=$\frac{1}{2}[(\Delta_1+\Delta_2)-(\Delta_3+\Delta_4)]$.
$\Delta_i$=$\omega_{31}$-$\omega_i$($i$=$1,2$), and
$\Delta_j$=$\omega_{32}$-$\omega_j$($j$=$3,4$) are the detunings of
two pairs of laser angular frequencies $\omega_i$($i$=$1,2$) and
$\omega_j$($j$=$3,4$) from the corresponding atomic transitions
$\omega_{31}$ and $\omega_{32}$. Real $\Omega_i$ and $\phi_i$
($i$=$1,2 ; 3,4$) are corresponding Rabi frequencies and phases of
two pairs of laser fields, and we assume that
$\Omega_1$=$\Omega_2$=$\Omega_{12}$,
$\Omega_3$=$\Omega_4$=$\Omega_{34}$,
$\phi_1$=$\phi_2$+$\pi$=$\phi_{12}$, and
$\phi_3$=$\phi_4$=$\phi_{34}$, with the difference
$\Delta\phi$=$\phi_{12}$-$\phi_{34}$.
\begin{figure}
\centering
\includegraphics[scale=0.20]{final-big2.eps}%
\caption{\label{final-big2} Two pairs of classical laser fields
$\Omega_1$ , $\Omega_2$ and $\Omega_3$ , $\Omega_4$ couple two lower
ground levels $|1\rangle$ , $|2\rangle$ to a single upper level
$|3\rangle$ , respectively, with their corresponding detunings
$\Delta_1$ , $\Delta_2$ and $\Delta_3$ , $\Delta_4$, and their fixed
phases: $\phi_1$=$\phi_2$+$\pi$, $\phi_3$=$\phi_4$, in a three-level
$\Lambda$ configuration atomic system.}
\end{figure}
The eigenstates of the interaction
Hamiltonian~(\ref{eq:interaction}), which can be generally described
in terms of the \textquotedblleft mixing angles\textquotedblright
$\theta$ and $\varphi$ , are given as
\begin{equation}
|a_+\rangle=-ie^{i\Delta\phi}\sin\theta\sin\varphi|1\rangle+\cos\theta\sin\varphi|2\rangle-\cos\varphi
e^{-i\phi_{34}}|3\rangle\;, \label{eq:a+}
\end{equation}
\begin{equation}
|a_0\rangle=\cos\theta|1\rangle-i\sin\theta
e^{-i\Delta\phi}|2\rangle\;,
\label{eq:a0}
\end{equation}
\begin{equation}
|a_-\rangle=-ie^{i\Delta\phi}\sin\theta\cos\varphi|1\rangle+\cos\theta\cos\varphi|2\rangle+\sin\varphi
e^{-i\phi_{34}}|3\rangle\;, \label{eq:a-}
\end{equation}
where the mixing angles that dependent upon time $t$ and Rabi
frequencies are written as
\begin{equation}
\tan\theta=\frac{\Omega_{12}\sin(\frac{\Delta_1-\Delta_2}{2}t)}
{\Omega_{34}\cos(\frac{\Delta_3-\Delta_4}{2}t)}\;, \label{eq:tan}
\end{equation}
\begin{equation}
\tan2\varphi=\frac{2\sqrt{|P|^2+|Q|^2}}{\Delta}\;,
\end{equation}
when we meet the condition $\Delta'=0$, i.e.
$(\Delta_1+\Delta_2)$=$(\Delta_3+\Delta_4)$. The eigenstate
$|a_0\rangle$ has no contribution from $|3\rangle$, and this ODS
corresponds the null eigenvalue of interaction
Hamiltonian~(\ref{eq:interaction}). If a system is in state
$|a_0\rangle$, there is no possibility of excitation to $|3\rangle$
and subsequent spontaneous emission, thus the atom is effectively
decoupled from laser fields. The eigenvalues of the pair of states
$|a_\pm\rangle$ which contain a component of all three bare atomic
states, are shifted up and down by an amount $\lambda_\pm$,
\begin{equation}
\lambda_\pm=\frac{1}{2}[\Delta\pm \sqrt{\Delta^2+4|P|^2+4|Q|^2} ]\;.
\label{eq:lambda-pm}
\end{equation}
The probability amplitudes of the system state can oscillate
periodically between two ground levels $|1\rangle$ and $|2\rangle$
upon time $t$, when the system is driven into the ODS by laser
fields with relatively fixed detunings and phases of four laser
fields.

\par Evolving adiabatic conditions, under which the system avoids
diabatic coupling to bright states and is always in some dark states
if it is initially in a certain one, is considered, since the ODS
follows the evolution of a time-dependent Hamiltonian
(\ref{eq:interaction}). The change rate of the mixing angle $\theta$
must be small compared with the separation of the corresponding
eigenvalues~\cite{Messiah1962,Gaubatz1990},
\begin{equation}
|\frac{\textit{d}}{\textit{dt}}\theta|^2\ll|\lambda_0-\lambda_\pm|^2\;.
\label{eq:totalconditon}
\end{equation}
Let $\Omega_{12}$=$\Omega_{34}$=$\Omega$ (fixing Rabi frequencies
when the system evolves) and
$\Delta_1$-$\Delta_2$=$\Delta_3$-$\Delta_4$=2$\delta$, thus we can
obtain evolving adiabatic conditions according to
(\ref{eq:totalconditon}) as
\begin{equation} |\delta|^2\ll \frac{1}{4}|\Delta\pm
\sqrt{\Delta^2+4\Omega^2}|^2\;, \label{eq:Evolving adiabatic
conditions}
\end{equation}
which changes into $|\delta|\ll$$\Omega$, when $\Delta=0$, i.e.
$\Delta_1$=-$\Delta_2$ and $\Delta_3$=-$\Delta_4$. Evolving
adiabatic conditions for the ODS ensure that the values of all
parameters can be established in advance, because the evolution
process is automated by laser-atom system itself, in contrast to
ordinary dark states~\cite{quantum-optics,bergmann1989}, where we
manually with experimental instruments change relative values of two
Rabi frequencies to attain different dark states and to utilize
STIRAP technique.

\section{ADIABATIC POPULATION TRANSFER AND QUANTUM STATE PREPARATION}

A process for adiabatic population transfer and quantum state
preparation is as follows: first, the initial instantaneous
eigenstate of the whole system, which is determined by chosen
interaction Hamiltonian, must be the atomic initial state, and here
we assume that the atomic initial state is $|1\rangle$. Then we
concurrently and rapidly upload two pairs of laser fields as shown
in Fig.~\ref{final-big2}, and the system periodically evolves
subsequently. Rapid unloading all laser fields simultaneously
completes a perfect population transfer from $|1\rangle$ to
$|2\rangle$ after ($\frac{1}{4}$+$\frac{n}{2}$) ($n$=$0,1,2,...$,
integral number) oscillation periods($T$=$2\pi/\delta$) if we want
to make an adiabatic population transfer, and retrieves a required
quantum superposition state after time of ($t_0$+$nT$)($t_0$ is
evolution time when the system firstly evolves from $|1\rangle$ to
our required state, which can be easily predicted from (\ref{eq:a0})
where a fixed $\Delta\phi$ should be calculated in advance) if we
use the same experimental setup for the preparation of an arbitrary
quantum superposition state.

When we upload/unload two pairs of laser fields, uploading/unloading
adiabatic conditions are necessary besides evolving adiabatic
conditions~(\ref{eq:Evolving adiabatic conditions}), in order to
keep the least diabatic coupling between bright states and dark
states of the system. We assume a concurrent increase/decrease of
four laser fields in a short time $\tau$($\tau$$\ll$$T$), i.e.
$\Omega_{1}$=$\Omega_{2}$=$\Omega_{3}$=$\Omega_{4}$=$\Omega$ and
$\frac{\textit{d}}{\textit{dt}}$$\Omega_{1}$=
$\frac{\textit{d}}{\textit{dt}}$$\Omega_{2}$=
$\frac{\textit{d}}{\textit{dt}}$$\Omega_{3}$=
$\frac{\textit{d}}{\textit{dt}}$$\Omega_{4}$. Thus, we can obtain
the same result as evolving adiabatic conditions in
(\ref{eq:Evolving adiabatic conditions}) for uploading/unloading
adiabatic conditions, according to (\ref{eq:totalconditon}).

Time-dependent interaction Hamiltonian in interaction picture that
describes atom-laser coupling can be written as
\begin{eqnarray}
H'_{int}=-\frac{1}{2}[\Omega_{1}e^{-i\phi_1}e^{i\Delta_1t}\sigma_{31}+\Omega_{2}e^{-i\phi_2}e^{i\Delta_2t}\sigma_{31}\nonumber\\
+\Omega_{3}e^{-i\phi_3}e^{i\Delta_3t}\sigma_{32}+\Omega_{4}e^{-i\phi_4}e^{i\Delta_4t}\sigma_{32}+H.c.]\;,
\label{eq:H-int'}
\end{eqnarray}
where $\sigma_{ij}$=$|i\rangle\langle j|$ ($i,j$=$1,2,3$) are atomic
projection operators. A master
equation~\cite{EIT-RMP-2005,cohen1992}, using density matrix $\rho$,
\begin{eqnarray}
\frac{\textit{d}\rho}{\textit{dt}}=-i[H'_{int},\rho]+\frac{\Gamma_{31}}{2}[2\sigma_{13}\rho\sigma_{31}-\sigma_{33}\rho-\rho\sigma_{33}]\nonumber\\
+\frac{\Gamma_{32}}{2}[2\sigma_{23}\rho\sigma_{32}-\sigma_{33}\rho-\rho\sigma_{33}]\nonumber\\
+\frac{\gamma_{3\textit{deph}}}{2}[2\sigma_{33}\rho\sigma_{33}-\sigma_{33}\rho-\rho\sigma_{33}]\nonumber\\
+\frac{\gamma_{2\textit{deph}}}{2}[2\sigma_{22}\rho\sigma_{22}-\sigma_{22}\rho-\rho\sigma_{22}]\nonumber\\
+\frac{\Gamma_{21}}{2}[2\sigma_{12}\rho\sigma_{21}-\sigma_{22}\rho-\rho\sigma_{22}]\;,
\label{eq:master-equation}
\end{eqnarray}
is required, because the decoherence of a practical atomic system
must be considered. $\Gamma_{31}$ and $\Gamma_{32}$ represent rates
of spontaneous emission from level $|3\rangle$ to $|1\rangle$ and
$|2\rangle$, respectively, and rates $\gamma_{2\textit{deph}}$ and
$\gamma_{3\textit{deph}}$ describe energy-conserving dephasing
processes. $\Gamma_{21}$ is longitudinal relaxation from level
$|2\rangle$ to $|1\rangle$, and it is much small compared with
$\gamma_{2\textit{deph}}$. We can define coherence decay rates as
$\gamma_{31}$=$\Gamma_{31}$+$\Gamma_{32}$+$\gamma_{3\textit{deph}}$,
$\gamma_{21}$=$\Gamma_{21}$+$\gamma_{2\textit{deph}}$, and estimate
$\Gamma_{21}$ at $1/10$ of $\gamma_{2\textit{deph}}$ in the
following analysis.

\begin{figure}
\includegraphics[scale=0.5]{2.eps}
\caption{\label{fig:2} Probability amplitudes evolution of initial
state $|1\rangle$ finds expression in the time-dependent populations
of three levels $\rho_{11}$, $\rho_{22}$, $\rho_{33}$, and coherence
term between two ground levels $|\rho_{21}|^2$ with
$\Delta_1=\Delta_3=0.3$, $\Delta_2=\Delta_4=0.2$,
$\gamma_{2\textit{deph}}=0.02$, $\Gamma_{21}=0.002$, for different
case: (a). $\Omega=2$, (b). $\Omega=0.2$, (c). $\Omega=0.08$,
respectively. Randomly chosen $\Delta \phi$ (adopted for different
required superposition state, but having no effect on coherent
population transfer process) gives the same periodical oscillation.
All parameters are in the units of $\gamma_{31}$, and total
evolution time $t$=4$T$ is marked.}
\end{figure}

The evolution of the ODS is simulated numerically as shown in
Fig.~\ref{fig:2}, when the initial system is in $|1\rangle$ and a
fixed $\Delta\phi$ should be calculated in advance (for a required
arbitrary final state), and probability amplitudes evolution of the
ODS with different Rabi frequencies $\Omega$ is analyzed. The
populations of two ground lower levels, $\rho_{11}$ and $\rho_{22}$,
oscillate periodically after we applying interaction Hamiltonian
(\ref{eq:H-int'}), and the population of upper level $\rho_{33}$
stays almost zero in total evolution stage, when evolving adiabatic
conditions~(\ref{eq:Evolving adiabatic conditions}) are well
satisfied for the case in Fig.~\ref{fig:2}(a). The same periodical
oscillation of the coherence term $|\rho_{21}|^2$ between two ground
levels as $\rho_{11}$ and $\rho_{22}$, and
$|\rho_{21}|^2$$\approx$$\rho_{11}$$\cdot$$\rho_{22}$, imply a
coherent evolution of quantum superposition states, i.e. the ODS. In
case (b), Rabi frequency $\Omega$ is decreased, and evolving
adiabatic conditions are not well met, which results in a slight
occupation of population $\rho_{33}$, but the oscillations of
population $\rho_{11}$, $\rho_{22}$ and of coherence term
$|\rho_{21}|^2$ still dominate the whole evolution process. When we
decrease $\Omega$ further, laser fields are too weak to make upper
level $|3\rangle$ occupied despite the violation of evolving
adiabatic conditions, and the coherence term of density matrix
$|\rho_{21}|^2$ keeps $0$ all the time in case (c), which implies
there is no coherent evolution of the system state.

From case (c) to (a), the value of Rabi frequencies $\Omega$ of
laser fields increased little by little, it can be seen that
coherent reconstruction mechanism of laser fields becomes superior
gradually to the decoherence effect of the system in case (c) and
(b), and entirely dominates periodical oscillation of the system in
case (a). A complete population transfer between two lower ground
levels $|1\rangle$ and $|2\rangle$, and the preparation of an
arbitrary quantum superposition state in an automatical manner,
using a coherent evolution process of the ODS in case (a), is
reliable. Here upper level $|3\rangle$ population is much little and
the oscillation of probability amplitudes of two ground levels
$|1\rangle$ and $|2\rangle$ is periodical, since the decoherence
between two ground levels $|1\rangle$ and $|2\rangle$ can be greatly
suppressed, owing to coherent reconstruction contribution from laser
fields. It can be seen that there are many optional retrieval time,
with retrieval interval $T/2$ for an adiabatic population transfer
process and with interval $T$ for the preparation of an arbitrary
quantum superposition state.

With the contribution from applied laser fields, the state of the
system evolves periodically and stably despite decay factor
$\gamma_{21}\cdot t$ expressed in coherence term of the system
density matrix $|\rho_{21}|^2$. When two pairs of laser fields are
applied, we can also consider the long time behavior of state
evolution of the system and analyze the fidelity of initial state
$|1\rangle$ and retrieved states of the system after different
evolution time $t$ in Fig.~\ref{fig:t_VS_fidelity}, where only
integral period evolution ($t$=$nT$) is observed. $\gamma_{21}$ is
fixed at $0.022$ in the units of $\gamma_{31}$ (reference to $D_1$
line of $^{87}Rubidium$ with
$\gamma_{31}$$\sim$$36.10\times10^6$$s^{-1}$ and estimated
$\gamma_{2\textit{deph}}$$=$$0.02$$\gamma_{31}$), and different
total evolution times $t$ is given. It can be seen that, when total
evolution time $t$ is up to long enough $1000$$T$, where
$\gamma_{21}\cdot t$$>$$2000$, high fidelity above $0.95$ is
surprisingly available, which implies almost a coherence preserved
and pure periodical evolution of the system. Fidelity $\textit{F}$
decreases slightly in the first one $T$, during which a certain
complicated amplitude decay of the ODS experiences, and stays
unchangeably after it, which also can be found in the very beginning
of evolution process in Fig.~\ref{fig:2}. High fidelity $\textit{F}$
mainly depends on relatively large value of $\Omega$ and always can
be obtained by controlling laser fields.

\begin{figure}
\includegraphics[scale=0.34]{3.eps}
\caption{\label{fig:t_VS_fidelity} Fidelity of initial state
$|1\rangle$ and retrieved states of the system after different
evolution time of the system (only integral period evolution
observed, and up to a maximum $1000T$ here), with $\Omega=2$,
$\Delta_1=\Delta_3=0.3$, $\Delta_2=\Delta_4=0.2$,
$\gamma_{21}=0.022$. After the first one $T$ (see inset), High
fidelity above $0.95$ is available and keep stable at all times,
even if $\gamma_{21}\cdot t$$\gg$1, where different initial states
of system almost give the same result.}
\end{figure}

\section{DISCUSSION AND CONCLUSION}

The stability of frequency difference $2\delta$, $\pi$ phase
difference, and $\Delta\phi$ of each pair of laser fields determine
mainly the precision of the ODS evolution. The population transfer
efficiency and an ideal required superposition state mainly depend
on the exact retrieval time that can be easily calculated according
to~(\ref{eq:a0}). Two pairs of laser fields should be
uploaded/unloaded spontaneously instead of counterintuitive and
delayed uploading of one pair Stokes and pump laser fields. A
concurrent uploading/unloading time $\tau$$\sim$$10^{-8}$$s$ (a
reasonable $0.01T$ in our analysis which is short enough to avoid a
certain decay of $\rho$ at this stage) can be carried out using an
optical switch, provided all the parameters are chose as ours, which
can be tuned to change the oscillation period $T$. Laser fields may
not be monochromatic, but the stability of relative phases of two
pairs of laser fields, which can be easily achieved by
acousto-optical modulation(AOM) device, can also guarantee almost
the same result, despite a slow collective phase drift. A chosen
$\delta$ whose value is larger than the linewidth of laser fields
can avoid undesirable interference of laser fields and optical
pumping. The ODS method is immune to Doppler frequency shift, pulse
shape of laser fields, and concurrent fluctuation of Rabi
frequencies $\Omega_i$ ($i$=$1,2 ; 3,4$). Well developed laser
techniques can meet all above requirements to obtain a perfect
adiabatic population transfer and the preparation of an arbitrary
quantum superposition state in atomic system. The greatly suppressed
decoherence of the system according to our numerical calculation
suggests the stability of the ODS evolution. Numerous optional
retrieval time owing to the property of the ODS evolution is
available, with retrieval interval $T/2$ for an adiabatic population
transfer process and with interval $T$ for the preparation of an
arbitrary quantum superposition state, respectively. The transfer
process is insensitive to the initial state of the atomic
system~\cite{ottaviani-simon}, i.e., we can carry out complete
coherent population transfer and obtain a required state at right
retrieval time based on precise ODS evolution of atom-laser system,
whatever the initial state of atomic system is, which is totally
different from the STIRAP and fractional STIRAP technique. Different
initial states of the system besides $|1\rangle$ here, to which the
corresponding initial Hamiltonian of the system can be chosen to
match the initial instantaneous eigenstate, give almost the same
good result. Optical pumping is not necessary any longer in contrast
to conventional STIRAP technique, since it is not essential that
initial state of the system must be any single lower state
$|1\rangle$ or $|2\rangle$ when the automated ODS method is adopted.
The ODS can be implemented in a three-level atomic system by using
at least three laser fields in an appropriate configuration, but we
use two pairs of laser fields in a symmetric scheme, in order to
obtain complete oscillation of probability amplitudes between two
ground levels. The ODS is similar to the famous Rabi resonance in a
two-level atomic system driven by one resonant laser
field~\cite{quantum-optics}, where we can not find the same
preserved coherence.

In conclusion, we propose an automated method using the ODS for a
perfect adiabatic population transfer and the preparation of an
arbitrary quantum superposition state in atomic system. The
preserved coherence, the automated manner, numerous optional
retrieval time, an arbitrary initial state of system, make the ODS a
potential method for physics process in atomic physics.

\begin{center}
\textbf{ACKNOWLEDGMENTS}
\end{center}

We thank Carlo Ottaviani and Christoph Simon for helpful discussion.
This work was funded by the National Fundamental Research Program
(Grant No. 2006CB921907), the National Natural Science Foundation of
China (Grant No. 10674127 and No. 60621064), the Innovation Funds
from the Chinese Academy of Sciences, Program for NCET and
International cooperation program from CAS.


\section{Introduction}

The maturation of laser cooling and trapping techniques has enabled
many recent advances with elements ranging over much of the periodic
table. Narrow intercombination lines in Sr provide ultra-precise
atomic frequency standards \cite{Takamoto:2003a,Boyd:2006a}; Yb has
also been proposed as a frequency standard
\cite{porsev:2004a:021403,hoyt:2005a:083003}, with quantum
degeneracy recently reached in both bosonic \cite{Takasu:2003a} and
fermionic isotopes \cite{fukuhara:2007a:030401}; a Ra
magneto-optical trap (MOT) promises new tests of fundamental
symmetries, while also exhibiting an interesting blackbody radiation
repumping mechanism \cite{guest:2007a:093001}; Cr has a large
ground-state magnetic moment (6$\mu_B$), which is expected to give
long-range order to a Bose-Einstein condensate
\cite{griesmaier:2005a:160401}. In an extension of laser cooling
techniques to a highly magnetic rare earth atom, one of us recently
reported an Er MOT where a novel recycling mechanism allows trapping
despite significant leakage into energy levels that are dark to the
cooling light \cite{mcclelland:2006a:143005}.

Along with the other rare-earth elements, erbium offers an exciting
combination of favorable properties for future cold atomic physics
experiments. All six stable isotopes have been trapped, with
significant populations in three bosonic ($^{166}$Er, $^{168}$Er,
$^{170}$Er) and one fermionic ($^{167}$Er, nuclear spin $I=7/2$)
species. Er exhibits a variety of accessible transitions that are
useful for spectroscopy or laser cooling, including a broad (36~MHz)
line at 401 nm, a narrow (8~kHz) line at 841 nm, and an ultra-narrow
(2~Hz) line in the telecommunication spectrum at 1299 nm
\cite{Martin:1978a,Ban:2005a,mcclelland:2006b:064502}. It has a very
large ground-state magnetic moment (7$\mu_B$), which produces a
strong magnetic dipole interaction at low temperatures. Its closed
$5s^2$ and $6s^2$ electron shells shield the $4f$ core electrons,
resulting in novel scattering properties while also suppressing
undesirable Zeeman relaxation rates \cite{hancox:2004a}.
Additionally, a deterministic ion implantation source based on an Er
MOT \cite{hanssen:063416} could provide a single-ion analog of a
solid-state laser while dramatically enhancing the prospects for
scalability in recent rare-earth quantum computing proposals
\cite{wesenberg:2007a:012304}. However, none of these avenues can be
pursued without first understanding the laser cooling and trapping
properties of this new system.

The primary result of this paper is an exhibition of cooling and
trapping properties, both magnetic and optical, which bolster the
growing list of erbium's desirable features for future cold atomic
physics. First, we show that efficient sub-Doppler cooling occurs in
our MOT where temperatures reach as low as 100~$\mu$K, nearly an
order of magnitude below the Doppler temperature ($T_D=860~\mu$K).
Second, we point out that for the 401~nm transition, the Land\'{e}
$g$ factors describing the Zeeman shift of the ground- and
excited-state magnetic sublevels are nearly equal. This near
degeneracy of Land\'{e} $g$ should soften the usual constraint that
$\sigma^+\sigma^-$ polarization gradient cooling can only occur in
negligible magnetic fields and may have interesting applications for
atomic beam collimation or atomic waveguiding along magnetic field
lines. Finally, we show that a magnetically trapped atom fraction in
our MOT is unexpectedly cold, with temperatures falling below 25
$\mu$K, lower than expectations based on simple thermodynamic
considerations. We believe these magneto-optical properties, coupled
with the unique atomic properties of the rare-earth elements, make
Er an exciting candidate for a wide range of future experiments.

\section{Sub-Doppler cooling in the MOT}
\begin{figure}[t]\includegraphics*[width=3.2in]{LevelDiagram}
\caption{\label{Fig:LevelDiagram}(Color online) Schematic level
diagram of the Er MOT and magnetic trap
\cite{mcclelland:2006a:143005}. The excited state of the 401~nm
transition decays into magnetically trapped metastable states at a
rate $R_1\approx 1700$~s$^{-1}$. Relaxation to the ground state at
rate $R_2$ and trap loss at rate $R_{\mathrm{Loss}}$ both occur with
rates of 4~s$^{-1}$ to 5~s$^{-1}$. $R_{\mathrm{Load}}$ varies
depending on MOT conditions.}\end{figure}

Our apparatus was described in \cite{mcclelland:2006a:143005}. An
effusive atomic beam is produced in the horizontal ($x$) direction
by heating a sample of Er to $1200^\circ$C in a crucible with a 1 mm
aperture. The beam is subsequently decelerated in a $\sigma^-$
Zeeman slower \cite{Barrett:1991a:PhysRevLett.67.3483} then loaded
into a MOT operating on the strong ($\Gamma/2\pi=36$~MHz) transition
at 401 nm. Cooling and trapping light is derived from a
frequency-doubled Ti:sapphire laser. The Zeeman slower beam, MOT
beams, and a probe beam are derived from the same laser, with
frequency offsets and relative intensities controlled by
acousto-optic modulators. All laser frequencies are calibrated and
locked to a fluorescence probe on the atomic beam, with an overall
uncertainty of $\pm 2.8$~MHz \footnote{Unless otherwise noted, all
uncertainties represent one standard deviation combined random and
systematic uncertainty}. The quadrupole magnetic field at the trap
position is the sum of fields from our MOT coils together with the
Zeeman slower and compensation coils, with the largest gradient
along the direction of gravity ($z$). Typical field gradients are
$(dB_x,dB_y,dB_z) = (0.20,0.11,-0.33)$~T/m with $\pm5\%$
uncertainty. All of the measurements described here were performed
on $^{166}$Er, with trapped atom populations between $10^4$ and
$2\times 10^5$.

\begin{figure}
\includegraphics*[width=3.2in]{MOTTvsDetuning_Origin}
\caption{\label{Fig:MOTTVsDetuning}(Color online) MOT temperature
vs. detuning for $I=14.3$~mW/cm$^2$. The upper curve is expected for
Doppler cooling a two-level atom, while the lower curve is a fit to
Eq.~\eqref{Eq:SubDopplerScaling}. The atom number ranged from a
minimum value of approximately $10^4$ at $\delta/2\pi = -5~$MHz to a
maximum of $1.3\times 10^5$ at $\delta/2\pi = -30~$MHz
$(-0.8\Gamma)$. (Inset) MOT temperature vs. trapping beam intensity
at a fixed detuning $\delta=-0.7\Gamma$. The shaded gray area is a
fit to Eq.~\eqref{Eq:SubDopplerScaling} together with the spread
resulting from uncertainty in $\delta$. Error bars are described in
the text.}\end{figure}

Er exhibits a complex level structure, with 110 electronic states
lying between the ground and excited state of the 401 nm
laser-cooling transition. As shown in Fig.~\ref{Fig:LevelDiagram},
the MOT remains functional because decay of excited atoms into
metastable intermediate states does not necessarily result in the
loss of those atoms. On the contrary, a large fraction of these
metastable, dark-state atoms remain magnetically trapped and
eventually relax back to the ground state. We begin by studying the
temperature of the MOT, which we measure by observing the ballistic
expansion of the atom cloud after all optical and magnetic fields
are extinguished. The cloud is imaged by pulsing the MOT beams for
200~$\mu$s and recording the resulting fluorescence image (in the
$xz$ plane) using a CCD camera. Fig.~\ref{Fig:MOTTVsDetuning} shows
the measured temperature $T$ as a function of the detuning of the
MOT beams, $\delta$. The error bars represent uncertainty estimates
from variation in the expansion velocity determined independently
along the $x$ and $z$ directions. The upper solid curve shows the
prediction of simple Doppler theory for a two-level atom
\cite{Metcalf:1999}, while the lower curve is a fit to a
characteristic sub-Doppler scaling law \cite{Drewsen:1994a,Xu:2003a}
\begin{equation}\label{Eq:SubDopplerScaling} T = T_0+C_{\sigma^+\sigma^-} \frac{ \hbar
\Gamma}{2 k_B}\left(\frac{\Gamma}{|\delta|}\right) \frac{I}{I_s}
\end{equation}where $I$ is the total intensity at the MOT
position and $I_s =72$~mW/cm$^2$ is the saturation intensity. The
resulting fit gives $C_{\sigma^+\sigma^-}~=~0.38\pm0.02$ and a
predicted minimum temperature $T_0=(77\pm10)~\mu$K in the
zero-intensity limit. We also varied the intensity at fixed detuning
$\delta=-0.7\Gamma$ (inset of Fig.~\ref{Fig:MOTTVsDetuning}), to
find $C_{\sigma^+\sigma^-}=0.15\pm0.08$ and $T_0=(125\pm21)\mu$K.
The difference in fit parameters between the data sets indicates
that the temperature may have a more complicated dependence on
intensity than the simple scaling law of
Eq.~\eqref{Eq:SubDopplerScaling}. This is not too surprising, since
in our trap, intensity-dependent rates of optical pumping into
magnetically trapped dark states may strongly affect the MOT
temperature (we show below that the magnetically trapped atom
fraction exhibits anomalous temperature behavior). Nevertheless, our
results show qualitatively and quantitatively that strong, efficient
sub-Doppler cooling occurs in the MOT. We emphasize that this is a
single-stage process (see also Ref.~\cite{Xu:2003a}), where
sub-Doppler cooling occurs in the MOT without an additional
``optical molasses" period of polarization-gradient cooling in a
nulled magnetic field.

\begin{figure}[t]\includegraphics*[width=3.2in]{RbVsEr} \caption{\label{Fig:SubDopplerCalcs}
(Color online) Semiclassical force $F(p)$ (in units of $\hbar
k\Gamma$), momentum diffusion coefficient $D(p)$ (in units of
$\hbar^2k^2\Gamma$), and equilibrium probability distribution $W(p)$
(arbitrary units) for $^{87}$Rb (dashed lines, red) and $^{166}$Er
(solid lines, blue) in a $\sigma^+\sigma^-$ optical molasses in a
longitudinal magnetic field $B=1$~mT. For both cases, the saturation
parameter $s=2$ and detuning $\delta=-2.5\Gamma$.}
\end{figure}

As discussed in Refs.~\cite{Walhout:1992a,Werner:1992a},
$\sigma^+\sigma^-$ polarization-gradient cooling typically breaks
down in a longitudinal magnetic field because the ``locking
velocity," where atomic-motion induced frequency shifts are balanced
by Zeeman shifts, is generally different for sub-Doppler and Doppler
cooling mechanisms. This difference arises because Doppler
mechanisms rely on absorption and spontaneous emission, so Zeeman
shifts of both the ground and excited states affect the cooling
efficiency; on the other hand, sub-Doppler mechanisms rely on
stimulated processes between ground state sub-levels, and thus
depend primarily on Zeeman shifts in the ground state. Because the
Land\'{e} $g$ factors of the ground and excited states of a cooling
transition are generally unequal, there is a critical magnetic field
value at which the two mechanisms ``unlock," with the larger capture
range of the Doppler cooling force dominating beyond this point
\cite{Walhout:1992a,Werner:1992a}. For most laser-cooled atoms, this
critical magnetic field value is around 0.1~mT. In Xe, for example,
the difference in $g$ factors between the ground and excited state
is $\Delta g_{eg}=-0.17$ and polarization-gradient cooling was
observed to fail at 0.2 mT~\cite{Walhout:1996a}. On the other hand,
the $g$ factor difference for the 401 nm transition in Er is $\Delta
g_{eg} = -0.004$ \cite{Martin:1978a}, suggesting that
$\sigma^+\sigma^-$ polarization-gradient cooling should remain
effective in longitudinal magnetic fields as large as 10~mT.

To investigate the possible effect of $g$ factor degeneracy on Er
sub-Doppler cooling in a magnetic field, we calculated the
semiclassical force $F(p)$, momentum diffusion coefficient $D(p)$,
and steady-state momentum distribution $W(p)$ for an atom with
momentum $\hbar k p$ in a one-dimensional $\sigma^+\sigma^-$ optical
molasses in a longitudinal magnetic field $B$. Our calculations
follow the method of Ref.~\cite{Walhout:1992a}, in which $p$ is
treated as a classical variable and a Fokker-Planck equation is
derived by expanding the optical Bloch equations to second order in
$\hbar k$. Results are shown in Fig.~\ref{Fig:SubDopplerCalcs} for
the 401~nm line of $^{166}$Er and for the D$_2$ line of $^{87}$Rb
($\Delta g_{eg}=-0.17$) in a magnetic field $B=1$~mT. Because of the
complicated interplay between optical coherences, Zeeman coherences,
and atomic momentum, the distributions $W(p)$ are clearly not
Gaussian. As a result, a temperature does not exist in the strict
sense, but we may still take the variances of $W(p)$ as a measure of
the momentum spread. These variances correspond to center-of-mass
``temperatures" of $25~\mu$K and $3.6$~mK for $^{166}$Er and
$^{87}$Rb, respectively. The same calculation gives respective
temperatures of $25~\mu$K and $17~\mu$K at zero magnetic field
$B=0$. For this simplified one-dimensional case, the $g$-factor
difference results in a temperature increase of more than 2 orders
of magnitude for $^{87}$Rb, with no change in temperature for
$^{166}$Er, at a magnetic field value of $B=1$~mT (10~G). Note that
this value of the magnetic field is attained in our setup at a
distance of a few millimeters from the MOT position, within the beam
waist of our trapping lasers. The near degeneracy of Er $g$-factors
may contribute to capture and loading processes in our MOT and
should also have interesting applications for laser cooling Er in
moderately large magnetic fields.

\begin{figure}[t]
  \centering\includegraphics*[width=3.2in]{BTrapImage2}\caption{\label{Fig:BTrapImage}
  (Color online) CCD image of magnetically trapped atoms after 0.3~s in the magnetic trap. The marginal
  distributions and fits to $p_x(x)$ and $p_z(z)$ are shown at bottom and left, respectively.
  For this image, the fit parameters are $\bar{x}=150~\mu$m, $\bar{z}=105~\mu$m,
  and $\bar{g}=0.23$, corresponding to an effective magnetic moment
  $\bar{\mu}=3.9~\mu_B$ and temperature $T=(44\pm7)~\mu$K.}
\end{figure}
\section{Magnetic Trapping}

\begin{figure}[t]
  \includegraphics*[width=3.2in]{MOTTvsBTrapT_Origin}\caption{\label{Fig:MOTTvsBTrapT}
  (Color online) Magnetic trap temperature vs. MOT temperature. The MOT temperature was
  varied by changing the detuning $\delta$.
  Blue circles and red squares represent data from
  separate experimental runs at different Zeeman slower and MOT beam intensities.
  Simple thermodynamic arguments predict that all
  points should lie above the dashed line where $T_B\geq T_M/3$.}
\end{figure}
We now consider the fraction of ground-state atoms that remain
magnetically trapped after the MOT beams have been turned off.
Assuming these atoms can be described by an \emph{effective}
magnetic moment $\bar{\mu}\geq 0$, then the potential energy for a
trapped atom is given by
\begin{equation*}\label{Eq:MagneticEnergy}U(x,y,z) = \bar{\mu} \sqrt{dB_x^2 x^2 +
dB_y^2 y^2+dB_z^2z^2}+M g z
\end{equation*} where $M$ is the atomic mass, and $g$ is the
gravitational acceleration. In equilibrium at temperature $T$, the
marginal distributions of atoms along the $x$ and $z$ directions are
found to be
\begin{equation*}\begin{split} &p_x(x) =N_x\exp\left(-
2\frac{|x|}{\bar{x}}\sqrt{1-\bar{g}^2}\right)\left(1+2\frac{|x|}{\bar{x}}\sqrt{1-\bar{g}^2}\right)\\
&p_z(z) = N_z \exp\left(- 2\frac{|z|}{\bar{z}}-
2\bar{g}\frac{z}{\bar{z}}\right)\left(1+2\frac{|z|}{\bar{z}}\right)
\end{split}\end{equation*} where $N_x$ and $N_z$ are normalization constants.
The length scales and effects of gravity are characterized by
$\bar{g}~=~Mg/(\bar{\mu}dB_z)$ and $\bar{x}=2k_B T/(\bar{\mu}
dB_{x})$ with a similar expression for $\bar{z}$. Such a trap is
stable against gravity only if $\bar{g}<1$. For $\bar{g}=0$,
$\bar{x}$ and $\bar{z}$ are the standard deviations of the cloud,
while for $\bar{g}>0$, the cloud stretches along gravity. We capture
images of magnetically trapped atomic clouds by loading the MOT for
2~s then extinguishing the loading and trapping lasers, leaving all
magnetic fields unchanged. After a variable hold time $\Delta t$, we
pulse the MOT beams and capture the resulting fluorescence image.
Such an image is shown in Fig.~\ref{Fig:BTrapImage}. By analyzing
these images, we can determine $T$, $\bar{\mu}$, and atom number as
the magnetic trap evolves in time \footnote{In our image analysis,
we do not account for Zeeman shifts over the spatial extent of the
cloud, which are expected to alter the excited state fraction by at
most 7\%.}.

\begin{figure}[t]
\includegraphics*[width=3.2in]{BTrapDecay_Origin}
\caption{\label{Fig:BTrapDecay}(Color online) Evolution of the
magnetically trapped atom number. After a time $\Delta t\approx
0.2$~s, all of the metastable atoms have returned to the ground
state and the trap decays exponentially with a time constant of
$(270\pm13)$~ms (solid line). (Inset) The temperature of the trapped
atoms increases over time with an error-weighted heating rate of
$(14\pm4)~\mu$K/s (solid line).}
\end{figure}
If the magnetic trap loading efficiency is independent of an atom's
energy, then thermodynamic arguments
\cite{Stuhler:2001a:PhysRevA.64.031405} predict that the
magnetically trapped atom temperature $T_B$ should be related to the
MOT temperature $T_M$ by $T_B\geq T_M/3$ (equality holds when the
spatial extent of the MOT is much smaller than that of the magnetic
trap). However, our magnetic trap temperatures systematically
violate this bound as shown in Fig.~\ref{Fig:MOTTvsBTrapT} where we
measure temperatures below 25~$\mu$K. Similar anomalously low
temperatures have been observed in MOT-loaded, magnetically trapped
Ca \cite{Hansen:2003a} and Sr \cite{Nagel:2003a} while unexpectedly
high temperatures were reported in Cr
\cite{Stuhler:2001a:PhysRevA.64.031405}. In Ref. \cite{Nagel:2003a},
a selection effect during trap loading was suggested to explain the
low temperatures. We suspect a similar phenomenon is at work here,
in which spatially-dependent optical pumping in the MOT tends to
correlate an atom's spin state with its energy; the ejection of
atoms in untrapped magnetic states ($m_J\leq 0$) may then reduce the
average energy of the sample when the MOT is turned off. This
hypothesis cannot be quantitatively evaluated without studying the
loading efficiency of the magnetic trap, which is difficult to
assess due to the presence of trapped, dark-state atoms that cannot
be probed directly.

An example of the time evolution of the magnetic trap population and
temperature is shown in Fig.~\ref{Fig:BTrapDecay}. Before $\Delta
t\approx 0.2~$s, the trap population appears to increase. This is
due to the relaxation of metastable atoms back to the ground
electronic state where they are resonant with the probe light (MOT
beams). After this time, we see single-exponential decay, most
likely arising from collisions with the background gas, with a
pressure around $2.3\times10^{-6}$~Pa \footnote{We find good fits to
the magnetic trap relaxation data of Fig.~\ref{Fig:BTrapDecay} for a
simple model in which a fraction $f$ of initially dark (metastable)
atoms relaxes to the ground state exponentially while atoms are
simultaneously lost from the trap at a second exponential rate.
However, this larger number of fit parameters cannot give an
independent indication of both the initial dark state fraction and
the metastable relaxation rate so we have only included a fit to the
long-time relaxation behavior here.}. This claim is further
supported by a strong correlation between background gas pressure
and trap lifetime and the similar lifetimes for $^{166}$Er and
$^{168}$Er (data not shown). We do not suspect that the trap decay
rate is due to density-dependent collisional effects within the
trap, because the number densities in the present measurements are
not particularly large [we observe peak densities of
$n=(2.1\pm0.3)\times 10^{9}$~cm$^{-3}$] and because we see no
evidence of non-exponential decay. On the other hand, we do not know
the origin of the observed heating rate, which may arise from
dipolar relaxation processes or collisional heating mechanisms.
Improvements to the vacuum system should eliminate our dominant loss
channel and allow us to study other relaxation processes. Finally,
our data indicate a phase-space density $\rho\approx
2.5\times10^{-8}$. It should be straightforward to increase this by
at least an order of magnitude by increasing the atomic beam flux.

\section{Conclusions}
In summary, we have presented a number of favorable cooling and
trapping properties for atomic erbium. We showed that efficient
sub-Doppler cooling occurs in our MOT without the need for an
additional optical molasses cooling stage. Furthermore, sub-Doppler
cooling should remain effective in large magnetic fields due to the
near degeneracy of Land\'{e} $g$ factors. This phenomenon may have
interesting applications for laser cooling in large magnetic
gradient traps, and might be useful for constructing atomic
waveguides along magnetic field lines. We also studied the
thermodynamics of a continuously-loaded magnetic trap, where
temperatures reach unexpectedly low values. The efficient cooling
and trapping methods discussed here should be straightforwardly
useful as a starting point for phase-space compression towards
quantum degeneracy. However, magnetic and collisional relaxation
processes must be studied before such a procedure can be considered.

\acknowledgements{The authors gratefully acknowledge J.~L. Hanssen,
N. Lundblad, and J.~V. Porto for stimulating discussions and D.
Rutter and A. Band for technical assistance. S.~L. acknowledges
support from NIST and NSF. A.~B. received financial support from the
National Research Council.}