text:
\begin{document}
\title[The Degenerate Parametric Oscillator ]{The Degenerate Parametric
Oscillator and Ince's Equation}
\author{Ricardo Cordero-Soto}
\address{Mathematical, Computational and Modeling Sciences Center, Arizona
State University, Tempe, AZ 85287--1904, U.S.A.}
\email{ricardojavier81@gmail.com}
\author{Sergei K. Suslov}
\address{Mathematical, Computational and Modeling Sciences Center, Arizona
State University, Tempe, AZ 85287--1904, U.S.A.}
\email{sks@asu.edu}
\urladdr{http://hahn.la.asu.edu/\symbol{126}suslov/index.html}
\date{\today }
\subjclass{Primary 81Q05, 35C05. Secondary 42A38}
\keywords{The time-dependent Schr\"{o}dinger equation, generalized harmonic
oscillators, degenerate parametric amplification, Green's function,
propagator, Ince's equation, Ermakov's equation.}

\begin{abstract}
We construct Green's function for the quantum degenerate parametric
oscillator in terms of standard solutions of Ince's equation in a framework
of a general approach to harmonic oscillators. Exact time-dependent wave
functions and their connections with dynamical invariants and $SU\left(
1,1\right) $ group are also discussed.
\end{abstract}

\maketitle

\section{The Degenerate Parametric Amplifier Hamiltonian}

We use standard annihilation and creation operators in the coordinate
representation \cite{Flu}, \cite{La:Lif}, \cite{Merz} given by
\begin{equation}
\widehat{a}=\sqrt{\frac{m}{2\hslash \omega }}\left( \omega x+\frac{ip}{m}
\right) ,\qquad \widehat{a}^{\dagger }=\sqrt{\frac{m}{2\hslash \omega }}
\left( \omega x-\frac{ip}{m}\right)   \label{at1}
\end{equation}
with $p=-i\hslash \partial /\partial x$ and $\widehat{a}\widehat{a}^{\dagger
}-\widehat{a}^{\dagger }\widehat{a}=1$ throughout this Letter. The following
Hamiltonian
\begin{equation}
H\left( t\right) =\frac{\hslash \omega }{2}\left( \widehat{a}\ \widehat{a}
^{\dagger }+\widehat{a}^{\dagger }\ \widehat{a}\right) -\frac{\hslash
\lambda }{2}\left( e^{2i\omega t}\ \widehat{a}^{2}+e^{-2i\omega t}\left(
\widehat{a}^{\dagger }\right) ^{2}\right)   \label{at2}
\end{equation}
of the time-dependent Schr\"{o}dinger equation
\begin{equation}
i\hslash \frac{\partial \psi }{\partial t}=H\left( t\right) \psi
\label{at2a}
\end{equation}
describes the process of \ degenerate parametric amplification in quantum
optics. The first term corresponds to the self-energy of the oscillator
representing the mode of interest, the second term describes the coupling of
the classical pump to that mode (with the phenomenological constant $\lambda
),$ giving rise to the parametric amplification process (the phase of the
pump is taken to be zero at $t=0$ for convenience). Takahasi \cite{Taka65}
wrote a detailed paper on this subject as far as in 1965 (see also \cite
{Louisel64}, \cite{Orszag}, \cite{Per:Hra:Jur}, \cite{Puri}, \cite{Raiford70}
, \cite{Raiford74}, \cite{Shen}, \cite{Stoler74}, \cite{Walls:Milburn} and
references therein). The same Hamiltonian had also been considered by
Angelow and Trifonov in the mid 1990s in order to describe the light
propagation in the nonlinear $Ti:LiNbO_{3}$ anisotropic waveguide (details
can be found in \cite{Angelow98} and \cite{Angelow:Trifonov95}). Another
form of this variable quadratic Hamiltonian is given by
\begin{eqnarray}
H &=&\frac{1}{2m}\left( 1+\frac{\lambda }{\omega }\cos \left( 2\omega
t\right) \right) p^{2}+\frac{m\omega ^{2}}{2}\left( 1-\frac{\lambda }{\omega
}\cos \left( 2\omega t\right) \right) x^{2}  \label{at3} \\
&&+\frac{\lambda }{2}\sin \left( 2\omega t\right) \left( px+xp\right)
,\qquad p=\frac{\hslash }{i}\frac{\partial }{\partial x}.  \notag
\end{eqnarray}
We refer to (\ref{at2})--(\ref{at2a}) and/or (\ref{at3}) as the degenerate
parametric oscillator \cite{Orszag}. It can be recognized as a special case
of so-called generalized harmonic oscillators that had attracted
considerable attention over the years in view of their great importance to
many advanced quantum problems including coherent states and uncertainty
relations, Berry's phase, asymptotic and numerical methods, quantization of
mechanical systems, Hamiltonian cosmology, charged particle traps and motion
in uniform magnetic fields, molecular spectroscopy and polyatomic molecules
in varying external fields. Quadratic Hamiltonians have particular
applications in quantum electrodynamics because the electromagnetic field
can be represented as a set of forced harmonic oscillators. See, for
example, \cite{Berry85}, \cite{Cor-Sot:Sua:SusInv}, \cite{Dod:Mal:Man75},
\cite{Dodonov:Man'koFIAN87}, \cite{FeynmanPhD}, \cite{Fey:Hib}, \cite
{GomezVergel:Villasenor09}, \cite{Hannay85}, \cite{Leach90}, \cite
{Malkin:Man'ko79}, \cite{Menouar:Maamache:Choi10}, \cite{Suslov10}, \cite
{Wolf81}, \cite{Yeon:Lee:Um:George:Pandey93} and references therein.
Nonlinear oscillators play a central role in the theory of Bose--Einstein
condensation \cite{Dal:Giorg:Pitaevski:Str99} because the dynamics of gases
of cooled atoms in a magnetic trap at very low temperatures can be described
by an effective equation for the condensate wave function known as the
Gross--Pitaevskii (or nonlinear Schr\"{o}dinger) equation \cite
{Kagan:Surkov:Shlyap96}, \cite{Kagan:Surkov:Shlyap97}, \cite
{Kivsh:Alex:Tur01}, \cite{Per-G:Tor:Mont} and \cite
{Trall-Gin:Drake:Lop-Rich:Trall-Herr:Bir}.\medskip

A goal of this Letter is to construct Green's function of the Hamiltonian (
\ref{at3}) in the coordinate representation in terms of solutions of Ince's
equation studied in \cite{Mag:Win} and \cite{Menn68} (see also references
therein). Exact time-dependent wave functions and connections with linear
and quadratic dynamical invariants and $SU\left( 1,1\right) $ group are also
briefly discussed.

\section{Generalized Harmonic Oscillators}

The degenerate parametric oscillator Hamiltonian (\ref{at3}) belongs to
quantum systems described by the one-dimensional time-dependent Schr\"{o}
dinger equation
\begin{equation}
i\frac{\partial \psi }{\partial t}=H\psi  \label{Schroedinger}
\end{equation}
with variable quadratic Hamiltonians of the form
\begin{equation}
H=a\left( t\right) p^{2}+b\left( t\right) x^{2}+c\left( t\right) px+d\left(
t\right) xp,\quad p=-i\frac{\partial }{\partial x},  \label{GenHam}
\end{equation}
where $a\left( t\right) ,$ $b\left( t\right) ,$ $c\left( t\right) ,$ and $
d\left( t\right) $ are real-valued functions of time, $t,$ only (see, for
example, \cite{Cor-Sot:Lop:Sua:Sus}, \cite{Cor-Sot:Sua:Sus}, \cite
{Cor-Sot:Sua:SusInv}, \cite{Cor-Sot:Sus}, \cite{Dod:Mal:Man75}, \cite
{FeynmanPhD}, \cite{Feynman}, \cite{Fey:Hib}, \cite{Lo93}, \cite{Me:Co:Su},
\cite{Suaz:Sus}, \cite{Wolf81} and \cite{Yeon:Lee:Um:George:Pandey93} for a
general approach and known elementary solutions; a case related to Airy
functions is discussed in \cite{Lan:Sus} and our note deals with another
special case of transcendental solutions).\medskip

The corresponding Green function, or Feynman's propagator, can be found as
follows \cite{Cor-Sot:Lop:Sua:Sus}, \cite{Suaz:Sus}:
\begin{equation}
\psi =G\left( x,y,t\right) =\frac{1}{\sqrt{2\pi i\mu _{0}\left( t\right) }}\
e^{i\left( \alpha _{0}\left( t\right) x^{2}+\beta _{0}\left( t\right)
xy+\gamma _{0}\left( t\right) y^{2}\right) },  \label{in2}
\end{equation}
where
\begin{eqnarray}
&&\alpha _{0}\left( t\right) =\frac{1}{4a\left( t\right) }\frac{\mu
_{0}^{\prime }\left( t\right) }{\mu _{0}\left( t\right) }-\frac{c\left(
t\right) }{2a\left( t\right) },  \label{in3} \\
&&\beta _{0}\left( t\right) =-\frac{\lambda \left( t\right) }{\mu _{0}\left(
t\right) },\qquad \lambda \left( t\right) =\exp \left( \int_{0}^{t}\left(
c\left( s\right) -d\left( s\right) \right) \ ds\right) ,  \label{in4} \\
&&\gamma _{0}\left( t\right) =\frac{1}{2\mu _{1}\left( 0\right) }\frac{\mu
_{1}\left( t\right) }{\mu _{0}\left( t\right) }+\frac{c\left( 0\right) }{
2a\left( 0\right) }  \label{in5}
\end{eqnarray}
and $\mu _{0}\left( t\right) $ and $\mu _{1}\left( t\right) $ are two
linearly independent solutions of the\ characteristic equation
\begin{equation}
\mu ^{\prime \prime }-\tau \left( t\right) \mu ^{\prime }+4\sigma \left(
t\right) \mu =0  \label{in6}
\end{equation}
with
\begin{equation}
\tau \left( t\right) =\frac{a^{\prime }}{a}+2c-2d,\qquad \sigma \left(
t\right) =ab-cd+\frac{c}{2}\left( \frac{a^{\prime }}{a}-\frac{c^{\prime }}{c}
\right)  \label{in7}
\end{equation}
subject to initial conditions
\begin{eqnarray}
&&\mu _{0}\left( 0\right) =0,\qquad \mu _{0}^{\prime }\left( 0\right)
=2a\left( 0\right) \neq 0,  \label{in8} \\
&&\mu _{1}\left( 0\right) \neq 0,\qquad \mu _{1}^{\prime }\left( 0\right) =0.
\label{in9}
\end{eqnarray}
In this note, we present the time-dependent coefficient $\gamma _{0}\left(
t\right) $ in terms of two standard solutions $\mu _{0}$ and $\mu _{1}$ of
the characteristic equation (\ref{in6})--(\ref{in7}) with initial data (\ref
{in8})--(\ref{in9}), respectively, instead of a previously used formula in
terms of only one solution $\mu _{0}$ (see, for example, \cite
{Cor-Sot:Lop:Sua:Sus}, \cite{Cor-Sot:Sua:SusInv} and \cite{Suslov10}). Our
expression (\ref{in5}) can be verified by a direct differentiation:
\begin{equation}
\left( \frac{\mu _{1}}{\mu _{0}}\right) ^{\prime }=\frac{\mu _{1}^{\prime
}\mu _{0}-\mu _{1}\mu _{0}^{\prime }}{\mu _{0}^{2}}=\frac{W\left( \mu
_{0},\mu _{1}\right) }{\mu _{0}^{2}},  \label{in10}
\end{equation}
where the Wronskian can be found with the help of Abel's theorem as follows
\begin{equation}
W\left( \mu _{0},\mu _{1}\right) =\text{constant\ }\lambda ^{2}\left(
t\right) a\left( t\right) .  \label{in11}
\end{equation}
Then, as required,
\begin{equation}
\frac{d\gamma _{0}}{dt}+a\left( t\right) \beta _{0}^{2}\left( t\right) =0
\label{in12}
\end{equation}
and the constant term in (\ref{in5}) gives a correct asymptotic as $
t\rightarrow 0^{+}$ (cf.~Refs.~\cite{Cor-Sot:Sus} and \cite{Suaz:Sus}). The
Green function (\ref{in2}) is an eigenfunction of the linear dynamical
invariant of Dodonov, Malkin, Man'ko and Trifonov \cite{Dod:Mal:Man75}, \cite
{Dodonov:Man'koFIAN87}, \cite{Malkin:Man'ko79} and \cite{Malk:Man:Trif73}
(see also Theorem~1 of \cite{Suslov10}).\medskip

By the superposition principle the solution of the Cauchy initial value
problem can be presented in an integral form
\begin{equation}
\psi \left( x,t\right) =\int_{-\infty }^{\infty }G\left( x,y,t\right) \ \chi
\left( y\right) \ dy,\quad \lim_{t\rightarrow 0^{+}}\psi \left( x,t\right)
=\chi \left( x\right)  \label{CauchyInVProb}
\end{equation}
for a suitable initial function $\chi $ on $

\mathbb{R}

$ (a rigorous proof is given in \cite{Suaz:Sus} and uniqueness is analyzed
in \cite{Cor-Sot:Sua:SusInv}; another form of solution is provided by an
eigenfunction expansion \cite{Suslov10}; generalized coherent states and
transition amplitudes are discussed in a forthcoming paper \cite{Suslov10a}
). In the next sections, we apply these general results to the case of the
degenerate parametric oscillator (\ref{at3}).

\section{Green's Function and Ince Equation}

For the Hamiltonian (\ref{at3}), variable coefficients are
\begin{eqnarray}
&&a=\frac{\hslash }{2m}\left( 1+\frac{\lambda }{\omega }\cos \left( 2\omega
t\right) \right) ,  \label{ATa} \\
&&b=\frac{m\omega ^{2}}{2\hslash }\left( 1-\frac{\lambda }{\omega }\cos
\left( 2\omega t\right) \right) ,  \label{ATb} \\
&&c=d=\frac{\lambda }{2}\sin \left( 2\omega t\right)  \label{ATcd}
\end{eqnarray}
and general expressions for Green's function (\ref{in2})--(\ref{in5}) can be
simplified by letting $\lambda \left( t\right) \equiv 1$ and $c\left(
0\right) =0.$ The characteristic equation is given by
\begin{equation}
\mu ^{\prime \prime }+\frac{2\lambda \omega \sin \left( 2\omega t\right) }{
\omega +\lambda \cos \left( 2\omega t\right) }\mu ^{\prime }+\frac{\omega
\left( \omega ^{2}-3\lambda ^{2}\right) -\lambda \left( \omega ^{2}+\lambda
^{2}\right) \cos \left( 2\omega t\right) }{\omega +\lambda \cos \left(
2\omega t\right) }\mu =0,  \label{InceEqCh}
\end{equation}
which can be identified as a special case of the Ince equation \cite{Mag:Win}
, \cite{Menn68}:
\begin{equation}
\left( 1+a_{0}\cos 2s\right) y^{\prime \prime }\left( s\right) +b_{0}\sin
2s\ y^{\prime }\left( s\right) +\left( c_{0}+d_{0}\cos 2s\right) y\left(
s\right) =0,  \label{InceEq}
\end{equation}
when $s=\omega t$ and parameters are given by
\begin{equation}
a_{0}=\frac{\lambda }{\omega },\quad b_{0}=2a_{0}=2\frac{\lambda }{\omega }
,\quad c_{0}=1-3\frac{\lambda ^{2}}{\omega ^{2}},\quad d_{0}=-\frac{\lambda
}{\omega }\left( 1+\frac{\lambda ^{2}}{\omega ^{2}}\right) .
\label{InceParsCD}
\end{equation}

Traditionally, a special question which arises in the theory of Ince's
equation is the problem of the existence of periodic solutions. By
Theorem~7.1 on p.~93 of \cite{Mag:Win}, if Ince's equation (\ref{InceEq})
has two linearly independent solutions of period $\pi ,$ then the polynomial
\begin{equation}
P\left( \xi \right) =2a_{0}\xi ^{2}-b_{0}\xi -d_{0}/2  \label{PolQ}
\end{equation}
has a zero at one of the points $\xi =0,$ $\pm 1,$ $\pm 2,$ $...\ .$ If (\ref
{InceEq}) has two linearly independent solutions of period $2\pi ,$ then
\begin{equation}
Q\left( \xi \right) =2P\left( \xi -1/2\right)  \label{PlQSt}
\end{equation}
vanishes for one of the values $\xi =0,$ $\pm 1,$ $\pm 2,$ $...\ .$\medskip

In the case of the degenerate parametric oscillator Hamiltonian (\ref
{InceParsCD}), both of these polynomials are strictly positive:
\begin{equation}
P\left( \xi \right) =\frac{2\lambda }{\omega }\left( \left( \xi -\frac{1}{2}
\right) ^{2}+\frac{\lambda ^{2}}{4\omega ^{2}}\right) >0  \label{PolQin}
\end{equation}
and
\begin{equation}
Q\left( \xi \right) =\frac{4\lambda }{\omega }\left( \left( \xi -1\right)
^{2}+\frac{\lambda ^{2}}{4\omega ^{2}}\right) >0.  \label{PolQinSt}
\end{equation}
Thus the corresponding Ince's equation (\ref{InceEqCh}) may not have two
periodic solutions of periods $\pi /\omega $ or $2\pi /\omega $ and our
standard solutions $\mu _{0}$ and $\mu _{1}$ that satisfy initial conditions
(\ref{in8})--(\ref{in9}) cannot be constructed in terms of simple variants
of the Fourier sine and cosine series:
\begin{equation}
\mu _{1}\left( t\right) =\sum_{n=0}^{\infty }A_{2n}\cos \left( 2\omega
nt\right) ,\qquad \mu _{0}\left( t\right) =\sum_{n=1}^{\infty }B_{2n}\sin
\left( 2\omega nt\right) ,  \label{StSolI}
\end{equation}
or
\begin{equation}
\mu _{1}\left( t\right) =\sum_{n=0}^{\infty }A_{2n+1}\cos \left( \omega
\left( 2n+1\right) t\right) ,\qquad \mu _{0}\left( t\right)
=\sum_{n=0}^{\infty }B_{2n+1}\sin \left( \omega \left( 2n+1\right) t\right) ,
\label{StSolII}
\end{equation}
when coefficients $A_{m}$ and $B_{m}$ satisfy certain three-term recurrence
relations (more details can be found in \cite{Mag:Win}). Therefore, the
degenerate parametric oscillator (\ref{at3}) motivates a detailed study of
non-periodic solutions of Ince's equation.

\section{The Dynamical Invariant and Time-Dependent Wave Functions}

Exact time-dependent wave functions of generalized harmonic oscillators (\ref
{Schroedinger})--(\ref{GenHam}) can also be obtained as eigenfunctions of
the quadratic invariant in terms Hermite polynomials and certain solutions
of Ermakov-type equation (see Ref.~\cite{Suslov10} for a modern introduction
and references therein). The dynamical invariant of the Hamiltonian (\ref
{at3}) can be derived from a general form presented in \cite
{Cor-Sot:Sua:SusInv}, \cite{Suslov10} as follows:
\begin{equation}
E\left( t\right) =\left( \mu \ p+\frac{m}{\hslash }\frac{\lambda \sin \left(
2\omega t\right) \mu -\mu ^{\prime }}{1+\left( \lambda /\omega \right) \cos
\left( 2\omega t\right) }\ x\right) ^{2}+\frac{C_{0}}{\mu ^{2}}\ x^{2},\quad
p=-i\frac{\partial }{\partial x},  \label{QuadInv}
\end{equation}
where $C_{0}>0$ is a constant and function $\mu \left( t\right) $ satisfies
the nonlinear auxiliary equation:
\begin{eqnarray}
&&\mu ^{\prime \prime }+\frac{2\lambda \omega \sin \left( 2\omega t\right) }{
\omega +\lambda \cos \left( 2\omega t\right) }\mu ^{\prime }+\frac{\omega
\left( \omega ^{2}-3\lambda ^{2}\right) -\lambda \left( \omega ^{2}+\lambda
^{2}\right) \cos \left( 2\omega t\right) }{\omega +\lambda \cos \left(
2\omega t\right) }\mu  \notag \\
&&\qquad =C_{0}\left( \frac{\hslash }{m}\right) ^{2}\frac{\left( 1+\left(
\lambda /\omega \right) \cos \left( 2\omega t\right) \right) ^{2}}{\mu ^{3}}.
\label{InceErmEq}
\end{eqnarray}
A general solution of this nonlinear equation can be found by the following
`law of cosines':
\begin{equation}
\mu ^{2}\left( t\right) =Au^{2}\left( t\right) +Bv^{2}\left( t\right)
+2Cu\left( t\right) v\left( t\right)  \label{AuxSol}
\end{equation}
in terms of two linearly independent solutions $u$ and $v$ of the
homogeneous equation. The constant $C_{0}$ is related to the Wronskian of
two linearly independent solutions $u$ and $v:$
\begin{equation}
AB-C^{2}=C_{0}\frac{\left( 2a\right) ^{2}}{W^{2}\left( u,v\right) },\quad
W\left( u,v\right) =uv^{\prime }-u^{\prime }v  \label{AuxSolWronskian}
\end{equation}
(more details can be found in \cite{Suslov10}; see also references therein
regarding this so-called Pinney's solution). Thus solution of the
corresponding initial value problem is given by
\begin{equation}
\mu ^{2}\left( t\right) =\left( \frac{\mu ^{\prime }\left( 0\right) }{
2a\left( 0\right) }\mu _{0}\left( t\right) +\frac{\mu \left( 0\right) }{\mu
_{1}\left( 0\right) }\mu _{1}\left( t\right) \right) ^{2}+\frac{C_{0}}{\mu
^{2}\left( 0\right) }\mu _{0}^{2}\left( t\right) ,\quad \mu \left( 0\right)
\neq 0  \label{AuxSolIvp}
\end{equation}
in terms of standard solutions $\mu _{0}$ and $\mu _{1}$ corresponding to
initial data (\ref{in8})--(\ref{in9}) of the homogeneous equation (\ref
{InceErmEq}).\medskip

Time-dependent wave functions can be presented in the form
\begin{equation}
\psi _{n}\left( x,t\right) =e^{-i\left( n+1/2\right) \varphi \left( t\right)
}\ \Psi _{n}\left( x,t\right)   \label{WaveFuncs}
\end{equation}
(see, for example, \cite{Suslov10}, \cite{Yeon:Lee:Um:George:Pandey93} and
references therein), where
\begin{equation}
\frac{d\varphi }{dt}=\sqrt{C_{0}}\frac{\hslash }{m}\frac{1+\left( \lambda
/\omega \right) \cos \left( 2\omega t\right) }{\mu ^{2}},\quad \varphi
\left( 0\right) =0  \label{WavePhase}
\end{equation}
and orthonormal time-dependent eigenfunctions $\Psi _{n}\left( x,t\right) $
of the quadratic invariant (\ref{QuadInv}) are expressed in terms of Hermite
polynomials \cite{Ni:Su:Uv}:
\begin{eqnarray}
&&\Psi _{n}\left( x,t\right) =D_{n}\exp \left( ix^{2}\frac{m}{2\hslash }
\frac{\left( \mu ^{\prime }/\mu \right) -\lambda \sin \left( 2\omega
t\right) }{1+\left( \lambda /\omega \right) \cos \left( 2\omega t\right) }
\right)   \label{InceInvEigfs} \\
&&\qquad \qquad \times e^{-x^{2}\sqrt{C_{0}}/\left( 2\mu ^{2}\right)
}H_{n}\left( x\frac{C_{0}^{1/4}}{\mu }\right) ,\quad \left\vert
D_{n}\right\vert ^{2}=\frac{C_{0}^{1/4}}{\sqrt{\pi }2^{n}n!\mu }.  \notag
\end{eqnarray}
Then
\begin{equation}
i\frac{\partial \psi _{n}}{\partial t}=H\psi _{n},\qquad E\psi _{n}=2\sqrt{
C_{0}}\left( n+\frac{1}{2}\right) \psi _{n}.
\label{Schroedingereigenfunctions}
\end{equation}
If the orthonormal initial wave function is given by
\begin{equation}
\psi _{n}\left( x,0\right) =\sqrt{\frac{\varepsilon }{\sqrt{\pi }2^{n}n!}}\
e^{\left( i\delta -\varepsilon ^{2}/2\right) x^{2}}H_{n}\left( \varepsilon
x\right) ,  \label{InitFunc}
\end{equation}
where $\varepsilon $ and $\delta $ are constants, then initial conditions
for the auxiliary equation (\ref{InceErmEq}) are given by
\begin{equation}
\mu \left( 0\right) =\frac{C_{0}^{1/4}}{\varepsilon },\qquad \mu ^{\prime
}\left( 0\right) =2C_{0}^{1/4}\frac{\hslash }{m}\left( 1+\frac{\lambda }{
\omega }\right) \frac{\delta }{\varepsilon }  \label{InitConds}
\end{equation}
(the integral $C_{0}>0$ can be chosen at reader's convenience). Then
integration of Ermakov's equation (\ref{InceErmEq}), say with the help of
Pinney's solution (\ref{AuxSolIvp}), determines complete dynamics of the
oscillator-type wave function through the explicit formula (\ref
{InceInvEigfs}). Further relations of these particular solutions of the Schr
\"{o}dinger equation with the Cauchy initial value problem are discussed in
\cite{Suslov10} (see Theorem~2).

\section{The $SU\left( 1,1\right) $ Symmetry}

The Hamiltonian (\ref{at2}) can also be rewritten as follows
\begin{equation}
H\left( t\right) =\hslash \omega J_{0}-\hslash \lambda \left( e^{2i\omega
t}J_{+}+e^{-2i\omega t}J_{-}\right) ,  \label{sym1}
\end{equation}
where generators of a non-compact $SU\left( 1,1\right) $ group are given by
\begin{equation}
J_{+}=\frac{1}{2}\left( a^{\dagger }\right) ^{2},\qquad J_{-}=\frac{1}{2}
a^{2},\qquad J_{0}=\frac{1}{4}\left( a^{\dagger }a+a^{\dagger }a\right) .
\label{sym2}
\end{equation}
Then a use can be made of the group properties of the corresponding discrete
positive series $\mathcal{D}_{+}^{j}$ for further investigation of the
degenerate parametric oscillator. This is a `standard procedure' --- more
details can be found in Refs.~\cite{Lo93}, \cite{Malkin:Man'ko79}, \cite
{Me:Co:Su}, \cite{Ni:Su:Uv} and \cite{Smir:Shit} and/or elsewhere.

\section{A Special Case}

An integrable special case $\lambda =\omega =m=\hslash =1$ of the degenerate
parametric oscillator has been recently considered by Meiler, Cordero-Soto
and Suslov \cite{Me:Co:Su}, \cite{Cor-Sot:Sus} (this Hamiltonian is a
simplest time-dependent quadratic integral of motion for the linear harmonic
oscillator \cite{Cor-Sot:Sua:SusInv}). Here, the Ince equation (\ref
{InceEqCh}) simplifies to
\begin{equation}
\mu ^{\prime \prime }+2\tan t\ \mu ^{\prime }-2\mu =0.  \label{MeCo-SoSuEq}
\end{equation}
It has two elementary non-periodic standard solutions:
\begin{eqnarray}
&&\mu _{0}=\cos t\sinh t+\sin t\cosh t,\qquad \mu _{0}\left( 0\right)
=0,\quad \mu _{0}^{\prime }\left( 0\right) =2,  \label{SolZero} \\
&&\mu _{1}=\cos t\cosh t+\sin t\sinh t,\qquad \mu _{1}\left( 0\right)
=1,\quad \mu _{1}^{\prime }\left( 0\right) =0  \label{SolOne}
\end{eqnarray}
and Green's function is given in terms of trigonometric and hyperbolic
functions as follows
\begin{eqnarray}
G\left( x,y,t\right)  &=&\frac{1}{\sqrt{2\pi i\left( \cos t\sinh t+\sin
t\cosh t\right) }}  \label{modGreen} \\
&&\times \exp \left( \frac{\left( x^{2}-y^{2}\right) \sin t\sinh
t+2xy-\left( x^{2}+y^{2}\right) \cos t\cosh t}{2i\left( \cos t\sinh t+\sin
t\cosh t\right) }\right)   \notag
\end{eqnarray}
as a simple consequence of expressions (\ref{in2})--(\ref{in5}). It is worth
noting that formula (\ref{modGreen}) has been obtained in \cite{Me:Co:Su} by
a totally different approach using the $SU\left( 1,1\right) $-symmetry of $n$
-dimensional oscillator wave functions and properties of the
Meixner--Pollaczek polynomials. More details can be found in\ \cite{Me:Co:Su}
, \cite{Cor-Sot:Sus}.

\section{An Extension}

The degenerate parametric amplification with time-dependent amplitude and
phase had been considered by Raiford \cite{Raiford74}. The corresponding
Hamiltonian, without damping and neglecting high-frequency terms, has the
form
\begin{equation}
H\left( t\right) =\frac{\hslash \omega }{2}\left( \widehat{a}\ \widehat{a}
^{\dagger }+\widehat{a}^{\dagger }\ \widehat{a}\right) -\frac{\hslash
\lambda \left( t\right) }{2}\left( e^{i\left( 2\omega t+\delta \left(
t\right) \right) }\ \widehat{a}^{2}+e^{-i\left( 2\omega t+\delta \left(
t\right) \right) }\left( \widehat{a}^{\dagger }\right) ^{2}\right) .
\label{RaiHam}
\end{equation}
In this model, the phenomenological coupling parameter $\lambda \left(
t\right) ,$ which describes the strength of the interaction between the
quantized signal of frequency $\omega $ and the classical pump of frequency $
2\omega ,$ and the pump phase $\delta \left( t\right) $ are in general
functions of time as indicated. It includes the special case of the pump and
signal being off-resonance by a given amount $\epsilon ,$ i.~e., the pump
frequency being $2\omega +\epsilon ,$ by letting $\delta \left( t\right)
=\epsilon t$ and $\lambda \left( t\right) =\lambda ,$ a constant. More
details, including Heisenberg's equation of motion for the annihilation and
creation operators for the signal photons, can be found in the original
paper \cite{Raiford74}.

In the coordinate representation, the Hamiltonian (\ref{RaiHam}) is
rewritten as
\begin{eqnarray}
&&H=\frac{1}{2m}\left( 1+\frac{\lambda \left( t\right) }{\omega }\cos \left(
2\omega t+\delta \left( t\right) \right) \right) p^{2}\label{RaiHamPX} \\
&&\qquad +\frac{m\omega ^{2}}{2}\left( 1-\frac{\lambda \left( t\right) }{
\omega }\cos \left( 2\omega t+\delta \left( t\right) \right) \right) x^{2}
\notag \\
&&\qquad \quad +\frac{\lambda }{2}\sin \left( 2\omega t+\delta \left(
t\right) \right) \left( px+xp\right) ,\qquad p=\frac{\hslash }{i}\frac{
\partial }{\partial x}  \notag
\end{eqnarray}
and Green's function can be found by using the method described in
section~2. Further details are left to the reader.

\section{A Conclusion}

In this Letter, we have constructed Green's function of the degenerate
parametric oscillator (\ref{at3}) in terms of standard non-periodic
solutions of Ince's equation. To the best of our knowledge, this oscillator
was introduced for the first time by Takahasi \cite{Taka65} in order to
describe the process of \ degenerate parametric amplification in quantum
optics (see also \cite{Louisel64}, \cite{Orszag}, \cite{Per:Hra:Jur}, \cite
{Raiford70}, \cite{Raiford74}, \cite{Stoler74}, \cite{Walls:Milburn} and
references therein). The corresponding Hamiltonian had also been considered
later by Angelow and Trifonov \cite{Angelow:Trifonov95}, \cite{Angelow98} in
order to describe the light propagation in a nonlinear anisotropic
waveguide. Our observation motivates further investigation of properties of
the degenerate parametric oscillator including a systematic study of
corresponding non-periodic solutions of Ince's equation that seems to be
missing in mathematical literature. Moreover, the Hamiltonian (\ref{RaiHam}
), which was considered by Raiford \cite{Raiford74} in the case of the
degenerate parametric amplification with time-dependent amplitude and phase,
requires a generalization of Ince's equation.

\noindent \textbf{Acknowledgments.\/} We thank Professor Carlos Castillo-Ch
\'{a}vez, Professor Vladimir~I.~Man'ko, Professor Svetlana Roudenko and
Professor Dimiter A. Trifonov for support, valuable discussions and
encouragement. The authors are indebted to Professor Victor V.~Dodonov for
kindly pointing out papers \cite{Taka65}, \cite{Raiford74} and \cite
{Stoler74} on the degenerate parametric amplification to our attention ---
his valuable comments have essentially helped to improve the paper. This
work is done as a part of the summer 2010 program on analysis of
Mathematical and Theoretical Biology Institute (MTBI) and Mathematical,
Computational and Modeling Sciences Center (MCMSC) at Arizona State
University. The MTBI/SUMS Summer Research Program is supported by The
National Science Foundation (DMS-0502349), The National Security Agency
(DOD-H982300710096), The Sloan Foundation and Arizona State University. One
of the authors (RCS) is supported by the following National Science
Foundation programs: Louis Stokes Alliances for Minority Participation
(LSAMP): NSF Cooperative Agreement No. HRD-0602425 (WAESO LSAMP Phase IV);
Alliances for Graduate Education and the Professoriate (AGEP): NSF
Cooperative Agreement No. HRD-0450137 (MGE@MSA AGEP Phase II).


TITLE: What does directional derivative zeros imply when directional vector is not zero?
QUESTION [8 upvotes]: This question might sound stupid but I want to confirm an answer from it.
I saw somewhere online that it means that when the directional derivative of function $f$ along the none zero vector $v$ at certain point is equal to $0$, it means that the function $f$ is constant in that direction. But what does "constant in direction" mean? can anyone give me an example of it such as $f(x,y)$ to explain this?
Thanks!

REPLY [1 votes]: Suppose that $f$ is differentiable and that the directional derivative of $f$ along the vector $v=(a_1,a_2,...,a_n)$ is zero i.e. $\nabla f\cdot v=0$ or $a_1f_1+a_2f_2+...+a_nf_n=0$.
Claim: $f$ is constant along any line having direction $v$
$L$ be any such line which passes through $p=(b_1,b_2,...,b_n)$. Now the parametric equations of $L$ are: $x_1(s)=b_1+a_1s,x_2(s)=b_2+a_2s,...,x_n(s)=b_n+a_ns$. Keeping this in mind, we define $F(s)=f(x_1(s),x_2(s),...,x_n(s))$. The derivative $F'(s)=f_1\dot x_1(s)+f_2\dot x_2(s)+...+f_n\dot x_n(s)=f_1a_1+f_2a_2+...+f_na_n=0$. So $F$ must be a constant function. This shows that $f$ is constant along $L$.


TITLE: A variation of conformal vector field
QUESTION [2 upvotes]: Let $H$ be a smooth hypersurface in $\mathbb{R}^d$. I'm interested in vector field $X \in \mathfrak{X}(H)$ for which there is $f \in C^{\infty}(H)$ such that:
\begin{align}
\forall (u,v) \in \mathfrak{X}(H)^2, \quad L_Xg(u,v)=fII(u,v),
\end{align}
where $L$ is the Lie derivative, $g$ the induced metric on $H$ and $II$ it's second fundamental form. In the case of the euclidian sphere, which was my first interest, $II=g$ so that such vector fields are conformal. Is this relation well known in the general case? Can we say anything about $X$?

REPLY [5 votes]: Using that
$$ L_X g = \mathrm{Sym}(\nabla X^\flat) $$
where $\nabla$ is the induced Levi-Civita connection and
$$ II(u,v) = - \langle D_u n,v\rangle_e $$
where $D$ is the Euclidean connection on $\mathbb{R}^d$,
we get that your relation is equivalent to
$$ \mathrm{proj}_H \mathrm{Sym} ( D (X + fn)^{\flat_e}) = 0 $$
where $\flat_e$ is the lowering operation relative to the Euclidean ambient metric, and $\mathrm{proj}_H$ is the projection from the tensor bundle over $\mathbb{R}^d$ to that over $H$. Therefore a sufficient condition for $X$ is that that $X + fn$ is a Euclidean Killing vector field, and since $X$ is by definition tangential to $H$, this means that $X$ is the projection of a Euclidean Killing field to $H$. In the case of the sphere, this recovers the conformal fields of $\mathrm{S}^{d-1}$.
At the moment I am not sure whether there are other solutions, but am tempted to guess that there are not.


TITLE: Why can the deformation of a spring at a given point of the spring be considered directly proportional to the relative distance of the point?
QUESTION [1 upvotes]: Hello i have been studying differential equations and in one example my professor tries to deduce the partial differential equation that describes the longitudinal displacement on a elastic, homogeneous, bar of constant transversal sections. In the text he gives a stiffness coefficient $K$ that apparently works as a type of normalized hook's constant in relation to the length of the bar.
Here begins my actual question: in the deduction of the partial differential equation he uses that if you take a elastic bar in its relaxed position as in the following figure:

And compress it with a force $F$ generating a displacement $S$ as in the figure below:

But then he argues that if we take the middle point of the bar in its relaxed position, that is the point of distance $L/2$, in the compressed state this point "must have" a displacement of "$S/2$" as in the figure bellow:

So my question is why could he make the statement that the deformation is directly proportional to the ratio between the distance of a point and the length of the bar ?, implying that $K=Lk$, where $k$ is the hook's constant, is it always true or only in this specific case ? how can one demonstrate this fact ?

REPLY [2 votes]: The formal way to treat a problem like this is to consider a small slice of the bar at a location $x$. The deformation amount as a function of location is $u(x)$ with the end deformation $\delta=u(\ell)$.
Take a slice of width ${\rm d}x$ and balance the forces, assuming the tension changes by ${\rm d}T$ across the slice. The amount of strain $\epsilon$ of the slice is given by Hook's law $$\left. \epsilon =\frac{\sigma}{E} = \frac{T}{E A} = \frac{u(x+{\rm d}x)-u(x)}{{\rm d}x} = \frac{\partial}{\partial x} u(x) \right\}\;\; T = E A \,\frac{\partial u}{\partial x} $$
Where $A$ is the sectional area and $E$ the Young's Modulus (or Modulus of Elasticity).

Now for the balance of forces:
$$ \left. (T + {\rm d}T) -T = 0 \right\} \;\; {\rm d}T = 0 $$
If ${\rm d}T=0$ then $$\left. T=E A \,\frac{\partial u}{\partial x}=\mbox{(const.)} \right\} \;\;u(x) \mbox{ is linear with }x$$
Given the end conditions in the problem, the solution is
$$ \boxed{ u(x) = \delta \,\frac{x}{\ell} }$$
Or, as you stated, the local deformation is a linear function of position. But for other end conditions, the solution might be different. What to take from this is if there is no external forces applied throughout the bar, the tension is constant, and the deformation is linear with position.


TITLE: The codimension of a parabolic subalgebra of a semisimple Lie algebra
QUESTION [8 upvotes]: Given a complex semisimple Lie algebra $\mathfrak g$ and a subalgebra $\mathfrak h$. If  we are given that the complex vector space $\mathfrak g/\mathfrak h$ has dimension $1$ over $\mathbb C$. Is $\mathfrak h$ a parabolic subalgebra. i.e., contains a Borel subalgebra?
Is the statement above still true when $\dim_{\mathbb C}\mathfrak g/\mathfrak h=2$?

REPLY [4 votes]: The first assertion is true and actually the only possibility (up to direct product of both $\mathfrak{g}$ and $\mathfrak{h}$ by some other semisimple algebra) is when $\mathfrak{g}$ is $\mathfrak{sl}_2$.
(Edit: switched to an algebraic proof)
In $\mathfrak{sl}_2$ it is easy to check that all codimension 1 subalgebras are conjugate to the Borel (parabolic) one.
It is enough to prove that if $\mathbf{g}$ is simple of rank $\ge 2$, then it has no codimension 1 subalgebra $\mathfrak{h}$.
Choose a Cartan subalgebra $\mathfrak{h}_0$ of $\mathfrak{h}$. It induces a grading $(\mathfrak{g}_\alpha)$ of $\mathfrak{g}$, which has to be a quotient of its own Cartan grading.
If $\mathfrak{h}_0=\mathfrak{g}_0$, then $\mathfrak{h}_0$ is a Cartan subalgebra of $\mathfrak{g}$, so $(\mathfrak{g}_\alpha)$ is the Cartan grading of $\mathfrak{g}$. In this case, it follows that $\mathfrak{h}$ is a graded subalgebra containing $\mathfrak{g}_0$, so there exists a nonzero root $\alpha$ such that $\mathfrak{h}=\bigoplus_{\beta\neq\alpha}\mathfrak{g}_\beta$.
Using that $\mathfrak{g}$ has rank $\ge 2$ and is simple, there exist two nonzero roots summing to $\alpha$, and this implies that such $\mathfrak{h}$ is not a subalgebra, contradiction.
Next, if $\mathfrak{h}_0\neq\mathfrak{g}_0$, then being a quotient of the Cartan grading $(\mathfrak{g}_{(\gamma)})$ of $\mathfrak{g}$, we have $\mathfrak{g}_0$ reductive and containing a Cartan subalgebra of $\mathfrak{g}$. If $\mathfrak{g}\neq\mathfrak{g}_0$, we can argue as follows: $\mathfrak{g}_0$ is the sum of $\mathfrak{g}_{(\gamma)}$ where $\gamma$ ranges over some proper subspace $M$ of the space of roots. Since $\mathfrak{g}$ is simple, the set of roots $\gamma$ not in $M$ generates the space of roots (the set of roots is not contained in the union of two proper subspaces), and for each such $\gamma$, we have $\mathfrak{g}_{(\pm\gamma)}\in\mathfrak{h}$ and hence $h_\gamma\in\mathfrak{h}$. Hence $\mathfrak{h}$ contains a Cartan subalgebra of $\mathfrak{g}$, and this implies $\mathfrak{h}_0=\mathfrak{g}_0$ contradiction.
So we have $\mathfrak{g}=\mathfrak{g}_0$: in particular $\mathfrak{h}=\mathfrak{h}_0$. Hence $\mathrm{ad}(h)$ is nilpotent for every $h\in\mathfrak{h}=\mathfrak{h}_0$. Letting $\mathfrak{c}$ be a Cartan subalgebra of $\mathfrak{g}$, this implies that every element of $\mathfrak{h}\cap\mathfrak{c}$ is in the kernel of every root. Since the intersection of kernels of roots is zero in a semisimple Lie algebra, this forces $\mathfrak{c}$ to have dimension $\le 1$. Hence $\mathfrak{g}$ has rank $\le 1$, a contradiction.

Edit: I was a bit frustrated to make such a proof for such a weak result, but indeed it adapts to the following more natural and stronger (and classical) statement:

Let $\mathfrak{g}$ be a absolutely simple Lie algebra over a field of characteristic zero, of (absolute) rank $r$. Then $\mathfrak{g}$ has no proper subalgebra of codimension $<r$.

Lemma: let $\Phi$ be an irreducible root system in dimension $r\ge 1$. Then $\Phi$ is not contained in the union of two proper subspaces.
This follows from:
Sublemma: let $\Phi$ be a root system in dimension $r$ (not necessarily generating). Suppose that $\Phi\subset V_1\cup V_2$ where $V_i$ are subspaces. Then there exist subsets $\Phi_1,\Phi_2$ such that $\Phi_i\subset V_i$, $\Phi_1\cup\Phi_2=\Phi$, and $\langle\Phi_1,\Phi_2\rangle=0$.
Proof of sublemma. This is vacuously true in dimension $0$, and more generally if $V_2=V$. In dimension $r\ge 1$, write $\Psi_1=\Phi\smallsetminus V_2$, $\Psi_2=\Phi\smallsetminus V_1$, and $\Psi_{12}=\Psi\cap V_1\cap V_2$. Clearly $\Psi$ is the disjoint union $\Psi_1\sqcup\Psi_2\sqcup\Psi_{12}$. Also, $\Psi_1,\Psi_2$ are orthogonal: indeed otherwise, we can find $\alpha\in\Psi_1$, $\beta\in\Psi_2$ with $\langle\alpha,\beta\rangle<0$, so $\alpha+\beta\in\Phi$ and this is a contradiction because $\alpha+\beta$ belongs to neither $V_1$ nor $V_2$.
Next, we consider the subspace $V_2$, and its two subspaces $W_1=V_1\cap V_2$, and $W_2$ the orthogonal of $V_1$ in $V_2$, and $\Phi'=\Phi\cap V_2=\Psi_2\sqcup \Psi_{12}$, with $\Psi_2\subset W_2$ and $\Psi_{12}\subset W_1$. We argue by induction inside $V_2$ (the trivial case $V_2=V$ being excluded), to infer that we can write $\Phi\cap V_2=\Phi'_1\cup\Phi'_2$ with $\langle\Phi'_1,\Phi'_2\rangle =0$ and $\Phi'_i\subset W_i$. Then $\Phi=\Psi_1\sqcup\Phi'_1\sqcup\Phi'_2$, with $\Phi'_2\subset V_2$ orthogonal to $\Psi_1\sqcup\Phi'_1\subset V_1$. This finishes the induction.$\Box$
Now let us proceed to the proof of the result. It's an adaptation of the previous proof. Only the first case requires a modification, which is the reason for the above lemma. Namely, let $\mathfrak{h}$ have codimension $<r$ and suppose, with the previous notation that $\mathfrak{h}_0=\mathfrak{g}_0$. In this case the grading is the Cartan grading of $\mathfrak{g}$, so $\mathfrak{h}=\mathfrak{g}_0\oplus\bigoplus_{\alpha\in\Phi\smallsetminus F}\mathfrak{g}_\alpha$, where $F$ is a subset of the root system $\Phi$ of $\mathfrak{g}$, of cardinal $<r$.
Let $V_1$ be the subspace spanned by $F$ (a proper subspace of $\mathfrak{g}_0^*$). Fix a root $\alpha\in F$, and $V_2$ its orthogonal. Then by the lemma, there exists a root $\beta\notin V_1\cup V_2$. Let $P$ be the plane generated by $\alpha$ and $\beta$. So $\Phi\cap P$ is an irreducible root system in $P$, and we can find in $P$ two roots, not collinear to $\alpha$, and avoiding the $V_2\cap P$ (which has dimension $\le 1$), with negative scalar product and summing to $\alpha$. This shows that $\mathfrak{g}_\alpha\subset\mathfrak{h}$, a contradiction.
In the other case $\mathfrak{g}_0\neq\mathfrak{h}_0$, we almost only need to copy the previous proof.


text:
\begin{document}
\begin{frontmatter}
\title{2D force constraints in the method of regularized Stokeslets}

\author[OM]{Ondrej Maxian} \ead{om759@cims.nyu.edu}
\address[OM]{Courant Institute of Mathematical Sciences, New York, NY 10012 U.S.A.}
\author[WS]{Wanda Strychalski\corref{cor}}  \ead{wis6@case.edu}
\cortext[cor]{Corresponding author.}
\address[WS]{Department of Mathematics, Applied Mathematics, and Statistics, Case Western Reserve University, Cleveland, OH  44106, U.S.A.}

\begin{abstract}
\input{abstract.tex}
\end{abstract}

\begin{keyword}
fluid-structure interaction \sep Stokes flow \sep Stokes' paradox \sep regularized Stokeslets

\MSC 65M80 \sep 74F10 \sep 92C37
\end{keyword}
\end{frontmatter}

\section{Introduction}
\input{Intro.tex}

\section{Mathematical framework}
\input{UniformViscosity.tex}

\section{Discretization}
\input{Discretization.tex}

\section{Examples}

\label{sec:unifex}
\input{Examples.tex}

\section{Conclusion}
\input{Conclusion.tex}

\section*{Acknowledgements}
The authors thank Ricardo Cortez, Saverio Spagnolie, Charles Peskin, Aleksander Donev, and Yoichiro Mori for helpful discussions. We thank Alex Mogilner for discussions regarding the cell migration model. OM is supported by the NSF Graduate Research Fellowship DGE-1342536 and the Henry MacCracken fellowship. This work was also supported by a grant from the Simons Foundation (\#429808, Wanda Strychalski).

\section*{References}

\bibliographystyle{model5-names}
\bibliography{StokesBib}

\end{document}

config: arxiv


TITLE: Continuity problem involving piecewise function and glue points
QUESTION [0 upvotes]: Determine all values of the constant a such that the function f(x) below is continuous at the glue point.
$$f(x)= \{ \frac{ax}{tanx}, x>0; a^2-2,x\leqslant0$$
I tried to do many things like setting them equal to each other, and graphing the functions with different values of a

REPLY [0 votes]: For the function to be continuous at $x=0$, you need the left-hand and right-hand limits as $x\rightarrow0$ to be equal. Well...
$$\lim_{x\rightarrow 0^+}f(x)=\lim_{x\rightarrow 0^+}(a^2-2)=a^2-2$$
$$\lim_{x\rightarrow 0^-}f(x)=\lim_{x\rightarrow 0^-}\frac{ax}{\tan x}=\lim_{x\rightarrow 0^-}(a\cdot\frac{x}{\sin x}\cdot\cos x)=a$$
So find the values of $a$ that make $a^2-2=a$.


TITLE: On uniqueness of sums of prime powers
QUESTION [11 upvotes]: An exercise in number theory led to me to the following problem:

Find all solutions $(p,n,q,m)$ of the following equation:
  $$\sum_{k=0}^n p^k = \sum_{h=0}^m q^h,$$
  where $p<q$ are distinct primes, and $1 \le m < n$ are indeterminates.

Numerical evidence gives me the only solution $$(p,n,q,m)=(2,4,5,2).$$
There might be no other solution: I have no idea on how to show it. For those interested of the source of this equation, here it comes the interesting exercise in number theory.

Find all numbers $A$ such that the sum of divisors of $A$ divisible by
  $5$ equals the sum of divisors of $A$ divisible by $2$: $$\sum_{2|d|A}
> d = \sum_{5|e|A} e. $$

Clearly this is equivalent to the condition that the sum of divisors of $A$ NOT divisible by $5$ equals the sum of divisors of $A$ NOT divisible by $2$.
Let's factorize $A= 2^n \cdot 5^m \cdot w$, with $n,m \ge 0$ and $\gcd(w,10)=1$. Without loss fo generality we can consider $n,m\neq 0$. Then we get
$$\sum_{d|2^nw} d = \sum_{e|5^mw} e$$
which is equivalent (using multiplicativity of the sum-of-divisors-function) to
$$\sum_{k=0}^n 2^k = \sum_{h=0}^m 5^h$$
Giving us $A=2^4 \cdot 5^2 \cdot w=400w$ (with $\gcd(w,10)=1$).
Thus, I am looking for simple generalizations with arbirary primes $p \neq q$.

REPLY [4 votes]: Promoting my comment to an answer: the question here is almost exactly the Goormaghtigh Conjecture — the only difference is that the Goormaghtigh Conjecture is just slightly more general and drops the primality condition on the two bases.  Unfortunately, none of the references I'm able to find for it suggests that the special case is any simpler than the general conjecture (it's not at all clear how to make use of the primality of the bases to simplify the problem), and given the sparsity of the results, it seems as though the conjecture isn't attracting too much active research.


text:
chapter \<open>Instantiation for \<open>HOL-ex.Unification\<close> from session \<open>HOL-ex\<close>\<close>

theory Unification_Compat
imports
  "HOL-ex.Unification"
  Term_Class
begin

text \<open>
  The Isabelle library provides a unification algorithm on lambda-free terms. To illustrate
  flexibility of the term algebra, I instantiate my class with that term type. The major issue is
  that those terms are parameterized over the constant and variable type, which cannot easily be
  supported by the classy approach, where those types are fixed to @{typ name}. As a workaround, I
  introduce a class that requires the constant and variable type to be isomorphic to @{typ name}.
\<close>

hide_const (open) Unification.subst

class is_name =
  fixes of_name :: "name \<Rightarrow> 'a"
  assumes bij: "bij of_name"
begin

definition to_name :: "'a \<Rightarrow> name" where
"to_name = inv of_name"

lemma to_of_name[simp]: "to_name (of_name a) = a"
unfolding to_name_def using bij by (metis bij_inv_eq_iff)

lemma of_to_name[simp]: "of_name (to_name a) = a"
unfolding to_name_def using bij by (meson bij_inv_eq_iff)

lemma of_name_inj: "of_name name\<^sub>1 = of_name name\<^sub>2 \<Longrightarrow> name\<^sub>1 = name\<^sub>2"
using bij by (metis to_of_name)

end

instantiation name :: is_name begin

definition of_name_name :: "name \<Rightarrow> name" where
[code_unfold]: "of_name_name x = x"

instance by standard (auto simp: of_name_name_def bij_betw_def inj_on_def)

end

lemma [simp, code_unfold]: "(to_name :: name \<Rightarrow> name) = id"
unfolding to_name_def of_name_name_def by auto

instantiation trm :: (is_name) "pre_term" begin

definition app_trm where
"app_trm = Comb"

definition unapp_trm where
"unapp_trm t = (case t of Comb t u \<Rightarrow> Some (t, u) | _ \<Rightarrow> None)"

definition const_trm where
"const_trm n = trm.Const (of_name n)"

definition unconst_trm where
"unconst_trm t = (case t of trm.Const a \<Rightarrow> Some (to_name a) | _ \<Rightarrow> None)"

definition free_trm where
"free_trm n = Var (of_name n)"

definition unfree_trm where
"unfree_trm t = (case t of Var a \<Rightarrow> Some (to_name a) | _ \<Rightarrow> None)"

primrec consts_trm :: "'a trm \<Rightarrow> name fset" where
"consts_trm (Var _) = {||}" |
"consts_trm (trm.Const c) = {| to_name c |}" |
"consts_trm (M \<cdot> N) = consts_trm M |\<union>| consts_trm N"

context
  includes fset.lifting
begin

lift_definition frees_trm :: "'a trm \<Rightarrow> name fset" is "\<lambda>t. to_name ` vars_of t"
  by auto

end

lemma frees_trm[code, simp]:
  "frees (Var v) = {| to_name v |}"
  "frees (trm.Const c) = {||}"
  "frees (M \<cdot> N) = frees M |\<union>| frees N"
including fset.lifting
by (transfer; auto)+

primrec subst_trm :: "'a trm \<Rightarrow> (name, 'a trm) fmap \<Rightarrow> 'a trm" where
"subst_trm (Var v) env = (case fmlookup env (to_name v) of Some v' \<Rightarrow> v' | _ \<Rightarrow> Var v)" |
"subst_trm (trm.Const c) _ = trm.Const c" |
"subst_trm (M \<cdot> N) env = subst_trm M env \<cdot> subst_trm N env"

instance
by standard
   (auto
      simp: app_trm_def unapp_trm_def const_trm_def unconst_trm_def free_trm_def unfree_trm_def of_name_inj
      split: trm.splits option.splits)

end

instantiation trm :: (is_name) "term" begin

definition abs_pred_trm :: "('a trm \<Rightarrow> bool) \<Rightarrow> 'a trm \<Rightarrow> bool" where
"abs_pred_trm P t \<longleftrightarrow> True"

instance proof (standard, goal_cases)
  case (1 P t)
  then show ?case
    proof (induction t)
      case Var
      then show ?case
        unfolding free_trm_def
        by (metis of_to_name)
    next
      case Const
      then show ?case
        unfolding const_trm_def
        by (metis of_to_name)
    qed (auto simp: app_trm_def)
qed (auto simp: abs_pred_trm_def)

end

lemma assoc_alt_def[simp]:
  "assoc x y t = (case map_of t x of Some y' \<Rightarrow> y' | _ \<Rightarrow> y)"
by (induction t) auto

lemma subst_eq: "Unification.subst t s = subst t (fmap_of_list s)"
by (induction t) (auto split: option.splits simp: fmlookup_of_list)

end

subset_name: curated


TITLE: Non-uniform scaling
QUESTION [0 upvotes]: I have 10 numbers $x_1, x_2, \dots, x_9, x_{10}$ which sum to a total of $1,000$.
I want to scale these numbers so the total is equal to $10,000$, however I don't want them to scale exactly.
I'm not sure exactly (will depend on the output) of the variance required in the scaling but I'm thinking somewhere $\approx10\%$.
Does this question make sense? Can anyone provide a solution on how to achieve this?
Thanks

REPLY [1 votes]: Multiply each number by a randomly chosen real number between $9$ and $11$.
Scale the resulting numbers by the same factor to make their sum equal to $10000$.

If the result must be integer, add step (3) rounding and (4) adjusting some numbers by $\pm 1$ to bring the sum to $10000$ again, as necessary.


TITLE: Proving $\eqalign{n\choose k} = {n\choose k-1} \frac{n-k+1}{k}$
QUESTION [3 upvotes]: I want to prove that
$$\eqalign{n\choose k} = {n\choose k-1} \frac{n-k+1}{k}$$
So I can just write
$$\eqalign{\frac{n!}{(n-k)!\cdot k!}} = {\frac{n!}{(n-(k-1))! \cdot(k-1)!}} \frac{n-k+1}{k}$$
$$\eqalign{\frac{n!}{(n-k)!\cdot k!}} = {\frac{n!}{(n-k+1)! \cdot (k-1)!}} \frac{n-k+1}{k}$$
$$\eqalign{\frac{n!}{(n-k)!\cdot k!}} = {\frac{n!}{(n-k)! \cdot (k-1)!}}$$
I think I did a mistake in the last line, because I cancelled $(n-k+1)!$ since $(n-k+1)! = (n-k+1) \cdot (n-k)!$
So, is my mistake that I cancelled faculty wrongly?

REPLY [2 votes]: $$ \frac {n!}{(n−k)!⋅k!} = \frac {n!(n-k+1)}{(n−k+1)!⋅(k−1)!k} $$
At this point you are done
$$ = \frac {n!}{(n−k)!⋅k!}$$


TITLE: Find the eigenvalues of the matrices.
QUESTION [0 upvotes]: The characteristic equations for the two matrices are:
$x^3-8x-7=0$
and $x^3-6x^2+11x-6=0$
I know that in order to find the eigenvalues, I need to factor these two equations out. I'm just having a brain freeze on how to factor cubic polynomials.
Can anyone refresh my memory on solving these?

REPLY [7 votes]: A very useful tool in factoring $n$ degree polynomials is the rational root theorem. It states that all rational roots of a polynomial are a ratio between individual prime factors of the constant term, and the leading coefficient:
$$x^3-6x^2+11x-6=0$$
$$x=\pm\frac{3,2,1}{1}$$
So assuming there are rational roots, they are either $x=\pm3$, $x=\pm2$, or $x=\pm1$. By trial and error (plugging them in to see if any of them evaluate the polynomial to $0$), you can find that the roots are $x=1,2,3$.


text:
\begin{document}

\def\tr{{\rm Tr}}

\title{Exponential approach to, and properties of, a
    non-equilibrium steady state in a dilute gas}

\author{\vspace{5pt} Eric A. Carlen$^1$, Joel L. Lebowitz$^{1,2}$ and
  Cl\'ement Mouhot$^{3}$ \\
  \vspace{5pt}\small{$1.$ Department of Mathematics, $2.$ Department
    of Physics,}\\
  [-6pt]\small{Rutgers University, 110 Frelinghuysen Road, Piscataway NJ 08854-8019 USA}\\
  \vspace{5pt}\small{$3.$  DPMMS, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK}\\
}

\date{\version}

\maketitle

\medskip

\centerline{\it Dedicated to Errico Presutti, friend and mentor.}

\let\thefootnote\relax\footnote{ \copyright \, 2014 by the
authors. This paper may be reproduced, in its entirety, for
non-commercial purposes.}

\begin{abstract} We investigate a kinetic model of a system in contact
  with several thermal reservoirs at different temperatures
  $T_\alpha$.  Our system is a spatially uniform dilute gas whose
  internal dynamics is described by the nonlinear Boltzmann equation
  with Maxwellian collisions.  Similarly, the interaction with
  reservoir $\alpha$ is represented by a Markovian process that has
  the Maxwellian $M_{T_\alpha}$ as its stationary state. We prove
  existence and uniqueness of a non-equilibrium steady state (NESS)
  and show exponential convergence to this NESS in a metric on
  probability measures introduced into the study of Maxwellian
  collisions by Gabetta, Toscani and Wenberg (GTW). This
  shows that the GTW distance between the current velocity
  distribution to the steady-state velocity distribution is a Lyapunov
  functional for the system.  We also derive expressions for the
  entropy production in the system plus the reservoirs which is always
  positive.
\end{abstract}

\section{Introduction} \label{intro}

The existence, uniqueness and nature of a non-equilibrium steady state
(NESS) of a system in contact with several reservoirs at different
temperatures and/or chemical potentials continues to be one of the
central problems in statistical mechanics, as is the approach to such
a state.  There are only a few models in which the isolated system
evolves according to classical Hamiltonian mechanics or according to
quantum mechanics for which we have even partial answers to these
questions \cite{BLR,D,LLP}.  In addition to the rather unphysical
models corresponding to harmonic crystals and ideal gases, existence
and uniqueness was proven for systems interacting with soft potentials
in contact with thermal walls \cite{KGI}. The resulting NESS is
spatially non-uniform, and we have little information about its
structure. This is true even for cases in which the system is
described mesoscopically by a one particle distribution $f(x,v,t)$, as
in kinetic theory, where correlatons between particles are neglible.

More is possible to prove for NESS of kinetic systems that are
spatially uniform. Such a system, with one reservoir, but acted upon
by an electric field, is investigated in \cite{CELMR1,CELMR2} and will
be discussed later in this paper.

Here we extend this investigation to the case in which the system is
coupled to several thermal reservoirs at different
temperatures. Remarkably, we find, for the first time we believe, a
Lyapunov functional for such systems. This is described in Sections 2
and 3. Then in Section 4 we consider entropy production for such
systems, and in Section 5 we consider the presence of an external
electric field. Finally, in Section 6 we consider the possibility of
deriving such kinetic models from more microscopic descriptions.

\subsection{Description of the basic model}

As already indicated, we deal in this note with a system

described by the one-particle probability density $f(v,t)$.  We are
interested in particular in the case in which the evolution is given
by the non-linear Boltzmann equation with pseudo-Maxellian molecular
collisions, for which the collision kernel is
$$
Q(f,f) := \int_{\R^3}\int_{\mathbb S^2}
b\left(\frac{v-v_*}{|v-v_*|}\cdot \sigma\right) \big[ f(v_*')f(v') -
f(v_*)f(v) \big]\dd v_* \dd \sigma\ .
$$
Here $\dd \sigma$ denotes the
uniform probability measure on the sphere, and
$$
v' := \frac{v+ v_* + |v-v_*|\sigma}{2}\qquad{\rm and}\qquad v_*' :=
\frac{v+ v_* - |v-v_*|\sigma}{2}\ .
$$
We assume \emph{Grad's angular cut-off} with $\int_{\mathbb S^2}
b(u\cdot \sigma)\dd \sigma =1$ for any unit vector $u$, so that
\begin{equation}\label{norm} \frac12 \int_{-1}^1 b(s)\dd s =1\ .
\end{equation}

Then we can separate $Q(f,f)$ into its gain and loss terms
$$
Q(f,f) = Q^+(f,f) - Q^-(f,f)
$$
where
\begin{equation}\label{res0}
Q^+(f,g)(v) := \int_{\R^3}\int_{\mathbb S^2}
b\left(\frac{v-v_*}{|v-v_*|}\cdot \sigma\right) g(v_*')f(v')\dd v_*
\dd \sigma \ , \quad Q^- (f,g) := f \ .
\end{equation}

Observe that the equilibria cancelling this collision operator are
given by the so-called Maxwellian density
\begin{equation*}
  M_{u,T} (v) := \frac{1}{(2 T \pi)^{3/2}} \exp \left( -
    \frac{|v-u|^2}{2T} \right)
\end{equation*}
with bulk velocity $u \in \R^3$, and temperature $T >0$. We denote
simply $M_T = M_{0,T}$ when $u=0$.

The system is spatially homogeneous, and is coupled to several thermal
reservoirs, indexed by $\alpha$, at temperatures $T_\alpha$.  To
describe the interaction with the reservoirs, we include in the
evolution equation a term of the form $Q(f,R)$. The simplest example
would have two reservoirs at temperatures $T_1$ and $T_2$ with the
same coupling, in which case
\begin{equation}\label{res2}
R := \frac12 M_{T_1} + \frac12 M_{T_2} \ .
\end{equation}

However, our methods and results do not depend very much on this
particular form of $R$, and for this reason we leave the distribution
$R$ unspecified in much of our discussion.

It will be convenient to choose the time scale so that the total loss
term coming from both $Q(f,f)$ and $Q(f,R)$ is simply $-f$. We can do
this whatever the relative strength of the two collision mechanisms by
making an appropriate choice of the time scale so that in terms of a
parameter $\gamma\in (0,1)$, the evolution equation can be written as
\begin{equation}\label{res}
\frac{\partial f}{\partial t} = (1-\gamma)
Q(f,f) + \gamma Q(f,R)\
\end{equation} where $R$ be any given probability density on $\R^3$.
Of course, if $R = M_T$ for some $T$, then $M_T$ is the unique steady
state solution of \eqref{res}.  However, if $R$ is given by
\eqref{res2} for $T_1\neq T_2$, then we have no simple expression for
any steady state.

Without loss of generality, we can scale the energy such that
\begin{equation}\label{normal} \int_{\R^3} v R(v)\dd v= 0 \qquad{\rm
and}\qquad \int_{\R^3} |v|^2 R(v)\dd v= 1\ .
\end{equation}

\begin{lm} Let $f_\infty$ be any steady state probability density of
  \eqref{res}. Then, assuming \eqref{normal}, we have
$$\int_{\R^3} v f_\infty (v)\dd v= 0 \qquad{\rm and}
\qquad    \int_{\R^3} |v|^2 f_\infty (v)\dd v= 1\ .$$
\end{lm}

\begin{proof} Let $f$ be a solution of \eqref{res}. For any test
function $\varphi(v)$, and any two probability densities $f$ and $g$,
we have, by a standard computation
$$
\int_{\R^3} Q(f,g)\varphi (v) \dd v =
\int_{\R^3\times \R^3 \times \mathbb S^2}
b(\sigma\cdot k) f(v)g(v_*) \big[ \varphi(v') - \varphi(v) \big]\dd v
\dd v_* \dd \sigma
$$
where $k := |v-v_*|^{-1} (v-v_*)$.

For $\varphi(v) = v$, $\varphi(v') - \varphi(v) = \frac12(v_*-v +
|v_*-v|\sigma)$.  Decomposing $\sigma = (\sigma\cdot k)k +
\sigma^\perp$, we have that
\begin{multline*}
\int_{\R^3\times \R^3 \times \mathbb S^2} b(\sigma\cdot k)
f(v)f(v_*) \big[\varphi(v') - \varphi(v)\big]\dd v \dd v_* \dd \sigma = \\
\left[1 -\frac12 \int_{-1}^1s b(s) \dd s\right] \int_{\R^3\times
\R^3}\left( \frac{v_*-v}{2} \right) f(v)g(v_*)\dd v \dd v_* \ .
\end{multline*}
Therefore,
\begin{eqnarray*}
  \frac{{\rm d}}{{\rm d}t}\left( \int_{\R^3} f(t,v)v\dd
    v \right) &=& \gamma \int_{\R^3} Q( f,R)v\dd v\nonumber \\ &=& -
  \frac\gamma2\left[1 -\frac12 \int_{-1}^1s b(s) \dd s\right]\int_{\R^3}
  f(t,v)v\dd v \ .
\end{eqnarray*}
By \eqref{norm}, $\left[1 -\frac12 \int_{-1}^1s
b(s) \dd s\right] > 0$, and so the first moment relaxes to zero
exponentially fast.

Likewise, for $\varphi(v) = |v|^2$,
$$\varphi(v') - \varphi(v) = \frac{|v_*|^2 - |v|^2}{2} - \sigma\cdot (v+v_*)\ .$$
Decomposing $\sigma$ as before,
\begin{multline*}
\int_{\R^3\times \R^3 \times \mathbb S^2} b(\sigma\cdot k)
f(v)f(v_*) \big[\varphi(v') - \varphi(v)\big]\dd v \dd v_* \dd \sigma = \\
\left[1 -\frac12 \int_{-1}^1s b(s) \dd s\right] \int_{\R^3\times
\R^3}\left( \frac{|v_*|^2-|v|^2}{2}\right) f(v)g(v_*)\dd v \dd v_* \end{multline*}
Therefore,
\begin{eqnarray}
\frac{{\rm d}}{{\rm d}t}\left( \int_{\R^3}
f(t,v)|v|^2\dd v \right) &=& \gamma \int_{\R^3} Q( f,R)v\dd v\nonumber \\ &=&
- \frac\gamma2\left[1 -\frac12 \int_{-1}^1s b(s)\dd s\right]\left(
\int_{\R^3} f(t,v)|v|^2 \dd v -1\right)\ .\nonumber
\end{eqnarray} It follows that $\left( \int_{\R^3} f(t,v)|v|^2 \dd v
-1\right)$ relaxes to zero exponentially fast. In any steady state,
these moments must have the limiting value.
\end{proof}

\subsection{The fixed-point equation}

Because we have fixed the time scale so that the total loss term is
simply $f$, the steady state equation can be written as
$$
f = (1-\gamma)Q^+(f,f) + \gamma Q^+(f,R).
$$
We now follow a method introduced in \cite{CCG} to solve this
equation.

Define the function $\Phi$ from the space of probability densities on
$\R^3$ into itself by
\begin{equation}\label{Phi} \Phi(f) = (1-\gamma)Q^+(f,f) + \gamma
Q^+(f,R)
\end{equation} so that the steady state equation is simply
\begin{equation}\label{fixed} f = \Phi(f)\ .
\end{equation}

We shall show that $\Phi$ is contractive in the
Gabetta-Toscani-Wennberg metric \cite{GTW}, which is the metric
defined as follows: Let $f$ and $g$ be two probability densities on
$\R^3$ with finite second moments such that the first and second
moments are identical. Let $\widehat{f}$ and $\widehat{g}$ denote
their Fourier transforms.  Then
$$
d_{{\rm GTW}}(f,g) := \sup_{\xi\neq 0}\frac{|\widehat f(\xi) -
  \widehat g(\xi)|}{|\xi|^2}\ .
$$

The following contraction lemma gives us existence and uniqueness of
steady states for \eqref{res}.

\begin{lm}
  For all probability densities $f$ and $g$ with the same first and
  second moments as $R$,
$$
\dgtw(\Phi(f),\Phi(g)) \leq \left(1- \gamma \left[ \frac12 -\frac14
    \int_{-1}^1 sb(s)\dd s\right]\right) \dgtw(f,g)\ .
$$
In particular, if $b$ is even,
$$
\dgtw(\Phi(f),\Phi(g)) \leq \left(1 -\frac{\gamma}{2}\right) \dgtw(f,g)\ .
$$
\end{lm}

\begin{proof} Using the Bobylev formula \cite{Boby75,Bobylev},
\begin{equation}\label{bobylev} \widehat{Q^+}(f,g) = \int_{\mathbb S^2}
f(\xi_+)g(\xi_-)b\left( \sigma \cdot \frac{\xi}{|\xi|} \right)\dd \sigma
\end{equation} where
\begin{equation}\label{bobylev2} \xi_\pm = \frac{\xi \pm
|\xi|\sigma}{2}\ .
\end{equation} Note that $|\xi_+|^2+|\xi_-|^2 = |\xi|^2$.

Then we decompose
$$
\widehat{Q^+}(f,f) - \widehat{Q^+}(g,g) = \widehat{Q^+}(f-g,f) +
\widehat{Q^+}(g,f-g)$$
and we deduce
\begin{eqnarray*}
  \dgtw\left(Q^+(f,f), Q^+(g,g)\right) &=& \sup_{\xi\neq 0}
  \frac{\left|\widehat{Q^+}(f-g,f) + \widehat{Q^+}(g,f-g)\right|}{|\xi|^2}\nonumber\\
  &\leq& \sup_{\xi\neq 0} \int_{\mathbb S^2}
  \frac{|\widehat{f}-\widehat{g}|(\xi_+)|\widehat{g}|(\xi_-)}{|\xi|^2}
  b\left (\sigma\cdot \frac{\xi}{|\xi|}\right)\dd \sigma\nonumber\\
  & +& \sup_{\xi\neq 0} \int_{\mathbb S^2}
  \frac{|\widehat{g}|(\xi_+)|\widehat{f}-\widehat{g}|(\xi_-)}{|\xi|^2}b\left(\sigma\cdot
  \frac{\xi}{|\xi|}\right)\dd \sigma \ .
 \end{eqnarray*}
 Next, using the definition of $\dgtw(f,g)$ and the fact that
 $\|\widehat g\|_\infty \le 1$,
 \begin{multline*} \int_{\mathbb S^2}
\frac{|\widehat{f}-\widehat{g}|(\xi_+)|\widehat{g}|(\xi_-)}{|\xi|^2}b\left(\sigma\cdot
\frac{\xi}{|\xi|}\right)\dd \sigma \\ = \int_{\mathbb S^2}
\frac{|\widehat{f}-\widehat{g}|(\xi_+)|\widehat{g}|(\xi_-)}{|\xi_+|^2}
\frac{|\xi_+|^2}{|\xi|^2} b\left(\sigma\cdot \frac{\xi}{|\xi|}\right)\dd
\sigma \\ \leq  \dgtw(f,g) \int_{\mathbb S^2}
\frac{|\xi_+|^2}{|\xi|^2} b\left(\sigma\cdot \frac{\xi}{|\xi|}\right)\dd \sigma\
.
 \end{multline*}
Likewise,
 $$
 \sup_{\xi\neq 0} \int_{\mathbb S^2}
 \frac{|\widehat{g}|(\xi_+)|\widehat{f}-\widehat{g}|(\xi_-)}{|\xi|^2} b\left(\sigma\cdot
 \frac{\xi}{|\xi|}\right)\dd \sigma \leq \dgtw(f,g) \int_{\mathbb S^2}
 \frac{|\xi_-|^2}{|\xi|^2} b\left(\sigma\cdot
   \frac{\xi}{|\xi|}\right)\dd \sigma\ .
$$
Then since $|\xi_+|^2+|\xi_-|^2 = |\xi|^2$, we have that
  $$
\dgtw\left(Q^+(f,f),Q^+(g,g)\right) \leq \dgtw(f,g)\ .
$$

  Next, by essentially the same calculation,
 \begin{eqnarray*}
\dgtw\left(Q^+(f,R), Q^+(g,R)\right) &\leq& \int_{\mathbb S^2}
\frac{|\widehat{f}-\widehat{g}|(\xi_+)|\widehat{R}|(\xi_-)}{|\xi_+|^2}
\frac{|\xi_+|^2}{|\xi|^2} b\left(\sigma\cdot \frac{\xi}{|\xi|}\right)\dd \sigma
\nonumber\\ &\leq& \dgtw(f,g) \int_{\mathbb S^2} \frac{|\xi_+|^2}{|\xi|^2}
b\left(\sigma\cdot \frac{\xi}{|\xi|}\right)\dd \sigma \ .
 \end{eqnarray*}
 Since $|\xi_+|^2 = \tfrac12(|\xi|^2 + |\xi|(\xi\cdot \sigma))$,
 $$
 \int_{\mathbb S^2} \frac{|\xi_+|^2}{|\xi|^2} b\left(\sigma\cdot
   \frac{\xi}{|\xi|}\right)\dd \sigma = \frac14 \int_{-1}^1 \left(1 +
   s\right)b(s)\dd s \ .
$$
Altogether, by the triangle inequality, we have
$$
\dgtw\left(\Phi(f,f), \Phi(g,g)\right) \leq \left((1-\gamma) + \gamma
  \left[\frac14 \int_{-1}^1 \left(1 + s\right)b(s)\dd s\right]\right)
\dgtw(f,g)\ ,
$$
which gives the result.

\end{proof}

\begin{thm} Suppose that $R$ satisfies \eqref{normal}. Then there is a
unique steady state solution $f_\infty$ of \eqref{res}. Moreover, if
we define a sequence by $f_0 = R$ and $f_{n} = \Phi(f_{n-1})$ for all
$n\in \N$, then $f_\infty = \lim_{n\to\infty} f_n$, and
\begin{equation}\label{iter} \dgtw(f_n,f_\infty) \leq
\frac{\lambda^n}{1-\lambda} \dgtw(\Phi(R),R)\
\end{equation} where
$$\lambda = \left((1-\gamma) + \gamma \left[\frac14 \int_{-1}^1 \left(1 + s\right)b(s)\dd s\right]\right)  < 1\ .$$
\end{thm}

\begin{proof} This is a direct consequence of the previous lemma and
the contraction mapping theorem. Recall from the proof that
$$\dgtw(f_{n+1},f_{n}) \leq \lambda^n \dgtw(\Phi(R),R)$$
Then by the triangle inequality we obtain the final estimate.

\end{proof}

We remark that the equation \eqref{iter} allows the effective
computation of $f_\infty$. Once can certainly evaluate it numerically,
especially when $b$ is even so that $\lambda = 1-\gamma/2$.

\section{Exponential convergence}

Since for small $h$, and any two solutions $f$ and $g$ of \eqref{res},
\begin{equation*}
  f(t+h) - g(t+h) = h[\Phi(f(t)) - \Phi(g(t)] +
  (1-h)[f(t)-g(t)] + o(h)\ ,
\end{equation*} it follows from our contraction estimate in the
previous section that provided $f(0)$ and $g(0)$ have the same first
and second moments as $R$,
\begin{multline*}
\frac{{\rm d}}{{\rm d}t} \dgtw(f(t),g(t)) \\ \le -
\left[1 - \left((1-\gamma) + \gamma \left[\frac14 \int_{-1}^1 \left(1
        + s\right)b(s)\dd s\right]\right)\right] \dgtw(f(t),g(t))\ .
\end{multline*}

 In particular, taking $g(t) = f_\infty$, we see that
$\dgtw(f(t),f_\infty)$ decreases to zero exponentially fast.

We can dispense with the requirement that the initial data $f(0)$ has
the same first and second moments as $R$ by using the correction
technique introduced in \cite{CGT}. Let us describe briefly this
argument. Let us denote
\begin{align*}
  \lambda_0 &:= \frac12 \left[ 1 - \frac12 \int_{-1} ^{+1} s b(s) \dd s
  \right] >0, \\
  \lambda_1 &:= \left[1 - \left((1-\gamma) + \gamma \left[\frac14 \int_{-1}^1 \left(1
          + s\right)b(s)\dd s\right]\right)\right] >0 \ .
\end{align*}
We define (in Fourier variables)
\begin{equation*}
  \widehat{\mathcal{M}}[f] := \chi(\xi) \sum_{|\alpha| \le 2} \left( \int_{\R^3}
    v^\alpha f(v) \dd v \right) \frac{\xi^\alpha}{\alpha !}
\end{equation*}
where we use the standard notation for multi-indeces $\alpha =
(\alpha_1,\alpha_2,\alpha_3) \in \mathbb N^3$, $|\alpha| = \alpha_1 +
\alpha_2 + \alpha_3$, $v^\alpha = v_1^{\alpha_1} v_2^{\alpha_2}
v_3^{\alpha_3}$, $\alpha ! = \alpha_1 ! \alpha_2 ! \alpha_3 !$, and
where $\chi$ is a compactly supported smooth function that is equal to
one around $\xi=0$. Then if we consider two solutions $f$ and $g$ with
possibly different momentum and energy, we write $D = f- g -
\mathcal M[f-g]$, $S = f+g$,
and obtain
\begin{equation*}
  \partial_t D = (1-\gamma) Q (D,S) + (1-\gamma) Q(S,D) + \gamma Q
  (D,R) - W
\end{equation*}
with
\begin{multline*}
  W:= \Big[ \partial_t
    \mathcal M[f-g] + (1-\gamma) Q (\mathcal M[f-g],S) \\ + (1-\gamma)
    Q(S,\mathcal M[f-g]) +
    \gamma Q (\mathcal M[f-g],R) \Big] \ ,
\end{multline*}
and one checks by similar moment estimates as above that
\begin{equation*}
  \left| \widehat W (\xi,t)\right|
  \le C |\xi|^2 \left( \sum_{|\alpha| \le 2} \left| \int_{\R^3}
      v^\alpha \left[f(v,t) - g(v,t) \right]\dd v \right| \right) \le C' |\xi|^2 e^{-\lambda_0
    t}
\end{equation*}
for some constants $C,C' >0$. We then perform the same contraction
estimate as before since $D$ is now the Fourier transform of a
centered zero-energy function, and obtain
\begin{equation*}
  \sup_{\xi \in \R^3} \frac{|\widehat D(\xi,t)|}{|\xi|^2} \le \sup_{\xi \in
    \R^3} \frac{|\widehat D(\xi,0)|}{|\xi|^2} e^{-\lambda_1 t} +
  C'' e^{-\min(\lambda_0,\lambda_1) t}
\end{equation*}
for some constant $C_1$. Finally we deduce by taking $g=f_\infty$ that
$f$ is converging to the equilibrium $f_\infty$ with exponential
rate, measured in the distance
\begin{equation*}
\dgtw'(f,g) = \sup_{\xi \in \R^3} \frac{|\widehat f(\xi)-\widehat g(\xi)- \hat{\mathcal
  M}[f-g](\xi)|}{|\xi|^2} + \left|\mathcal M[f-g]\right| \ ,
\end{equation*}
which writes $\dgtw'(f,f_\infty) \le C''' e^{-
  \min(\lambda_0,\lambda_1)t}$ for some constant $C'''>0$.

\section{Diffusive thermal reservoirs}

In some physical situations it is more appropriate to model the
interaction of a system with reservoirs by an Ornstein-Uhlenbeck
continuous time process rather than the discrete time collision model
that we have considered above.  This leads to a kinetic equation of
the form
\begin{equation}\label{OUR} \frac{\partial}{\partial t}f(v,t) = Q(f,f)
+ \sum_\alpha \eta_\alpha T_\alpha \frac{\partial}{\partial v}
\left[M_{T_\alpha} \frac{\partial}{\partial
v}\left(\frac{f}{M_{T_\alpha}}\right)\right]\
\end{equation} The constant $\eta_\alpha$ sets the strength of the
interaction with the $\alpha$th reservoir.

Note that in this setting, the evolution equation for several
reservoirs reduces to the evolution equation for a single reservoir,
since
$$
\sum_\alpha \eta_\alpha T_\alpha \frac{\partial}{\partial v}
\left[M_{T_\alpha} \frac{\partial}{\partial
v}\left(\frac{f}{M_{T_\alpha}}\right)\right] = \eta T
\frac{\partial}{\partial v} \left[M_{T} \frac{\partial}{\partial
v}\left(\frac{f}{M_{T}}\right)\right]\ $$ where
  $$\eta = \sum_\alpha \eta_\alpha \qquad{\rm and}\qquad T = \frac{1}{\eta} \sum_\alpha \eta_\alpha T_\alpha\ .$$

The effective evolution equation
\begin{equation} \label{effective}
\frac{\partial f}{\partial t} = Q(f,f) + \eta T
\frac{\partial}{\partial v} \left[M_{T} \frac{\partial}{\partial
v}\left(\frac{f}{M_{T}}\right)\right]
\end{equation}
is then easy to analyze since in this case the unique stationary state
is $M_T$, and the relative entropy with respect to $M_T$, i.e.,
$\int_{\R^3} f [ \log f - \log M_T]\dd v$ decreases to zero
exponentially fast.  Likewise,
$$\left[ 3T -  \int_{\R^3}  v^2 f(v,t)\dd v\right] = e^{-2\eta t} \left[ 3T -  \int_{\R^3}  v^2 f(v,0)\dd v\right]\ .$$

In fact, more is true: This evolution also has the contractive
property proved in the previous section for \eqref{res}, now using a
different metric, but one that is equivalent to the GTW metric
\cite{GTW}, namely the $2$-Wasserstein metric.

A theorem of Tanaka \cite{T1,T2} says that the evolution described by
the spatially homogeneous Boltzmann equation for Maxwellian molecules
is contractive in this metric. (We do not describe this metric here,
other than to say that like the GTW metric, it metrizes the topology
of weak convergence of probability measures together with convergence
of second moments, and we refer to the book Villani \cite{V} for
the definition and the proof of this fact.)

As Otto has shown \cite{O}, the evolution described by
\begin{equation*}
  \frac{\partial f}{\partial t} = \eta T
  \frac{\partial}{\partial v} \left[M_{T} \frac{\partial}{\partial
      v}\left(\frac{f}{M_{T}}\right)\right]
\end{equation*}
 is exponentially contractive in the $2$-Wasserstein
metric: If $f$ and $g$ are any two solutions of this equation
$$d_{W_2}(f(\cdot,t),g(\cdot,t)) \leq e^{-\eta (t-s)} d_{W_2}(f(\cdot,s),g(\cdot,s))\ .$$
Together with Tanaka's Theorem for the equation with the collision
operator $Q$ only and a splitting argument (i.e., a non-linear Trotter
product argument), one easily establishes that this same estimate is
valid for solutions of \eqref{effective} and therefore \eqref{OUR}. In
particular, it follows that for all solutions $f(v,t)$ of \eqref{OUR},
\begin{multline} d_{W_2}(f(\cdot,t),M_T) \leq e^{-\eta t}
d_{W_2}(f(\cdot,0),M_T) \leq\\ e^{-(n/T)t} \left( \tfrac12 \int
v^2f(v,0)\dd v + T\right)\ .\end{multline} where we have used a simple
estimate (see \cite{V}) for $d_{W_2}(f(\cdot,0),M_T)$ in the last
inequality.

In conclusion, for the diffusive reservoirs, we have not only a ``free
energy type'' Lyapunov functional, namely the relative entropy with
respect to $M_T$, but also a a different one that is similar in nature
to the one we found in the previous section for \eqref{res}.

\section{Entropy production for thermal reservoirs}

We now consider the entropy production when the reservoirs are
thermal. More precisely, we assume that the time evolution of $f(v,t)$
is given by
\begin{equation}\label{thermal}
  \frac{\partial f}{\partial t} =
  Q(f,f) + \sum_\alpha \int_{\R^3} [K_\alpha(v,v')f(v') -
  K_\alpha(v',v)f(v)]\dd v' \ .
\end{equation}
The Markovian rates $K_\alpha(v,v')$ describe collisions with the
thermal reservoir at temperature $T_\alpha = \beta_\alpha^{-1}$
resulting in a transition from $v'$ to $v$, and $Q(f,f)$ is a general
Boltzmann type collision operator, not necessarily of the Maxwellian
type. We assume detailed balance for each reservoir; i.e.,
\begin{equation}\label{detailed}
\forall \, \alpha, \quad K(v,v')M_{T_\alpha}(v') =
K(v',v)M_{T_\alpha}(v) \ .
\end{equation}

Equations \eqref{thermal} and \eqref{detailed} thus include
\eqref{res} and \eqref{res2} as special cases. However, the
exponential approach proved for the latter may not hold in this more
general case. In fact, even the existence and uniqueness of a
stationary state for \eqref{thermal} and \eqref{detailed} is not
guaranteed; see e.g. \cite{CELMR1,CELMR2}. On the other hand, the
analysis below applies to the broader class of models for which there
is existence and uniqueness of the NESS, and also carries over to the
case in which $f$ and the $K_\alpha$ depend on position $x\in \R^3$,
though we shall not pursue this here.

The rate of change of the system's Boltzmann gas entropy is given by
\begin{equation}\label{enrate}
\dot S = -\frac{{\rm d}}{{\rm d}t} \int_{\R^3}
f \log f \dd v = \sigma_B + \sum_\alpha \sigma_\alpha - \sigma_R
\end{equation}
where $\sigma_B$ is the usual rate of change of the
entropy due to the Boltzmann collision term, which is non-negative and
equal to zero if and only if $f$ is a Maxwellian.  The $\sigma_\alpha$
are given by
\begin{equation}\label{siga}
  \sigma_\alpha := \frac12\int_{\R^3} \int_{\R^3}
  K_\alpha(v,v')M_\alpha(v')[\nu_\alpha(v,t) - \nu_\alpha(v',t)] \log
  \frac{\nu_\alpha(v,t)}{\nu_\alpha(v',t)} \dd v\dd v' \geq 0
\end{equation} where
$$\nu_\alpha(v,t) = \frac{f(v,t)}{M_\alpha(v)}\ .$$
Finally, $\sigma_R$ is the rate of production of entropy in the
reservoirs
\begin{equation}\label{resrate}
  \sigma_R := \sum_\alpha \beta_\alpha J_\alpha
\end{equation} with
\begin{equation}\label{resrate2}
  J_\alpha(t) := \frac12 \int \int
  K_\alpha(v,v')f(v',t)\big[v'^2 - v^2\big]\dd v \dd v'\
\end{equation}
being the flux of energy into the $\alpha$th reservoir at time $t$.

The total rate of entropy production in the system plus reservoirs is
given by
\begin{equation}\label{totalrate} \sigma = \dot S + \sum_\alpha
\beta_\alpha J_\alpha \geq 0\ .
\end{equation} In the stationary state, $\dot S = 0$, and
\begin{equation}\label{stat}\overline{\sigma} = \sum_\alpha
\beta_\alpha \overline{J_\alpha} = \overline{\sigma}_B +
\sum_\alpha\overline{\sigma}_\alpha\geq 0
\end{equation} where the bars denote quantities computed in the
stationary state $f_\infty$.

There is equality in \eqref{stat} if and only if $\overline{\sigma}_B
= 0$ and $\overline{\sigma}_\alpha = 0$ for each $\alpha$, and this is
the case if and only if the stationary state is a Maxwellian, and all
of the reservoirs have the same temperature.

In the case of equal temperature $\beta_\alpha = \beta$ for all
$\alpha$, $\sigma$ in \eqref{totalrate} is given by
$$
\sigma = \frac{{\rm d}}{{\rm d}t}[ S - \beta \mathcal{E}] = - \beta
\frac{{\rm d}}{{\rm d}t} \mathcal{F}$$ where $ \mathcal{E} = \tfrac12
\langle v^2\rangle$ is the average energy of the system and
$\mathcal{F}$ is a kind of free energy.  $ \mathcal{F}$ is thus a
Lyapunov functional achieving its minimum when $f = M_T$. In fact, it
is just the relative entropy of $f$ with respect to $M_T$.

When the $\beta_\alpha$ are unequal, $\sigma$ is not a time
derivative, and the only Lyapunov functionals we know are $d_{\rm
  GTW}(f,f_\infty)$ (or $d_{W_2}(f,f_\infty)$ in the case of diffusive
reservoirs), and only in the case in which the time evolution is given
by \eqref{res}.

In the stationary state we must of course have $\sum_\alpha
\overline{J}_\alpha =0$. Hence if we have only two reservoirs, then
\begin{equation}\label{hotcold} \overline{\sigma} =
(\beta_1-\beta_2)\overline{J}_1 \geq 0\ ;
\end{equation} i.e., heat flows from the hot to the cold reservoir.

We note that when the system is coupled to reservoirs (with equal or
unequal temperatures), then $S$ need not be monotone non-decreasing as
is evident from the fact that we can start from an initial state with
an entropy that is higher than that of the stationary sate; e.g., a
Maxwellian with a sufficiently high temperature. It is only when
$\sigma_R \leq 0$ that $\dot S$ must be nonnegative.

The above considerations remain valid when the Boltzmann collision
kernel $Q(f,f)$ is replaced by the modified Enskog collision kernel
which is generally considered to be a good approximation for a
moderately dense gas; see \cite{GL} and references provided there.

As noted above, in some physical situations it is more appropriate to
model the interaction of a system with reservoirs by an
Ornstein-Uhlenbeck continuous time diffusion process rather than a
discrete time jump process as in \eqref{thermal}.  Our analysis of
entropy production allows for a more general class of diffusive
reservoirs than did our discussion of the contraction property. In
particular, we can allow velocity dependent diffusion coefficients.
Similarly, the collision term may be modified when the system is not a
dilute gas. We shall therefore write a general kinetic equation in the
form
\begin{equation}\label{ias}
\frac{\partial f}{\partial t} = Q(f) +
\sum_\alpha \frac{\partial}{\partial v} \left[\eta_\alpha(v)T_\alpha
M_{T_\alpha} \frac{\partial}{\partial
v}\left(\frac{f}{M_{T_\alpha}}\right)\right]\
\end{equation}
requiring only that $Q(f)$ conserve the energy, momentum and mass of
$f$, and that
$$
-\int_{\R^3} Q(f) \log f \dd v \geq 0
$$
with equality if and only if $f = M_T$ for some $T$. The above
analysis can now be repeated with $\sigma_\alpha$ replaced by
$$
\sigma_\alpha ' := \int_{\R^3} f \eta_\alpha(v)T_\alpha \left|
  \frac{\partial}{\partial
    v}\left(\frac{f}{M_{T_\alpha}}\right)\right|^2 \dd v\ .
$$

\if false In this case the stationary state is given by a Maxwellian
with a reciprocal temperature $\overline{\beta}$ that is a weighted
average of the reciprocal temperatures $\beta_\alpha$ of the
reservoirs.  \fi

 \section{Systems driven by an electric field}

 In the systems considered so far in this paper, a non-equilibrium
 steady state has been maintained by at least two reservoirs with
 energy flowing out of some and into others.  A different sort of
 model is investigated in \cite{CELMR1,CELMR2} which concerns a weakly
 ionized plasma with energy supplied by an electric field $E$, and
 removed into a reservoir by a damping mechanism. In this model, the
 probability density $f(v,t)$ evolves according to
\begin{equation}\label{plasma}
\frac{\partial f}{\partial t} =
\frac{1}{\epsilon}Q(f) -E\cdot \frac{\partial }{\partial v} f + \nu \left[
\tilde f - f\right] + \frac{\partial }{\partial v}\left[ D(v) M
\frac{\partial }{\partial v} \left(\frac{f}{M}\right)\right]\ .
\end{equation}
In \eqref{plasma}, $E$ is a constant electric field, $\epsilon$ is a
small parameter setting the rate of internal collisions, $\tilde{f}$
is the spherical average of $f$ (hence radial, but with the same
energy distribution as $f$), $\nu$ is a constant, and $M$ is a
centered Maxwellian with $T=1$.  See (2.1-4) in \cite{CELMR1}.

When $E= 0$, the unique equilibrium is $M$, and the free energy
$\mathcal{F}$ is a Lyapunov function.  Also when $D(v) = D$
independent of $v$, and $\nu = 0$, the NESS is given by the shifted
Maxwellian $M(v-u)$ where $u = ET/D$. In this case
$$\mathcal{F} = \frac12 \langle |v-u|^2\rangle -\beta S$$
serves as a Lyapunov function governing convergence to the NESS.

However, when $D(v)$ has certain properties, $E\neq 0$ and $\nu >0$,
it is shown in \cite{CELMR1,CELMR2} that for all $\epsilon$
sufficiently small, there is a range of the parameters for which there
are multiple stable stationary solutions of \eqref{plasma}. This means
in particular that there does not exist any global Lyapunov function
for \eqref{plasma} for general parameters.

Finally we have in this stationary state that
$$
\overline{\sigma} = \beta \overline{j}\cdot E \quad \mbox{ where}
\quad \overline{j} = \int_{\R^3} v f_\infty(v)\dd v
$$
is the steady state current.

\section{Microscopic models}

We have investigated here the approach to the NESS and some properties
of that state for models in contact with several thermal reservoirs at
different temperatures. We showed that there exist in some cases
Lyapunov functionals even when we do not know the NESS in explicit
form. It would be nice to show that these models correspond to
$N$-particle microscopic models in some suitable scaling limits.

A plausible conjecture is the following: Suppose that the dynamics of
the isolated system is given by a Hamiltonian microscopic dynamics,
and it yields, in some suitable limit, a Boltzmann equation for the
single particle distribution; e.g., hard-sphere collisions under the
Boltzmann-Grad scaling limit \cite{L,PSS}.  Adding stochastic
interactions with thermal reservoirs should then lead, in a suitable
limit, to \eqref{thermal}.

However, this is beyond our current reach, even for the case in which
the system is in contact with a single reservoir.  If one drops the
requirement that the microscopic dynamics be Hamiltonian, the
situation is much better. Such a microscopic derivation was proven
recently by Bonetto, Loss and Vaidyanathan \cite{BLV} when the
isolated system dynamics is given by the Kac stochastic collision
model and there is a single thermal reservoir. One may expect a
similar result to be valid for the Kac system in contact with several
reservoirs.

\medskip

\noindent{\bf Acknowledgements} Work of E.A.C. is partially supported
by N.S.F. grant DMS 1201354.  J.L.L. wishes to thank the I.A.S. its
hospitality during the course of the work, and his work is partially
supported by N.S.F. grant DMR 1104500 and AFOSR grant FA9550.
C.M. wishes to thank the IAS for its support during his visit in may
2014 when this work was done. His research is also partially funded by
the ERC Starting Grant MATKIT.


text:
\begin{document}

\title{A Benson-Type Algorithm for Bounded Convex Vector Optimization Problems with Vertex Selection}

\author{ \name{Daniel Dörfler\textsuperscript{a}\thanks{CONTACT
D. Dörfler. \Letter \ \href{}{daniel.doerfler@uni-jena.de}, \protect
\includegraphics[keepaspectratio, width=0.8em]{Figures/ORCID.png}
\href{}{https://orcid.org/0000-0002-9503-3619}}, Andreas
Löhne\textsuperscript{a}, Christopher Schneider\textsuperscript{b} and
Benjamin Weißing\textsuperscript{a}} \affil{\textsuperscript{a}
Friedrich Schiller University Jena, Germany; \textsuperscript{b} Ernst
Abbe University of Applied Sciences Jena, Germany} }

\maketitle

\begin{abstract}
  We present an algorithm for approximately solving bounded convex
  vector optimization problems. The algorithm provides both an outer
  and an inner polyhedral approximation of the upper image. It is a
  modification of the primal algorithm presented by Löhne, Rudloff,
  and Ulus in 2014. There, vertices of an already known outer
  approximation are successively cut off to improve the approximation
  error. We propose a new and efficient selection rule for deciding
  which vertex to cut off. Numerical examples are provided which
  illustrate that this method may solve fewer scalar problems overall
  and therefore may be faster while achieving the same approximation
  quality.
\end{abstract}

\begin{keywords}
  Vector optimization, multiple objective optimization, polyhedral
  approximation, convex programming, algorithms
\end{keywords}

\begin{amscode}
  90C29, 90C25, 90-08, 90C59
\end{amscode}

\section{Introduction}
\input{intro}

\section{Preliminaries}\label{sec:2}
\input{prelim}

\section{Vector Convex Programs}\label{sec:3}
\input{vcp}

\section{An Algorithm for Bounded VCPs with Vertex
  Selection}\label{sec:4}
\input{algrthm}

\section{Numerical Examples}\label{sec:5}
\input{nmrcl}

\section{Conclusion}
\input{cncl}

\bibliographystyle{bib/tfs}
\bibliography{bib/references}

\end{document}


text:
\section{Introduction}\label{sec:intro}
We consider the system of stiff (nonlinear) differential equations,
\begin{equation}\label{eq:ode}
y' =f(y)\vcentcolon = \ff(y)+\fs(y), \qquad\qquad
y(0)=y_0,
\end{equation}
where $f:\Rn\rightarrow\Rn$ splits into an expensive but only mildly stiff part, $\fs$, associated with relatively slow (S) time-scales, and a cheap but severely stiff part, $\ff$, associated with fast (F) time-scales.
Typical applications include chemical reactions and electrical circuits with disparate time-scales, but also
spatial discretizations of diffusion dominated (parabolic) partial differential equations (PDEs) with local mesh refinement. Semi-discrete parabolic PDEs, in particular, lead to large systems of stiff ordinary differential equations,
where the eigenvalues of the Jacobian matrix, $\partial f / \partial y$, lie in a narrow strip along the negative real axis whose extent scales as $H^{-2}$ for a mesh size $H$. In the presence of local mesh refinement,
$\fs,\ff$ in \cref{eq:ode} then correspond to discrete diffusion operators in the coarse and locally refined region of the mesh, respectively. Although $\ff$ involves only a small number of degrees of freedom,
the extreme eigenvalue of its Jacobian will determine the spectral radius $\rho$ of $\partial f / \partial y$.

In contrast to multiscale methods \cite{AbP12,E03,EnT05,Van03},
{\it we do not assume any scale separation} in \cref{eq:ode}; hence, $\ff$ may contain {\it both fast and slow} scales.
In a situation of local mesh refinement, for instance, when $\fs$ and $\ff$ each represent the discrete Laplacian in the coarse and refined regions, both Jacobians in fact contain small eigenvalues in magnitude; hence, the spectrum of $\partial f / \partial y$ cannot simply be split into fast and slow modes, as in \cref{fig:noscalesep}. This stands in sharp contrast to the underlying assumption of recently introduced multiscale methods for stiff (dissipative) ODEs, such as the  heterogeneous multiscale method (HMM)
\cite{E03,EnT05} or  the projective method \cite{GIK03} which all require scale separation, as in \cref{fig:scalesep}.

\begin{figure}

		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\coordinate (a) at (0,0);
			\coordinate (b) at (6,0);
			\coordinate (c) at (6,0.5);
			\coordinate (d) at (6,-0.5);
			\coordinate (e) at (0,-1.3);
			\coordinate (f) at (6,-1.3);
			\coordinate (g) at (6,-0.8);
			\coordinate (h) at (6,-1.8);
			\coordinate (i) at (5,-1.3);
			\coordinate (l) at (3,0);
			\draw[thick,<-] (a)--(b);
			\draw[thick,<-] (c)--(d);
			\draw[thick,<-] (e)--(f);
			\draw[thick,<-] (g)--(h);
			\draw[\colorthree,fill=\colorthree,rounded corners=1mm,opacity=0.3] ([xshift=6mm,yshift=-3mm]a) rectangle ([yshift=3mm]b);
			\draw[\colorone,fill=\colorone,rounded corners=1mm,opacity=0.3] ([xshift=0mm,yshift=-3mm]i) rectangle ([yshift=3mm]f);
			\node[fill=white,rounded corners=2pt,inner sep=2pt] at ([xshift=0mm]l) {{$\ff$}};
			\node[fill=white,rounded corners=2pt,inner sep=2pt] at ([xshift=5.1mm]i) {{$\fs$}};
			\end{tikzpicture}
			\caption{Spectrum without scale separation.}
			\label{fig:noscalesep}
		\end{subfigure}\hfill
		\begin{subfigure}[t]{\subfigsize\textwidth}
			\centering
			\begin{tikzpicture}[scale=\plotimscale]
			\coordinate (a) at (0,0);
			\coordinate (b) at (6,0);
			\coordinate (c) at (6,0.5);
			\coordinate (d) at (6,-0.5);
			\coordinate (e) at (0,-1.3);
			\coordinate (f) at (6,-1.3);
			\coordinate (g) at (6,-0.8);
			\coordinate (h) at (6,-1.8);
			\coordinate (i) at (5,-1.3);
			\coordinate (l) at (3,0);
			\draw[thick,<-] (a)--(b);
			\draw[thick,<-] (c)--(d);
			\draw[thick,<-] (e)--(f);
			\draw[thick,<-] (g)--(h);
			\draw[\colorthree,fill=\colorthree,rounded corners=1mm,opacity=0.3] ([xshift=6mm,yshift=-3mm]a) rectangle ([xshift=-14mm,yshift=3mm]l);
			\draw[\colorone,fill=\colorone,rounded corners=1mm,opacity=0.3] ([xshift=0mm,yshift=-3mm]i) rectangle ([yshift=3mm]f);
			\node[fill=white,rounded corners=2pt,inner sep=2pt] at ([xshift=-18.9mm]l) {{$\ff$}};
			\node[fill=white,rounded corners=2pt,inner sep=2pt] at ([xshift=5.1mm]i) {{$\fs$}};
			\end{tikzpicture}
			\caption{Spectrum with scale separation.}
			\label{fig:scalesep}
		\end{subfigure}

	\caption{Stiff problems with identical spectral radii but different eigenvalues distribution.}
	\label{fig:spectrum}
\end{figure}

Standard explicit methods are notoriously inefficient for stiff differential equations due to their stringent stability constraint on the step size, $\tau$, which for parabolic problems must be proportional to $H^{2}$. Implicit methods, on the other hand, are unconditionally stable but require at every time step the solution of an $n \times n$ linear (or possibly nonlinear) system of equations, a high price to pay when $n$ is large.
Moreover, when sheer size calls for using iterative methods, the overall performance heavily relies on the availability of efficient preconditioners while the convergence of Newton-like nonlinear iterations is not even guaranteed for larger step sizes.

{\it Stabilized Runge--Kutta (RK) (or Chebyshev) methods} fall somewhere between explicit and implicit methods: they are explicit and thus avoid the solution of large systems of equations, while their stability interval on the negative real axis is proportional to $s^2$ for an $s$-stage method. Thanks to this remarkable quadratic dependency, the work load (number of stages $s$) per time step only needs to scale linearly with $H^{-1}$ for parabolic PDEs, in contrast to the quadratic increase in the number of time steps required by standard explicit integrators. Stabilized RK methods are thus particularly efficient for the
time integration of large-scale, possibly nonlinear, parabolic PDEs \cite{DDD13}.
Several stabilized RK methods have been proposed in the literature, such as DUMKA methods, based on the composition of Euler steps \cite{Leb94,LeM94,Med98}, Runge--Kutta--Chebyshev (RKC) methods, based on the linear combination of Chebyshev polynomials \cite{SSV98,HoS80,VHS90} and orthogonal Runge--Kutta--Chebyshev methods (ROCK), based on optimal orthogonal stabilized functions \cite{Abd02,AbM01}; note that ROCK and RKC methods differ only beyond order one.
Still, when applied to \cref{eq:ode}, the number of stages $s$ of any standard stabilized RK method will, yet again, be determined by the stiffest part $\ff$ and the method eventually become inefficient.

To overcome the stringent step size restriction due to the cheap but stiffer part, $\ff$, while retaining the efficiency of explicit time integration for $\fs$, \emph{multirate methods}
use a smaller step size, or even an entirely different scheme, for integrating $\ff$.
Since the early work of Rice \cite{Ric60} and the work of Gear and Wells \cite{GeW84} who proposed a number of multirate strategies for the interlaced time integration of the ``fast'' and ``slow'' components using classical multistep schemes,
various explicit, implicit or hybrid multirate schemes have been developed based on Runge--Kutta methods using splitting techniques for ``fast'' and ``slow'' components or extrapolation techniques \cite{And79,EnL97,GKR01,GuR93,Hof76,Kva99,SHV07,SaM10,SkA89}.
All these methods require a predictor step and either interpolate or extrapolate between ``fast'' and ``slow'' state variables, which is prone to instability.
Although some of the implicit-explicit (IMEX) methods are provably stable, they are more cumbersome to implement and rapidly become too expensive as the number of ``fast'' unknowns increases.

More recently, Günther and Sandu exploited the generalized additive Runge--Kutta (GARK) framework \cite{SaG15} to devise multirate GARK (MrGARK) methods \cite{GuS16}. Many explicit, implicit and hybrid schemes are developed (up to fourth order accuracy) in \cite{RLS21,SRS19}, although some degree of implicitness is typically required to achieve a larger stability domain. In \cite{RSS20,San19,SeR19},
multirate infinitesimal step (MIS) methods \cite{KnW98,WKG09}, which assume the fast variables integrated exactly, are recast into the GARK framework and further generalized. The resulting multirate infinitesimal (MRI)-GARK methods extend exponential integrators \cite{HoO10} to the nonlinear case, where again the fast (nonlinear) dynamics are integrated exactly. When these semi-discrete schemes are used in practice, that is, in a fully discrete setting, some implicitness or very small steps sizes are typically required for stability.

Local adaptivity and mesh refinement are certainly key to the efficient numerical solution of PDEs with
heterogeneous media or complex geometry. Locally refined meshes, however, also cause a severe
bottleneck for any standard explicit time integration, as the maximal time-step is
dictated by maybe a few small elements in the mesh. To overcome the crippling effect
of local mesh refinement, various multirate (or local time-stepping) methods \cite{GMM15} were
proposed following the original local adaptive mesh refinement (AMR) strategy for first-order hyperbolic conservation laws by Berger and Oliger \cite{BeO84} --- see \cite{GaH13} for a review.
For parabolic problems, Ewing et al. \cite{ELV94,ELV90} derived and analyzed implicit
finite difference schemes when local refinement is utilized in space and time.
Dawson, Du and Dupont \cite{DQD91} combined implicit time integration in subdomains with
an explicit treatment of the interfaces, which leads to a decoupled but conditionally stable system.
In \cite{MAM06,ShV00}, various predictor-corrector and domain decomposition methods were combined to iteratively correct the solution or its boundary values at artificial interfaces.
By using static-regridding, Trompert and Verwer \cite{TrV91,TrV93a,TrV93b,TrV93c} developed a number of multirate time-stepping strategies for local uniform grid refinement (LUGR), where a first integration is performed on a global coarse grid and the accuracy is iteratively improved locally on
nested and increasingly finer subgrids.

In contrast to the above implicit, or locally implicit, multirate strategies,
fully explicit RKC time integration was recently combined with the AMR approach \cite{AbR20a,Mir17}
to tackle diffusion dominated problems. Again, the mesh is divided into two distinct regions,
the “coarse region,” which contains the larger elements
and corresponds to the mildly stiff part $\fs$, and the “fine region,” which contains the smallest elements
and thus corresponds to the severly stiff component $\ff$. In either subregion, the
number of stages $s$ in chosen according to the local mesh size, while
``ghost cell'' values at the coarse-to-fine interface are obtained by interpolating in time between stage values.
For certain problems, however, time interpolation of missing stage values from the other RKC method
can cause numerical instabilities \cite{AbR20a}.

To overcome the stringent stability condition due to a few severely stiff degrees of freedom, we first introduce in
\cref{sec:afequ} a {\it modified equation}, where the
spectral radius of its right-hand side, or {\it averaged force}, is bounded by that of the slower term $\fs$ yet still remains a good approximation of \cref{eq:ode}.
Evaluation of this averaged force requires the solution to a stiff, but cheap, auxiliary problem over short time and forms the basis of our multirate strategy.
In fact, the numerical integration of the modified equation by any explicit method via such a multirate approach
will be more efficient than integrating \cref{eq:ode} directly with the same explicit method.
In \cref{sec:mrkca}, we devise a multirate Runge--Kutta--Chebyshev (RKC) method explicit in both the fast and slow scales by utilizing
two different RKC methods to integrate the modified equation and evaluate the averaged force.
The resulting multirate RKC (mRKC) method assumes no scale separation,
requires no interpolation between stages, and remains accurate even if the roles of $\ff$ and $\fs$
change in time.
The fully discrete stability and accuracy analysis of the  mRKC scheme is given in
\Cref{sec:stab_conv}.
Finally, in \cref{sec:numexp}, we apply our mRKC method to a
series of test problems from both stiff ordinary and partial differential equations to demonstrate its usefulness and efficiency.


TITLE: References dealing with function spaces in pde like $C^k(\mathbb R^n)$, $C^\infty_c(\mathbb R^n), \ldots$?
QUESTION [2 upvotes]: What would be nice references for function spaces like $C^k(\mathbb R^n$), $C_0(\mathbb R^n)$, $C^\infty_c(\mathbb R^n), \ldots$ and most common function spaces which are offen employed in partial differential equations?
I'm looking references that deal with completeness, density, etc., more generally, the functional analytic aspects of those kind of spaces.
References for function spaces on manifolds would be also wellcome.
Thanks!

REPLY [1 votes]: I believe Lars Hormander's first volume of The Analysis of Partial Differential Operators covers relevant material to your first reference request. You might try Hans Triebel's Theory of Function Spaces also.