\documentclass[twoside]{article} 
\usepackage{amsmath,amsfonts,graphicx} 
\pagestyle{myheadings}

\markboth{\hfil Monotone solution for a sedimenting sphere
\hfil EJDE--2001/62}
{EJDE--2001/62\hfil A. Belmonte, J. Jacobsen \& A. Jayaraman \hfil}

\begin{document} 

\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
{\sc  Electronic Journal of Differential Equations},
Vol. {\bf 2001}(2001), No. 62, pp. 1--17. \newline
ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp  ejde.math.swt.edu  (login: ftp)}
\vspace{\bigskipamount} \\
%
 Monotone solutions of a nonautonomous differential equation
 for a sedimenting sphere
%
\thanks{ {\em Mathematics Subject Classifications:} 34C60, 34D05, 76D03. 
\hfil\break\indent
{\em Key words:} sedimenting sphere, unsteady Stokes flow, 
\hfil\break\indent nonautonomous ordinary differential equations,
monotone solutions.
\hfil\break\indent
\copyright 2001 Southwest Texas State University. \hfil\break\indent
Submitted January 10, 2001. Published September 24, 2001.} }

\date{}

\author{Andrew Belmonte, Jon Jacobsen, \& Anandhan Jayaraman}

\maketitle 
 
\begin{abstract} 
  We study a class of integrodifferential  equations
 and related 
  ordinary differential equations for the  initial value 
problem
  of a rigid sphere falling through an infinite  fluid medium.  
  We prove that   for creeping \textit{Newtonian} flow,
  the motion of the sphere is monotone in its  approach to the
  steady state solution given by the Stokes drag. We discuss this
  property in
terms of a general nonautonomous second order 
  differential equation,
focusing on a decaying nonautonomous term
  motivated by the sedimenting sphere problem. 
\end{abstract} 

\newtheorem{definition}{Definition}[section]
\newtheorem{theorem}[definition]{Theorem}
\newtheorem{corollary}[definition]{Corollary}

\renewcommand{\theequation}{\arabic{section}.\arabic{equation}}
\catcode`@=11
\@addtoreset{equation}{section}
\catcode`@=12

\section{Introduction}
 
A rigid sphere falling through a viscous medium is a classic problem
in fluid dynamics, which was first solved in the steady state for 
the limit of vanishingly small Reynolds number in an infinite  domain
by G. G. Stokes in 1851 \cite{stokes:51}.   The time-dependent
approach to the steady state allows the partial differential equation
for the  sphere and fluid to be reduced to an integrodifferential
equation  for the motion of the sphere in  an infinite medium. The
physical effects included in this equation are the  buoyancy which
drives the motion, the inertia of the sphere, the  viscous drag, an
added mass term, and a memory or history  integral \cite{clift}. 



In the case of a {\it Newtonian fluid}, the main effect of the memory 
integral on the dynamics is to modify the approach to steady state 
from   exponential to algebraic.  The integral also makes the 
equation effectively second order, though it is generally accepted 
that no oscillations occur as the sphere reaches its steady state 
value   \cite{clift,walt92}. Physically it is clear that no
oscillations can occur due to the absence  of a restoring force
against gravity, and a sphere   released from  rest in a Newtonian
fluid at low Reynolds number is observed to reach its terminal 
velocity monotonically. However, it is not directly evident
mathematically that oscillating solutions are precluded,
particularly as the governing integrodifferential equation can be
transformed to a  nonautonomous second order ordinary differential
equation which has  the form of a harmonic oscillator
\cite{clift,villat:44}. 


In a {\it non-Newtonian fluid}, such as a polymer solution,  a falling
sphere is often observed to  undergo transient  oscillations before
reaching its terminal velocity  \cite{arigo97,walt92}. These 
oscillations occur due to the elasticity of the fluid, which provides
a restoring force \cite{bird87}. The steady state value is of primary
concern  in many applications, and much work focuses only on this
aspect of  the problem. The oscillations which occur during the
approach to  steady state have been reproduced in a linear
viscoelastic model by King and Waters \cite{king:ums72}. 
 
More recently, \textit{nontransient} oscillations of falling spheres
(and  rising bubbles)
have been observed in specific aqueous solutions
of surfactants (wormlike micellar solutions) 
\cite{belmonte:soc00,jayaraman:osf00}. These observations were
initially made for a bubble in the wormlike micellar fluid CTAB/NaSal
\cite{hu98,rehage82}, which showed oscillations in its position and
shape.  The shape
oscillations included an apparent cusp which momentarily appears at
the trailing end of the bubble. Such a cusplike tail is a well known
property of rising bubbles in non-Newtonian fluids
\cite{hass79,liu95}, which we initially believed to play an important
role in the micellar oscillations. 
 
\begin{figure}                                                                                                                    
\begin{center} 
\includegraphics[width=0.48\textwidth]{fig1a.eps}
\quad
\includegraphics[width=0.48\textwidth]{fig1b.eps}
\end{center} 
\caption{Motion of a 1/8'' diameter teflon sphere falling through an  
aqueous solution of 6 mM CTAB/NaSal: (left) position vs time; (right)  
calculated velocity vs time \cite{jayaraman:osf00}.} 
\label{f-expt}                                                                                                                   
\end{figure}                                                                     
 
 
Subsequent observations of solid spheres which also oscillate while  
falling through the same solutions made it   clear that the cusp is
not involved in the phenomenon, and that   another explanation must
be sought. Unlike a sedimenting sphere in a conventional
non-Newtonian fluid, these oscillations do not appear to be transient
\cite{jayaraman:osf00}.  An example is given
in Figure \ref{f-expt}, which shows the motion of a 1/8''   teflon
sphere falling through a tube ($L = 98$ cm, $R = 3.2$ cm)   filled
with a 6 mM 1:1 solution of CTAB/NaSal \cite{jayaraman:osf00}. 

 
Our attempts to model this phenomenon  brought to our attention the
unusual aspects of the  integrodifferential equation for a falling
sphere. We prove here that the equation for sedimenting sphere in a
Newtonian fluid in the limit of zero Reynolds number (creeping flow)
does not  admit   oscillating solutions, despite some appearances
that it  does.  This result is due to the special  properties of the
error   function when multiplied by oscillating  functions.  It is
ultimately   related to the stability of  nonautonomous ordinary
differential equations with monotone secular terms,  which is
appropriately  viewed as an initial value problem, and not in terms
of linear stability analysis around the terminal velocity. 

 
\section{The Motion of a Sedimenting Sphere }
 
We begin by reviewing some classical results for
the equation of motion governing a falling sphere in
a viscous Newtonian fluid of infinite extent (no boundaries).

\subsection{Equation of Motion of the Sphere}
An incompressible fluid in the absence of body forces is described 
by the equations
\begin{gather}
\label{ce2}
\rho \left( \frac{\partial \vec u}{\partial t} + (\vec u \cdot 
\nabla) \vec u \right)  =  -\nabla p +  \, \mathop{\rm div} \sigma, \\
                \mathop{\rm div} \vec u  =  0,
\label{ce3}
\end{gather}
where $\rho(\vec x,t)$ is
the density of the fluid, $p$ is the  pressure,
$\vec u(\vec x,t)$ is the velocity field for the fluid, and
$\sigma(\vec x,t)$ is the extra stress tensor, which measures 
force per unit area (other than pressure) in the
present configuration of the fluid.  
A Newtonian fluid is a fluid for which the stress tensor $\sigma$ is
linearly related to the rate of strain tensor $D$ through the relation
\begin{equation}
   \sigma = 2 \mu D,  
\label{st}
\end{equation}
where  $\mu$ is the viscosity of the fluid 
and $D = (\nabla \vec u + (\nabla \vec u)^T)/2$ is the symmetric part of
the velocity gradient $\nabla \vec u$.   
 From (\ref{ce2}),(\ref{ce3}), and (\ref{st}) one
obtains the Navier-Stokes equation:
\begin{equation}
\rho \left(\frac{\partial \vec u}{\partial t} + (\vec u \cdot 
\nabla) \vec u \right)  =  -\nabla p + \mu \Delta \vec u.
\label{NS}
\end{equation}

Non-Newtonian fluids are fluids for which the assumption (\ref{st}) is
invalid.  For instance, polymeric and viscoelastic fluids often fail to
conform to the instantaneous relation between stress and velocity gradients
implicit in (\ref{st}). In general $\sigma=\sigma(t,\vec x,D)$ will depend nonlinearly
on $D$ and on the past history of stress in the fluid.

By choosing an appropriate time and length scale, (\ref{NS}) can 
be written in a nondimensional form
\begin{gather}
 \frac{\partial \tilde u}{\partial \tau} +\mbox{Re}\; (\tilde u \cdot 
\nabla) \tilde u  =  -\nabla \tilde p +  \Delta \tilde u, \label{NS2} \\
                \mathop{\rm div} \tilde u =  0, \label{incomp2}
\end{gather}
where, $\tilde u$ and $\tilde p$ are the nondimensional velocity and pressure. 
The dimensionless constant $\mbox{Re}$ is called the Reynolds number and it 
measures the relative importance of inertial effects to that of viscous effects.
When the inertial effects are negligible ($\mbox{Re}=0$), equation
(\ref{NS}) is called
the Stokes equation. In this paper we restrict our analysis to this situation.

Stokes solved the steady version of  (\ref{NS2})-(\ref{incomp2}) for the case 
of sphere falling though the fluid for vanishing Reynolds number. The Stokes 
solution gives the steady state drag on the sphere of radius $R$ falling 
through a fluid with a steady speed $U$ to be $F=6 \pi \mu R U$.
However, in order to solve the transient problem of falling sphere, we first 
solve the problem of sphere oscillating with a frequency $\omega$ and compute 
the drag experienced by the sphere as a function of $\omega$.  The drag experienced 
by a sphere falling at a arbitrary speed $U(t)$ can then be computed as a Fourier 
integral of this drag.

Consider a sphere of radius $R$ and density $\rho_s$ in a Newtonian
fluid of density $\rho$ and viscosity $\mu$.  The exterior Stokes
flow driven by small oscillations of the sphere at a frequency
$\omega$ can be solved exactly \cite{basset,lamb}, which leads to a
hydrodynamic force dependent on both $U$ and $dU/dt$ : 
\begin{eqnarray} 
F=6 \pi \mu R\left(1+\frac{R}{\delta}\right) U + 3 \pi R^2 \rho \delta 
\left(1+\frac{2R}{9\delta}\right) \frac{dU}{dt},\label{eq:Fw-drag} 
\end{eqnarray} 
where $\delta=\sqrt{2\nu/\omega}$ is a diffusive lengthscale common  
to Stokes problems, and $\nu = \mu/\rho$ is the kinematic viscosity.
  Using this, the
general time-dependent problem of the motion of a falling sphere can
be reduced from a partial to an ordinary differential equation for
the speed $U(t)$ of the sphere, an exact equation which takes into
account the motion of the surrounding fluid \cite{clift,villat:44}.  
For a sphere moving with an arbitrary speed $U(t)$, the 
hydrodynamic drag it experiences can be calculated by representing 
$U(t)$ as a Fourier integral: 
$$ 
U(t)= \int^\infty_{-\infty} U_\omega e^{-i \omega t}\; d\omega.  
$$ 
The drag for each Fourier component is then given by  
(\ref{eq:Fw-drag}).  The total hydrodynamic drag on the sphere is  
obtained by integrating over all Fourier components, leading to 
\begin{eqnarray} 
F_{drag} = 6\pi\mu R U(t) + \frac{1}{2} \rho \mathcal{V} \frac{dU}{dt}  
+ 
 6\pi\rho R^2 \sqrt{\frac{\nu}{\pi}} \int^t_{-\infty}  
\frac{U'(s)}{\sqrt{t-s}}\;ds \label{eq:F-net} 
\end{eqnarray} 
where
the first term represents the  
steady state drag on a sphere falling with a velocity $U$,  
the second term represents the added mass term (the force  
required to accelerate the surrounding fluid), the  
third term is the Basset memory term, and
 $\mathcal{V}=4 \pi R^3/3$ is the volume of the sphere. 
If the sphere starts from rest, then the lower limit of the integral  
in (\ref{eq:F-net}) starts from $\tau=0$ instead of  
$\tau=-\infty$.  The expression for the unsteady drag force can then  
be substituted into the balance of force equation for the sphere: 
$$ 
\rho_s \mathcal{V} U'(t) = F_{buoy} - F_{drag}. 
$$ 
Thus the equation of motion for the sphere is 
\begin{eqnarray} 
\lefteqn{(\rho_s+\frac{\rho}{2}) \mathcal{V} U'(t) \; 
+ \; 6\pi \mu R U(t) \,+ \;6\pi R^{2}\sqrt{\frac{\rho\mu}{\pi}} \int^t_0 
\frac{U'(s)}{\sqrt{t-s}}\, ds  }\nonumber\\ 
&  = & (\rho_s-\rho) \mathcal{V} g,  \hspace{6cm}
\label{eq-forces} 
\end{eqnarray} 
which can be rewritten in the simpler form 
\begin{equation} 
U'(t) + B U(t) + Q \int^t_0 \frac{U'(s)}{\sqrt{t-s}}\, ds  
 =M, \label{eq:uprime} 
\end{equation} 
where 
\begin{gather} 
B=\frac{9\mu}{R^2\left(2\rho_s + \rho\right)}, \label{eq:B} \\ 
Q=\frac{9\rho}{R\left(2\rho_s+\rho\right)}  
\sqrt{\frac{\mu}{\rho\pi}}, \label{eq:Q} \\ 
M=\frac{2 g \Delta \rho}{2\rho_s + \rho}, \label{eq:M} 
\end{gather} 
and $\Delta\rho = \rho_s - \rho$ is the density difference which 
drives  the motion.  In this approach the motion of the sphere is 
described by an integrodifferential equation whose
integral term has the same singularity as Abel's equation 
\cite{keener:pam00,estrada:sie00}. 
Note that this equation is only valid 
in the limit of zero Reynolds number \cite{ock68,michael97}. 
  
Physically one expects the solution $U(t)$ to approach  
a terminal velocity.  It is clear from 
(\ref{eq:uprime}) that the only steady state solution ($U'=0$)  
possible is  
\begin{equation} 
U_0 = \frac{M}{B} = \frac{2\Delta\rho gR^2}{9\mu},  
\label{e-stokes} 
\end{equation} 
which is the classical result of balancing the Stokes drag with the  
buoyancy. 
 
\subsection{Solving the Integrodifferential Equation} \label{solnIDE} 
 
We first rewrite the  integrodifferential equation (IDE) for the 
sphere in a  nondimensional form using $U_0$ as the velocity scale, 
and $1/B$, the viscous diffusion time, as  the time scale. The 
variables are 
$$ 
\tau = Bt \qquad \text{and} \qquad u(\tau) = U(\tau/B )/U_0. 
$$ 
With this rescaling (\ref{eq:uprime}) becomes
\begin{equation} 
u'(\tau) + u(\tau) + \sqrt{\frac{\kappa}{\pi}}\int^{\tau}_0  
\frac{u'(s)}{\sqrt{\tau-s}}\, ds =1, \label{eq:maineq} 
\end{equation} 
where the control parameter $\kappa$ is given by 
\begin{equation} 
\kappa = \frac{\pi Q^2}{B}=\frac{9\rho}{2\rho_s+\rho}. 
\label{eq:kappa} 
\end{equation}
Thus the motion of the sphere depends only on the relative densities 
of  the sphere and the fluid through the parameter $\kappa$. 
The  density of the sphere $\rho_s$ can range from zero to infinity,
which implies a parameter range  $0 < \kappa < 9$. When the density
of the sphere is  equal to the density of the fluid,  $\kappa=3$ and
$U_{0}=0$. Here we will only be concerned with falling spheres, for
which $0 < \kappa < 3$. 

 
Although we will solve the IDE (\ref{eq:maineq}) directly, it is of  
interest to connect the problem to ordinary differential equations  
(ODE's) and discuss some important consequences therein, especially  
with regard to the stability of the terminal velocity solution.  
Following Villat \cite{villat:44}, we can  
rewrite (\ref{eq:maineq}) as an ODE using Abel's Theorem 
(see Appendix):
\begin{gather} 
u'' + (2-\kappa) u' + u = 1- \sqrt{\frac{\kappa}{\pi  
\tau}},  
\label{eq:non-dim1} \\ 
u(0) = 0, \quad u'(0) = 1.   
\label{eq:non-dim2} 
\end{gather} 
More general initial conditions $U(0) \neq 0$ and $U'(0)=M -  
BU(0)$ lead to a slightly different ODE: 
\begin{gather} 
u'' + (2-\kappa) u' + u  =  1 + \sqrt{\frac{\kappa}{\pi  
\tau}}(u(0) - 1), 
\label{eq:gen1} \\ 
u(0)  =  \xi, \quad u'(0)  =  1 - \xi. \label{eq:gen2}  
\end{gather} 
Note that it follows from equation (\ref{eq:maineq}) that $u'(0)$ is
prescribed in terms of $u(0)$, and thus the second order ODE we have
obtained requires only one initial condition.


Since we are investigating the possibility of steady-state
oscillations of a sedimenting sphere, we are primarily concerned with
the asymptotic behavior of (\ref{eq:non-dim1}). Moreover, since the
nonautonomous term tends to zero as $t\rightarrow\infty$, one might
expect the  stability of (\ref{eq:non-dim1}) to mimic the homogeneous
problem. With this in mind, let $\alpha$ and $\beta$ denote the roots
of the characteristic equation 
\begin{equation}
m^2+(2-\kappa)m+1 =0. 
\label{ce}
\end{equation}
It is readily verified that (\ref{ce}) has complex roots
for  $0 < \kappa < 4$.  Moreover, the roots have positive real parts
for   $2 < \kappa < 4$.  
Since the relevant range of $\kappa$ for a falling
sphere is $0 < \kappa < 3$, one sees that oscillations are
not a priori precluded.
In terms of the actual densities of the fluid and  sphere
the condition for complex roots corresponds to  $\rho_s > (5/8) \rho$, 
which is true in the case of a heavy sphere falling through a lighter 
liquid ($\rho_s > \rho$).  If additionally
$\rho_s < (7/4) \rho$, then the  complex roots have positive real parts.
If we rewrite (\ref{eq:non-dim1}) in the asymptotic limit ($t \rightarrow \infty$)
as a  first  order linear system
\begin{eqnarray} 
x' &=& y \\ 
y' &=& (1-x) + (\kappa-2) y, 
\end{eqnarray}
where $x=u$, and $y = u'$, then 
$(x,y)=(1,0)$ is the unique equilibrium point,
which corresponds to the terminal velocity.  The  
eigenvalues of this  system are precisely $\alpha$ and   $\beta$, whence
the equilibrium  point becomes unstable.

Nonetheless, as we will show, even in this range ($2 < \kappa < 3$), the solution 
to the full equation (\ref{eq:non-dim1}) \textit{monotonically} approaches
the value 1, corresponding to the monotonic 
approach to the steady Stokes value (\ref{e-stokes}) for the actual velocity $U$.    
Clearly the nonautonomous term continues to play a dominant role in the stability
of (\ref{eq:non-dim1}), despite its algebraic approach to zero.

To solve for $u=u(\tau)$, we return to the IDE (\ref{eq:maineq}) and apply the Laplace  
transform in the case $u(0) = 0$:  
\[ 
\mathcal{L}\{u\}(s) = \frac{1}{s(s + \sqrt{\kappa} \sqrt{s} + 1)} \quad   
\text{or} \quad 
\mathcal{L}\{u'\}(s) = s \mathcal{L}\{u\} = \frac{1}{s + \sqrt{\kappa} \sqrt{s} + 1} . 
\] 
Since 
 $\alpha \beta = 1$ and $ \sqrt{\alpha}+ \sqrt{\beta} = \sqrt{\kappa} $, we may 
express this last equation in the form 
\[   
\mathcal{L}\{u'\}(s)  = \frac{\sqrt{\kappa}}{\alpha - \beta}  
\left[\frac{\sqrt{\alpha}}{\sqrt{s}(\sqrt{s} + \sqrt{\alpha})} 
- \frac{\sqrt{\beta}}{\sqrt{s}(\sqrt{s} + \sqrt{\beta})}\right]. 
\] 
Moreover, the identity  
\[   
 \mathcal{L}\{e^{\alpha t} \text{Erfc} \sqrt{\alpha t} \} = \frac{1}{\sqrt{s}(\sqrt{s} +  
\sqrt{\alpha})},   
\] 
implies 
\begin{equation} 
   u'(\tau) = \frac{\sqrt{\kappa}}{\alpha - \beta} 
                \left[\sqrt{\alpha} e^{\alpha \tau} \text{Erfc} \sqrt{\alpha \tau} -  
                \sqrt{\beta} e^{\beta \tau} \text{Erfc} \sqrt{\beta \tau} \right].   
  \label{eq:uprime2} 
\end{equation} 
Finally, since  
\[  \sqrt{\alpha} \int e^{\alpha t} \text{Erfc} \sqrt{\alpha t} \, dt = 2 \sqrt{\frac{t}{  
\pi}} +  
\frac{1}{\sqrt{\alpha}} e^{\alpha t} \text{Erfc} \sqrt{\alpha t} + C,  
\] 
we find the solution  
\begin{eqnarray}  
u(\tau) = 1 + \frac{\sqrt{\kappa}}{\alpha -\beta} 
\left[ \frac{e^{\alpha \tau} \text{Erfc} \sqrt{\alpha \tau}}{\sqrt{\alpha}}  -   
\frac{e^{\beta \tau}\text{Erfc} \sqrt{\beta \tau}}{\sqrt{\beta}}  
\right]. 
\label{eq:soln} 
\end{eqnarray} 
 
This solution to the IDE (\ref{eq:maineq}) is 
also the solution to the ODE (\ref{eq:non-dim1})-(\ref{eq:non-dim2}). 
Applying transform methods to the more general set of equations defined by 
(\ref{eq:gen1})-(\ref{eq:gen2}) one finds the solution
\begin{equation}    
u(\tau) = (1-\epsilon) u_0(\tau) 
+ \epsilon\, \qquad \epsilon \equiv u(0),\, u'(0)=1-\epsilon,
\label{220}
\end{equation}
where $u_0(\tau)$ is the solution defined by (\ref{eq:soln}).  Note that
the solution for arbitrary initial velocity $\epsilon$ is a simple
rescaling of the solution for the sphere  initially at rest.  It is not
obvious from the form of $u$ in (\ref{eq:soln}) that the values approach 1 monotonically.
Let us first investigate the asymptotic behavior of this solution.

\subsection{Asymptotic Approach to the Steady Stokes Solution} 
\label{asympt} 
 
Finding the asymptotic behavior of the solution $u_0$ is straightforward.
We employ the asymptotic expansion of the error function
\cite{gautschi:eff92} 
\begin{eqnarray} 
e^{z^2} \text{Erfc} \, z  \sim \frac{1}{\sqrt{\pi}z}\left[1 + 
  \sum_{m=1}^\infty (-1)^m \frac{1\cdot3 \cdots  
(2m-1)}{(2z^2)^m}\right] 
\end{eqnarray} 
as $z\rightarrow\infty$, provided $|\arg (z)| < 3\pi/4$. 
Since $0 < \arg(\sqrt{\alpha t}) < \pi/2$ for $0 < \kappa < 3$, we  
see that asymptotically  
$$ 
\text{Erfc} \sqrt{\alpha \tau} \sim \frac{e^{-\alpha \tau}}{\sqrt{\alpha  
\tau}}. 
$$ 
 
Thus the product of the exponential term and the error function  
approaches zero in the limit $t\rightarrow\infty$. Using this  
expansion in (\ref{eq:soln}) we obtain 
\begin{eqnarray} 
 \lim_{\tau \rightarrow \infty} u_0(\tau) = 1 \quad \text{ or }  
 \quad \lim_{t \rightarrow \infty} U(t) = U_0. 
\end{eqnarray} 
As the solution for any initial condition is a 
rescaling of $u_0$, we see this limit applies for all values of $u(0)$. 
 
Although the asymptotics of this solution are clear, the transient
solution has some unusual properties. Numerical
simulation of the IDE (\ref{eq:maineq}), or even attempts to plot the
analytic solution (\ref{eq:soln}), eventually  blow up at
large $\tau$ when $\kappa$ is in the unstable range. 
Clearly the cancellation between the exponentially growing and 
decreasing terms are quite sensitive to numerical errors.  This is an  
indication that the product $e^{\alpha \tau}\text{Erfc} \sqrt{\alpha
\tau}$ should be considered as a special function with its own
properties.

 
\section{Monotonicity of the Transient Solution} 
\label{monosection}
 
We begin with the transient solution to the IDE or ODE considered
in   the previous section.  In addition to the insensitivity of the  
transient solution to the real part of the homogeneous roots, it is  
surprising that the nonzero complex part  does not lead to
{\it any} oscillations in the velocity of the sphere, although there
has sometimes been confusion on this point regarding transient
oscillations   \cite{villat:44}.  Experimentally the sphere in a
Newtonian fluid has never   been observed to oscillate, in
contradistinction to most non-Newtonian   (particularly elastic)
fluids \cite{arigo97,walt92}. We will show   that the solution $u$
defined by (\ref{eq:soln}) is monotone as a   function of $\tau$ for all
$\kappa \in (0,4)$.  Although this may be well   known, we have not yet
found a reference to any proof other than the case  
$\alpha,\beta \in \mathbb{R}$ \cite{wil78}, which 
corresponds to $\kappa > 4$. 

 
Let us define the function $\text{Vi}:\mathbb{C} \rightarrow \mathbb{C}$ by
\begin{equation} 
\text{Vi}(z) \doteq e^{z}\text{Erfc} \sqrt{z}  
= \frac{2e^{z}}{\sqrt{\pi}} \int_{\sqrt{z}}^{\infty} e^{-s^{2}}\,ds. 
\label{defvi}
\end{equation}
We shall refer to $\text{Vi}$ as the Villat function, since this combination 
appeared in the explicit solution of the differential equation for the falling 
sphere problem by Villat \cite{villat:44}.  Closely related to $\text{Vi}(z)$ is
the ``plasma  dispersion  function'' $w$, defined by
\cite{gautschi:eff92}:  \[   w(z) = e^{-z^2} \text{Erfc} (-i z).  \]  In
fact, $\text{Vi}(\alpha t) = w(i \sqrt{\alpha t})$.   Using the Villat function
we may now prove the main theorem.

\begin{theorem}[Monotonicity] For each $\kappa \in (0,4)$, the function
\begin{equation}  
u(t) = 1 + \frac{\sqrt{\kappa}}{\alpha -\beta} 
\left[ \frac{e^{\alpha t} \text{Erfc} \sqrt{\alpha t}}{\sqrt{\alpha}}  -   
\frac{e^{\beta t}\text{Erfc} \sqrt{\beta t}}{\sqrt{\beta}}  
\right] 
\end{equation}
approaches the limit 1 monotonically.
\end{theorem} 

\paragraph{Proof.}
We have shown that $u(t) \rightarrow 1$, thus it remains to show it does so  
monotonically.  We will demonstrate this by proving $u'(t) > 0$ for  
all $t > 0$.  To this end, fix $t > 0$, $\kappa \in (0,4)$ and recall that $\alpha$ and  
$\beta$ denote the conjugate pair of roots of the polynomial $m^2 +  
(2-\kappa)m + 1$.  Recall from (\ref{eq:uprime2}) that 
\[   
u'(t) = \frac{\sqrt{\kappa}}{\alpha - \beta} 
        \left[ \sqrt{\alpha} \, \text{Vi}(\alpha t) - \sqrt{\beta} \, \text{Vi}(\beta t)  
\right].    
\] 
Since $\text{Erfc}(\overline{z}) = \overline{\text{Erfc} \, z}$, it follows that  
$\text{Vi}(\beta t) = \overline{\text{Vi} (\alpha t)}$ and 
\begin{equation} 
u'(t) = \sqrt{\kappa} \,  \frac{\sqrt{\alpha} \, \text{Vi}(\alpha t) -  
\overline{\sqrt{\alpha} \, \text{Vi}(\alpha t)}}{\alpha - \overline{\alpha}} 
= \sqrt{\kappa} \,  \frac{\Im \{\sqrt{\alpha} \, \text{Vi}(\alpha t) \}}{\Im \{\alpha \}}.       
\label{eq:uprimeform2} 
\end{equation} 
Since $\Im \{ \alpha \} > 0$ for each $\kappa \in (0, 4)$, it is evident from (\ref{eq:uprimeform2})  
that the sign of $u'(t)$ is determined by the imaginary part of  
the function $\sqrt{\alpha} \, \text{Vi}(\alpha t)$.   For the  plasma  
dispersion function $w$ introduced above, the real and imaginary parts are given by 
(see e.g., \cite[7.4.13-7.4.14]{gautschi:eff92}) 
\[  
\Re(w(x + i y)) = \frac{1}{\pi} \int_{-\infty}^{\infty} \frac{y  
e^{-s^2}}{(x-s)^2 + y^2} \, ds \quad (x \in \mathbb{R}, y > 0) 
\] 
and  
\[  
\Im(w(x + i y)) = \frac{1}{\pi} \int_{-\infty}^{\infty} \frac{(x-s)  
e^{-s^2}}{(x-s)^2 + y^2} \, ds \quad (x \in \mathbb{R}, y > 0). 
\] 
 It is readily verified that  
$| \alpha | = 1$, thus in polar form we have $\alpha = e^{i \theta}$ for  
some fixed $\theta \in (0,\pi)$, in which case $\text{Vi}(\alpha  
t) = w(i \sqrt{\alpha t}) = w(x + i y)$, where 
$$ 
x = -\sqrt{t} \sin \left( \frac{\theta}{2} \right) 
\qquad \text{and} \qquad y = \sqrt{t}  
\cos \left( \frac{\theta}{2}\right) . 
$$   
Using this information we compute 
\begin{eqnarray} 
 \Im \{\sqrt{\alpha} \, \text{Vi}(\alpha t)\} & = & 
\cos \left( \frac{\theta}{2}\right) \, \Im  
\{ w(x + i y) \} 
                                        + 
\sin \left( \frac{\theta}{2} \right) \,  
\Re \{ w(x + i y)\}  
\nonumber  \\ 
                & = & \frac{1}{\pi} \cos \left( \frac{\theta}{2} \right)  
\int_{-\infty}^{\infty} \frac{(x-s) e^{-s^2}}{(x-s)^2 + 
y^2} \, ds \nonumber \\
& & 
 + \frac{1}{\pi} \sin \left( \frac{\theta}{2} \right) 
\int_{-\infty}^{\infty}  
\frac{y e^{-s^2}}{(x-s)^2 + y^2} \, ds  
 \\ 
         & = & \frac{1}{\pi} \cos \left(  \frac{\theta}{2} \right) 
\int_{-\infty}^{\infty} \frac{-s e^{-s^2}}{(x-s)^2 + y^2} \, ds. 
\label{eq:poseq} 
\end{eqnarray} 
In the last step  we have used the fact that 
$x \cos(\theta/2) + y  \sin(\theta/2) = 0$. 
Since $\sqrt{\alpha}$ lies in the first quadrant, the prefactor of the 
last integral above is positive and we may conclude 
that  $\Im \{\sqrt{\alpha} \, \text{Vi}(\alpha t)\} > 0$ 
provided 
\begin{equation} 
\int_{-\infty}^{\infty} \frac{s e^{-s^2}}{(x-s)^2 + y^2} \, ds  
=  \int_{-\infty}^{\infty} \frac{s e^{-s^2}}{(s + \sqrt{t}  
\sin \frac{\theta}{2})^2 + t \cos^2 \frac{\theta}{2} } \, ds < 0.   
\label{negeq} 
\end{equation} 
Let us denote the integrand as  
\[ F(s) = \frac{s e^{-s^2}}{P(s)} \quad \text{where} \quad 
 P(s) = \left(s + \sqrt{t} \sin \frac{\theta}{2} \right)^2 + t  
\cos^2 \frac{\theta}{2}.  \] 
Note that $P(s) > 0$ for $s \in  
\mathbb{R}$ (recall $t > 0 $ is fixed).  The proof is complete once the 
following two  observations are made: 

\begin{enumerate} 
\item[($a$)] \, $| F(-s) | > F(s)$ for $s > 0$;  
\item[($b$)] \, $\int_\mathbb{R}  F \, ds = \int_0^{\infty} F(s) \, ds - \int_0^{\infty} |  
F(-s) | \, ds$.
\end{enumerate} 
To see $(a)$, notice for $s > 0$ we have $0 < P(-s) < P(s)$, thus 
$$ 
|F(-s)| = \frac{s e^{-s^2}}{P(-s)} > \frac{s e^{-s^2}}{P(s)} =  
F(s),  
\qquad \text{for }  s > 0. 
$$ 
Observation $(b)$ follows from a standard change of variables 
\[ 
\int_\mathbb{R} F(s) \, ds = \int_{-\infty}^0 F(s) \, ds + \int_0^{\infty}  
F(s) \, ds = \int_0^{\infty} F(s) \, ds - \int_0^{\infty} |  
F(-s)| \, ds .   
\]  
The two observations above imply the inequality (\ref{negeq}) holds, 
in which case by (\ref{eq:uprimeform2}) and (\ref{eq:poseq}) we 
see $ u'(t) > 0 $.  
Since $t > 0$ was arbitrary,  the proof is complete. 
\hfill$\diamondsuit$\medskip

\begin{corollary}
The solution to the initial value problem (\ref{eq:gen1})-(\ref{eq:gen2}) 
monotonically approaches its steady state value $u=1$.  
\end{corollary}
\paragraph{Proof.} 
The proof follows from applying Theorem 3.1 to equation (\ref{220}).
\hfill$\diamondsuit$\medskip
 
 
\section{Related Aspects of the Newtonian Problem} 
 
To investigate the generality of the above result, consider
the nonautonomous linear damped harmonic oscillator equation
for $u=u(t)$
\begin{equation} 
u'' + bu' + u = 1 - G(t),\label{eq:sho} 
\end{equation}
as an initial value problem with arbitrary initial conditions
$u(0)$   and $u'(0)$.  We are specifically interested in the case
where $G(t) \rightarrow 0$ as $t \rightarrow \infty$, 
as opposed to the often studied case where $G(t)$ is periodic 
(see e.g. \cite{hale}).  The Newtonian  sphere problem 
(\ref{eq:non-dim1}) is a special case of (\ref{eq:sho}), with 
$b =   2-\kappa$ and $G(t) = \sqrt{\kappa/\pi t}.$
Making the change of   variables $v = u-1$, we may simplify the 
equation to 
\begin{equation} 
v'' + bv' + v = - G(t)\label{eq:shov} 
\end{equation} 
so that $v=0$ solves the homogeneous equation.  Note however that if  
$G(t) \neq 0$, then $v=0$ is not a solution to (\ref{eq:shov}) for
 any $t > 0$. 
We are interested in the following question:  what  conditions  
on $G(t)$ and $b$ are necessary for a solution   $v(t)$ to remain monotone,
even within the regime of instability for the homogeneous equation.

As a first step in this direction we consider the following initial value problem  
for $t\geq0:$ 
\begin{equation} 
 v'' + b \, v' + v = -\frac{A}{\sqrt{\pi(t + t_0)}}, 
\label{eq:odealg} 
\end{equation} 
where $b,A \in \mathbb{R}$ and $t_0 \geq 0$ are constants.  
The motivation for this form is to test the 
necessity of the singularity at $t=0$ in the monotonicity result of  
Section \ref{monosection}.  To ensure complex roots, 
we assume $b \in (-2,2)$.
 
Using variation of parameters one finds a particular solution of  
(\ref{eq:odealg}) to be 
 \begin{eqnarray} 
v_p(t) &=&  \frac{A}{\beta - \alpha} \left\{ \sqrt{\beta}  
 e^{\alpha (t+t_0)} \left( \text{Erfc} \sqrt{\alpha t_0} - \text{Erfc}\sqrt{\alpha(t+t_0)} \right) \right. \nonumber \\  
&& \hspace{0.35in} 
-\left.\sqrt{\alpha}   e^{\beta (t+t_0)} \left( \text{Erfc} \sqrt{\beta t_0} -  
\text{Erfc}\sqrt{\beta(t+t_0)}\right) \right\},   
\label{up} 
\end{eqnarray} 
where $\alpha$ and $\beta$ are the roots of the characteristic polynomial 
$ m^2 + b  m + 1.$ Employing the Villat function we may 
express equation (\ref{up}) as 
\[  
v_p(t) =  \frac{A}{\beta - \alpha} \left\{ \sqrt{\beta} \text{Vi}(\alpha  
t_0) e^{\alpha t} 
 - \sqrt{\alpha} \text{Vi}(\beta t_0) e^{\beta t} \right\} + A\,M(t+t_0),  \] 
where  the function $M$ is defined by
\begin{equation}  
M(t) = \frac{1}{\alpha - \beta} \left\{\sqrt{\beta} \text{Vi}(\alpha t) -  
\sqrt{\alpha} \text{Vi}(\beta t)\right\}.  
\end{equation}
In Section \ref{monosection} we proved $M$ approaches 0 monotonically for all 
$b \in (-2,2)$. The general solution to (\ref{eq:odealg}) may be expressed
as
\begin{eqnarray} 
 v(t) & = & C_1 e^{\alpha t} + C_2 e^{\beta t} + v_p(t)  \label{prefsoln}   
\\ 
  & = &    \left\{ C_1 +  \frac{\sqrt{\beta} A \text{Vi}(\alpha t_0)}{\beta - \alpha} 
\right\} e^{\alpha t} 
      + \left\{  C_2 -   \frac{\sqrt{\alpha} A \text{Vi}(\beta t_0)}{\beta - \alpha}  
\right\} e^{\beta t} \nonumber \\
  &  & 
 + \, A \, M(t+t_0).  \label{fsoln} 
\end{eqnarray} 
This equation clearly demonstrates how the long term dynamics of $v(t)$ depend on
the solution of the homogeneous problem.   In particular, it shows that the solution
will retain the 
 stability properties of the homogeneous solution unless the coefficients $C_1$ and 
$C_2$ are  chosen to zero out the first  two terms in (\ref{fsoln}).  The unique choice 
of $C_1$ and  $C_2$ for this to happen is
\begin{equation} 
   C_1 = -\frac{A\sqrt{\beta}}{\beta - \alpha} \text{Vi}(\alpha t_0) \qquad  
\text{and } \qquad 
        C_2 = \frac{A\sqrt{\alpha}}{\beta - \alpha} \text{Vi}(\beta t_0). \label{ab} 
\end{equation} 
Moreover, it is clear from (\ref{fsoln}) that the solution in this case is  
$v(t)=A\,M(t+t_0)$, with $v(0)=A\,M(t_0)$ and  
$v'(0)=A\,M'(t_0)$.  Thus, in this case, the solution is
a translate of the monotone solution.   
The coefficients $C_1$ and $C_2$ are related to the initial conditions  
$v(0)$ and $v'(0)$ via 
\begin{equation} 
  C_1 = \frac{\beta \, v(0) - v'(0)}{\beta - \alpha}  \qquad  
\text{and } \qquad   
C_2 = \frac{v'(0) - \alpha \, v(0)}{\beta - \alpha}. \label{a0b0} 
\end{equation} 
The values of $v(0)=A\,M(t_0)$ and $v'(0)=A\,M'(t_0)$ may also be  
obtained by solving equations (\ref{ab}) and (\ref{a0b0}).
 
In summary, given $b \in (-2,2)$, $A > 0$, and $t_0 \geq 0$, for the  
equation 
\begin{equation}
v'' + b \, v' + v = -\frac{A}{\sqrt{\pi(t + t_0)}}, 
\label{gensho}
\end{equation}
there exists a unique choice of initial values $v(0)=A\,M(t_0)$,  
$v'(0)=A \, M'(t_0)$ such that the solution $v(t)$ remains monotone for all  
$t > 0$. Therefore the presence of a singularity at $t=0$ in the 
nonhomogeneous term is not necessary to obtain a monotone solution. 

In light of the above analysis it becomes clear how the solution for  the 
sedimenting sphere remains monotone in its approach to terminal velocity for 
\textit{all} relevant values of $\kappa$ (i.e., sphere densities).  
 From equations (\ref{eq:non-dim1})-(\ref{eq:non-dim2}) we see
that for each value of $\kappa$, equation 
(\ref{gensho}) describes the dynamics for the dimensionless velocity
$v=u-1$, with $b =2-\kappa$, $t_0=0$, and $A=\sqrt{\kappa}$. 
Moreover, for each $\kappa \in (0,4)$  we have demonstrated that
equation (\ref{gensho}) with 
$t_0=0$, $b=2-\kappa$, and $A=\sqrt{\kappa}$, has
a unique initial value for which the solution remains monotone, namely,
\begin{equation}
  v(0)= A M(0) =  \sqrt{\kappa} \, \frac{\sqrt{\beta} - \sqrt{\alpha}}{\alpha - \beta}
 =-\frac{\sqrt{\kappa}}{\sqrt{\alpha}+\sqrt{\beta}},
\label{kd}
\end{equation}
where $\alpha$ and $\beta$ denote the roots of the polynomial $m^2 + (2-\kappa)m + 1$.
However, since $\sqrt{\alpha}$ lies in the first quadrant, $\sqrt{\alpha}
+ \sqrt{\beta} > 0$, and the computation
 \[ (\sqrt{\alpha} + \sqrt{\beta})^2 = \alpha + \beta + 2 = -b + 2 = \kappa,  \]
together with (\ref{kd}), implies 
\[  v(0)= -1.  \]
In other words,
the particular relation between the parameters $A$ and $b$
decouples $v(0)$ from all parameters, so that one obtains a monotone solution for all
values of the sphere density. 

We conclude this section with a geometric interpretation of the
monotonicity result.  In particular, we focus on  the interesting 
case of (\ref{eq:odealg}) when the parameter $b\in(-2,0)$. For 
these parameter values the solution $v=0$ of the 
homogeneous problem is unstable.   We have shown that
there is a unique set of initial conditions that defines a solution to the 
nonhomogeneous problem which approaches the unstable fixed point $v=0$
monotonically, despite the surrounding instabilities. 
Thus we return to the  fundamental puzzle posed in Section \ref{asympt}: 
How is it that the nonautonomous term in  (\ref{eq:odealg}), which 
decays to zero as $t \rightarrow \infty$, can ``stabilize'' a trajectory for 
all $ t > 0$, in the sense that this solution approaches 0 while all other
trajectories diverge due to 
the  instability of the linearized problem?
The following observation resolves the puzzle. 
First, recall that  the unique initial conditions for which the 
nonhomogeneous problem remains monotone  are  defined by
\[ v(0) = A \, M(t_0) \qquad \text{and} \qquad v'(0)= A M'(t_0).  \]  
Second, note that as the  amplitude $A$ of the nonhomogeneous term
tends to zero,  the  initial conditions  
$\left(A\,M(t_0),A\,M'(t_0)\right)$  approach $(0,0)$. This 
corresponds to the initial condition starting on the unstable equilibrium
point, which is the unique initial condition for the homogeneous
problem  whose solution does not diverge. 
In other words we have a correspondence between the trajectories of the
homogeneous equation   and the nonhomogeneous equation,  which is
continuous with respect to the  parameter $A$. The monotone
solution is then the image of the unstable fixed point under this
map. 

\subsection*{Summary and Conclusion}

In this paper we have studied the ODE model for a sphere falling 
through a Newtonian fluid.  We have proven that the equations do not admit 
oscillations, even in the transient, in agreement with   general
experimental observations.  

 From our analysis it appears that the lack of oscillations is
due to a delicate balance of terms.  It is tempting
to conclude that an oscillating motion could be produced with only a
slight modification to the equations.  However it is important that
the solution  still remain bounded, and as we have shown there is
only one trajectory which is insensitive to the linear instability
($\Re(\alpha)>0$) of the homogeneous equation. Transient oscillations
of a falling sphere have been successfully modeled by King \& Waters
using an elastic constitutive model  \cite{king:ums72}, for which a
final steady state velocity is approached.  

In principle, however, one cannot simply modify the differential
equation (\ref{eq:non-dim1}) or even (\ref{eq:maineq}) to address the
oscillations of a sedimenting sphere in a micellar fluid
\cite{belmonte:soc00,jayaraman:osf00}; one must return to the full
time-dependent partial differential equation. This was indeed how
King \& Waters obtained their result for a linear  viscoelastic
constitutive model   \cite{king:ums72}, but it is not clear that this
approach will continue to be fruitful as the complexity of the
problem increases. Self-assembling wormlike micellar solutions  are
thought to have a nonmonotonic stress/shear rate relation
\cite{spenley93,porte97}, based on the existence of an apparently
inaccessible range of shear
rates \cite{porte97,cappel97}. 
It may be that the dynamics of such a nonlinear fluid requires the
spatial information inherent in the PDEs, and that the ODE reduction
discussed here is practically limited to linear models. 

 
\section{Appendix: Derivation of Eqns. (\ref{eq:non-dim1})-(\ref{eq:non-dim2})} 
\label{app}

The  equation describing the transient motion of a falling sphere is
\begin{equation}
u' + u + \sqrt{\frac{\kappa}{\pi}} \int_0^t \frac{u' (s)}{\sqrt{t-s}} ds =1 \label{eq:1},
\end{equation}
where $u(t)$ is the velocity of the sphere and $\kappa$ is 
a non-dimensional parameter which depends  on the relative densities
of the sphere and the fluid. 
This integro-differential equation can be converted to a second order ODE through the 
following procedure. If
$$F(t) = \int_0^t  \frac{u' (s)}{\sqrt{t-s}} ds, $$
then Abel's theorem (see e.g., \cite[\S3.7]{keener:pam00}) implies
\begin{equation}
\int_0^t  \frac{F(\tau)}{\sqrt{t-\tau}} \, d\tau = \pi \left[u(t)-u(0)\right].
\label{abthr}
\end{equation}
Multiplying (\ref{eq:1}) by $1/\sqrt{t-\tau},$ integrating, and using (\ref{abthr})
yields the equation 
\begin{equation}
  \int_0^t \frac{u'(\tau)}{\sqrt{t-\tau}} \, d\tau +  \int_0^t \frac{ u(\tau)}{\sqrt{t-\tau}} \, 
d\tau + \pi \sqrt{\frac{\kappa}{\pi}} \left[u(t) - u(0)\right] =   
\int_0^t \frac{1}{\sqrt{t-\tau}} \, d\tau.
 \label{eq:2}
\end{equation}
 From (\ref{eq:1}) one observes
\begin{equation}
 \int_0^t \frac{u'(\tau)}{\sqrt{t-\tau}} \, d\tau = \sqrt{\frac{\pi}{\kappa}} 
\left(1- u - u'\right).\label{eq:n}
\end{equation}
Substituting this into (\ref{eq:2}) and rewriting yields
\begin{equation}
u' = \Big(1- 2\sqrt{\frac{\kappa t}{\pi}}\Big) +\left(\kappa -1\right) u - 
\kappa u(0)+ \sqrt{\frac{\kappa}{\pi}} \int_0^t \frac{ u(\tau)}{\sqrt{t-\tau}} \, d\tau. 
\label{eq:4}
\end{equation}
The desired second order differential equation is now obtained by differentiating
(\ref{eq:4}). In this regard, note that the substitution $\tau=t-x^2$ implies
\begin{equation}
I(t) = \sqrt{\frac{\kappa}{\pi}} \int_0^t \frac{u(\tau)}{\sqrt{t-\tau}} \, d\tau = 
2\sqrt{\frac{\kappa}{\pi}} \int_0^{\sqrt{t}} u(t-x^2) \, dx, \label{eq:5}
\end{equation}
thus
\begin{eqnarray}
\frac{dI}{dt}&=& \sqrt{\frac{\kappa}{\pi}}\frac{u(0)}{\sqrt{t}} + 
2\sqrt{\frac{\kappa}{\pi}}\int_0^{\sqrt{t}}  u'(t-x^2) \, dx \nonumber\\
&=& \sqrt{\frac{\kappa}{\pi}}\frac{u(0)}{\sqrt{t}} +  \sqrt{\frac{\kappa}{\pi}} 
\int_0^t \frac{u'}{\sqrt{t-\tau}} \, d\tau\nonumber \\
&=&  \sqrt{\frac{\kappa}{\pi}}\frac{u(0)}{\sqrt{t}} + (1- u - u'), \label{lst}
\end{eqnarray}
where again we have used (\ref{eq:n}).

Therefore differentiating (\ref{eq:4}) and using (\ref{lst}) yields the second order 
equation
\begin{equation}
u'' = (\kappa-2) u' - u + \Big[1+\sqrt{\frac{\kappa}{\pi t}}( u(0)-1)\Big].
\end{equation}
Note that from (\ref{eq:1}) the initial value of $u'$ is determined by the initial
value of $u$, i.e.,  $ u' (0) = 1 - u(0)$.
Therefore, the equation describing the transient motion of the sphere is
\begin{gather} 
u'' + (2-\kappa) u' + u  =  1 + \sqrt{\frac{\kappa}{\pi  
t}}(u(0) - 1), 
\label{eq:genn1} \\ 
u(0)  =  \xi, \quad u'(0)  =  1 - \xi. \label{eq:genn2}  
\end{gather} 
If the sphere starts from rest (i.e., $u(0)=0$) then the system reduces to
\begin{gather} 
u'' + (2-\kappa) u' + u = 1- \sqrt{\frac{\kappa}{\pi  
\tau}},   \nonumber \\ 
u(0) = 0, \quad u'(0) = 1,  \nonumber 
\end{gather} 
which is precisely (\ref{eq:non-dim1})-(\ref{eq:non-dim2}).

 
\paragraph{Acknowledgments} The authors would like to thank J. P. 
Keener, H. A. Stone, B. Ermentrout, and W. Zhang for helpful
discussions.   A.~B.~acknowledges the support of the Alfred P.~Sloan
Foundation.



\begin{thebibliography}{00} \frenchspacing

\bibitem{arigo97} {\sc M.~T. Arigo and G.~H. McKinley}, \textit{The effects of  
viscoelasticity on the transient motion of a sphere in a  
shear-thinning fluid}, J. Rheol., 41 (1997), pp. 103-128.

\bibitem{basset} {\sc A. B. Basset}, {\em  A Treatise on
Hydrodynamics}, Dover, New York, NY, 1961, Volume 2 (reprint of 1888
edition). 

 \bibitem{belmonte:soc00} {\sc A. Belmonte}, {\em Self-oscillations of  
a cusped bubble rising through a micellar solution}, Rheol. Acta, 39 (2000), pp. 554-559.
 
\bibitem{berret94} {\sc  J.~F. Berret, D. Roux, and G. Porte}, \textit{Isotropic-to-nematic  
transition in wormlike micelles under shear}, J. Physique II, 4 (1994), pp. 1261-1279. 
 
\bibitem{bird87} {\sc R. Bird, R. Armstrong, and O. Hassager},  
\textit{Dynamics of Polymeric Liquids}, Vol. 1, Second ed.,  Wiley and Sons,  
New York, NY, 1987.  
 
\bibitem{cappel97} {\sc E. Cappelaere, J.~F. Berret, J.~P. Decruppe, R. Cressely, and  
P. Lindner}, \textit{Rheology, birefringence, and small-angle neutron  
scattering in a charged micellar system: Evidence of a shear-induced  
phase transition}, Phys. Rev. E, 56 (1997), pp. 1869-1878. 
 
\bibitem{clift} {\sc R. Clift, J.~R. Grace, and M.~E. Weber}, {\em  
Bubbles, drops, and particles}, Academic Press, New York, NY, 1978. 
 
\bibitem{estrada:sie00} {\sc R. Estrada and R.~P. Kanwal}, {\em  
Singular Integral Equations}, Birkhauser, Boston, MA, 2000. 
 
\bibitem{gautschi:eff92} {\sc W. Gautschi}, \textit{Error Function
and Fresnel Integrals}, in Handbook of  Mathematical Functions with
Formulas, Graphs, and  Mathematical Tables, M. Abramowitz and I.~A.
Stegun, eds.,  Dover Publications, Inc.,  New York, NY,  1992
(reprint of the 1972 edition), pp. 295-329. 

\bibitem{hale} {\sc J.~K. Hale and H. Kocak}, {\em  
Dynamics and Bifurcations}, Springer-Verlag, New York, NY, 1991. 
 
\bibitem{hass79} {\sc O Hassager}, \textit{Negative wake behind bubbles in  
non-Newtonian liquids}, Nature, 279 (1979), pp.  402-403. 
 
\bibitem{hu98} {\sc Y.~T. Hu, P. Boltenhagen, and D.~J. Pine}, {\em 
Shear thickening in low-concentration solutions of wormlike micelles. 
I. Direct visualization of transient behavior and phase 
transitions},  J. Rheol., 42 (1998), pp. 1185--1208. 
 
\bibitem{jayaraman:osf00} {\sc A. Jayaraman and A. Belmonte}, {\em  
Oscillations of a sphere falling through a micellar solution},  
preprint. 
 
\bibitem{keener:pam00} {\sc J.~P. Keener}, {\em Principles of Applied  
Mathematics: Transformation and Approximation}, Second ed., Perseus  
Books, Cambridge, MA, 2000. 
 
\bibitem{king:ums72} {\sc M.~J. King and N.~D. Waters}, {\em The  
unsteady motion of a sphere in an elastico-viscous liquid}, J. Phys.   
D: Appl.  Phys., 5 (1972), pp. 141--150. 

\bibitem{lamb} {\sc H. Lamb}, {\em   Hydrodynamics}, Dover, New York,
NY, 1945 (reprint of 1932 edition). 

 
\bibitem{liu95} {\sc Y. Liu, T. Liao, and D.~D. Joseph}, \textit{A two-dimensional cusp at  
the trailing edge of an air bubble rising in a viscoelastic liquid},  
J. Fluid Mech., 304 (1995), p. 321. 
%pp. 321-??.
 

\bibitem{michael97} {\sc E.~E. Michaelides}, {\em The Transient 
Equation of Motion for Particles, Bubbles, and Droplets}, J. Fluids 
Eng: Trans. ASME, 119 (1997), pp. 233--247. 
 
\bibitem{ock68} {\sc J.~R. Ockendon}, {\em The unsteady motion of a 
small sphere in a viscous fluid}, J. Fluid Mech., 34 (1968), 
pp. 229--239. 

\bibitem{porte97} {\sc G. Porte, J.~F. Berret, and J.~L. Harden}, \textit{Inhomogeneous flows  
of complex fluids: Mechanical instability versus non-equilibrium phase  
transition}, J. Physique II, 7 (1997), pp. 459-472. 

\bibitem{rehage82} {\sc H. Rehage  and H.  Hoffmann}, \textit{Shear induced phase  
transitions in highly dilute aqueous detergent solutions},   
Rheol.  Acta,  21 (1982), p. 561. 
%pp. 561-??.
   
\bibitem{spenley93} {\sc N. Spenley, M. Cates, and T. McLeish}, \textit{Nonlinear  
rheology of wormlike micelles},  Phys.  Rev.  Lett.,  71,  
(1993), pp. 939-942. 

\bibitem{stokes:51} {\sc G.~G. Stokes}, \textit{On the effect of the
internal friction of fluids on the motion of a pendulum},   Trans.
Cambridge Philos. Soc.,  9 (1851), pp. 8-106.

\bibitem{villat:44} {\sc H. Villat},  \textit{Le\c{c}ons sur les 
Fluides Visqueux}, Gauthier-Villars, Paris, 1944. 


\bibitem{walt92} {\sc K. Walters and R.~I. Tanner}, {\em The
motion of a sphere through an elastic fluid}, in R. P. Chhabra
and D. De Kee, eds., Transport Processes in Bubbles, Drops, and
Particles, Hemisphere, New York, NY, 1992. 

\bibitem{wil78} {\sc D.~G. Wilson}, {\em Existence and uniqueness for  
similarity solutions of one dimensional multi-phase Stefan problems},  
SIAM J. Appl. Math., 35 (1978), pp. 135--147. 
\end{thebibliography} 

\noindent {\sc Andrew Belmonte} \\
The W. G. Pritchard Laboratories \\
Department of Mathematics \\
Penn State University \\
University Park, PA 16802 USA \\ 
e-mail: belmonte@math.psu.edu 
\smallskip

\noindent {\sc Jon Jacobsen}\\
The W. G. Pritchard Laboratories \\
Department of Mathematics\\
Penn State University \\
University Park, PA 16802 USA\\ 
e-mail: jacobsen@math.psu.edu 
\smallskip

\noindent {\sc Anandhan Jayaraman} \\
The W. G. Pritchard Laboratories \\
Department of Mathematics \\
Penn State University \\
University Park, PA 16802 USA \\ 
e-mail: anand@math.psu.edu 


\end{document}