\section{Introduction}\label{intro}

In a recent letter,~\cite{KE99} Ku and Eguiluz have shown that
band structure effects are responsible for the positive
dispersion of the plasmon line width in K. Such a behavior had been
previously found by vom Felde, Spr\"osser-Prou, and
Fink~\cite{FSF89} in electron energy loss experiments. They
interpreted the result as due to short range electron-electron
correlations beyond the random phase approximation (RPA). However,
Ku and Eguiluz explain the result within the RPA in terms of the
extra plasmon decay channels provided by a manifold of unoccupied
d-symmetry bands.  The effect of these localized bands is
overlooked when treating K as an homogeneous electron system in
the {\it jellium} model. Therefore, band structure effects turn
out to be decisive in determining, not only the
existence interband transitions, but also  the
properties of plasmonic excitations. Indeed, {\it ab initio}
calculations of the dynamical response  have
properly described the experimental plasmon dispersion in a
variety of simple metals.~\cite{QE93,AK94,FSE97,MG94}
On the other hand, the {\it jellium} model
still provides a qualitative, though
not accurate, description of the  energy loss spectrum of
these systems in the long wavelength limit.

 The energy loss spectrum can be obtained either by doing electron
energy loss spectroscopy~\cite{R80} (EELS), inelastic scattering of X-ray
photons~\cite{SSSK95}(IXSS), or indirectly by means of optical
measurements.~\cite{W72} In noble metals this spectrum has a
complicate structure, which bears little resemblance
to that of simple metals. It is known~\cite{EP62}
that such a structure has to do with the existence of a manifold of
$d$-like bands a few eV below the Fermi level. These lead to important
deviations from even the crudest predictions for an
homogeneous electron liquid. Thus, whereas for Ag a narrow
resonance is observed at an energy $\approx 3.8$ eV, no
well-defined plasmonic excitations appear in the spectra of Cu and Au.

 The purpose of this paper is to investigate
the structure of the energy loss spectrum in the long wavelength
limit, where plasmonic excitations are expected to be important.
In particular, we have focused on Ag for our study, but
some of our conclusions do also apply to Cu and Au,
and could be extended to similar metallic compounds with
fully occupied $d$-like bands close to the Fermi level. The spectrum
of Cu was studied from first principles in a previous
work.~\cite{CRP99} In the present work,
we show that the structure of the energy loss spectrum
of Ag up to $\sim 10$ eV can be understood in terms of a
simple model. In particular, we are concerned with the role played by
the relative values of the Drude plasma frequency and the
threshold energy for the excitation of the $d$-like bands.
After clarifying the relevant elements for the existence
of a plasmonic excitation below the interband excitation threshold,
numerical experiments in an {\it ab initio} framework
for calculating the dielectric matrix of Ag
are performed to confirm this picture. Furthermore,
comparison  of the {\it ab initio} results with experimental data
gives good agreement for excitations well above threshold,
being able to reproduce all the structure in the loss
spectrum up to $\approx 30$ eV. The underestimate
of the threshold energy in local density
calculations also shows up in our calculations,
giving a lower plasma frequency than the experiments
(see below for a more detailed discussion).

 The question of how band structure alters the properties
of plasmons in metals has been previously addressed by a number of
authors. Wilson~\cite{W60} has studied the shift in the Drude plasma
frequency by the presence of an optical absorption band that lies
above (below) the free electron plasma frequency. He has found
that the plasma frequency is shifted
downwards (upwards) by the effect of the absorption band. More recently,
Sturm~\cite{S82} has studied within the RPA the corrections to the dielectric
function of an homogeneous electron gas due to a weak lattice
potential. He has shown that the lattice potential can provide the
necessary momentum  for a plasmon to decay into an electron-hole
pair below the Landau cutoff.~\cite{landau_cutoff} In other words,
when taking into account the lattice structure by folding the free
electron bands the first Brillouin zone (BZ), a
plasmonic excitation can decay by producing an interband
transition. As a consequence, even in th RPA, plasmons
acquire a finite lifetime. Core polarization effects on the plasma frequency
have been investigated by Zaremba and Sturm.~\cite{ZS85} They treated the
core polarizability as the superposition of the polarizability of the
isolated lattice ions. The polarization of core electrons lowers the plasma
frequency from its Drude value,  bringing it closer to the experimental value.

 Recently, the dielectric matrix as well as the energy loss
spectrum have been evaluated numerically within a first
principles framework. In these calculations, one usually starts from a set
of  LDA~\cite{KS65} orbitals and energies,
using time-dependent density functional theory (TDDFT),~\cite{tddft}
an approximation to the density correlation function
can be calculated. Hence, The energy loss spectrum is obtained
from the imaginary  part of the density correlation function, as dictated by
the fluctuation-dissipation theorem.~\cite{PN66}
To mention a few examples of theses studies, besides
the above mentioned work by Ku and Eguiluz~\cite{KE99},
Quong and Eguiluz~\cite{QE93} have thus investigated the
anisotropy in the plasmon dispersion of
aluminum~\cite{UR76}. Aryasetiawan and Karlsson~\cite{AK94} have
studied the excitations in the energy loss spectra of Li. The
negative dispersion of bulk plasmons in Cs has been explained by
Fleszar and co-workers.~\cite{FSE97} In general, specific details
of the band structure of metals as the presence of band gaps in
the conduction band region are needed to describe the experimental
energy loss spectra in some simple~\cite{MG94,F95} and noble
metals.~\cite{CRP99}

  The behavior of noble metals has been considered by
Ehrenreich and Phillipp.~\cite{EP62} These authors used the data
obtained  in experimental measurements of the reflectivity from Cu
and Ag~\cite{TP61} to obtain the optical and energy loss spectra.
They separated the Drude and
interband contributions to the optical response, thus shedding
light on the reason why the energy of the plasma resonance
in Ag is shifted down in energy. In a later work, Cooper,
Ehrenreich and Phillipp~\cite{CEP65} extended this analysis to Au.
They touched upon the question of assigning interband
transitions to the most relevant features observed in the optical
spectrum. To further clarify this point, Mueller and
Phillips~\cite{MP67} performed,  within the RPA and based on
band structure calculations  by  Burdick,~\cite{B63} a numerical
calculation of the imaginary part of interband contribution to the
optical response of Cu. Thus, they could give a more
correct identification of the interband transitions
that contribute to the optical spectrum at a given energy. Their
interpretation is consistent with the optical data obtained in a
number of reflectivity measurements from Ag and Ag alloys.~\cite{ML68,IHW70}

 As we have mentioned above, the energy loss function of Cu has been
recently  evaluated from first principles by Campillo, Rubio, and
Pitarke.~\cite{CRP99} These authors solved the
Kohn-Sham~\cite{KS65} equations on a plane wave basis using a norm
conserving pseudopotential.~\cite{TM91} Subsequently, they evaluated the
dielectric and energy loss functions. The method is shown to
give an accurate description of the energy loss spectrum of Cu once
the fully occupied $d$-like bands are taken into
account as part of the valence electron complex.
The remaining core electrons
do not take an active part in the excitations and  would only
act as a global dielectric background as long as the excitation
energy remains small compared to their binding energy.
In the present work, we obtain the energy loss
spectrum for Ag using the same {\it ab initio}
techniques but, this time, we focus on the plasmonic
excitations of noble metals. In the case of Ag,
the appearance of a plasma resonance at $\approx 3.8$
eV can be also qualitatively understood with help of a simple model,
which provides a new perspective on how a subthreshold plasmonic mode
can appear, and  what are the parameters controlling
the existence of this type of excitations.

  The outline of the paper is as follows. In the following
section,  we discuss a model which provides a qualitative
explanation for some features in the energy loss spectrum up to
$\sim 10$ eV. The details of the {\it ab initio} computation of
the spectrum can be found in Sect.~\ref{abinitio}. A
description and discussion of the results is given in Sect.\ref{results}.
Finally, the main conclusions of the present paper are summarized in
Sect.~\ref{conclusions}. Atomic units ($\hbar = e^2 = m_e = 1$) are
used in all mathematical expressions.

\section{A simple model}\label{model}

 In order to estimate the plasma frequency in Ag  the classical
expression $\omega_p = \sqrt{4 \pi n}$ can be used. If we set
$n = n_c$, the density of conduction electrons, i.e., those
coming from the atomic $5s$ orbital, we find that $\omega^c_p = 8.98$ eV.
However, if we furthermore include those coming from the atomic
$4d$ orbitals it turns out that $\omega_p = 32.6$ eV.
In both cases, we obtain values clear disagreement with the
experiment, which gives  a value of $3.8$ eV. Thus, as
we have remarked in the previous section,  Ag does not seem to
behave as a ``simple'' metal whose plasmonic excitations can
approximately understood within the {\it jellium} model.
Ehrenreich and Phillipp~\cite{EP62} considered this issue,
concluding that it is also necessary to take
into account the  effect of electrons coming from the atomic $4d$ orbitals in a more
subtle way. In fact, it is necessary to account for the full band
structure of Ag in order to explain what is observed in the
experiments. In what follows, we  shall study how the interband
transitions affect the collective motion of conduction electrons,
leading to the existence of a subthreshold plasmonic excitation.

 The usual way to find out the energy plasmonic excitations goes as follows.
Disregarding crystal local field effects,
plasmonic excitations correspond to zeroes in the
macroscopic dielectric function. In the $|{\bf q}| = 0$ limit
we have that:
\begin{equation}\label{meq1}
\epsilon(\Omega) = 0,
\end{equation}
for some complex $\Omega$. Using the Sellmeyer-Drude
expression~\cite{P65} for $\epsilon(\omega)$, and
neglecting damping effects,
\begin{equation}\label{meq2}
\epsilon(\omega) = 1 - \frac{(\omega^{D}_p)^2}{\omega^2} +
 \delta\epsilon^{(ib)}(\omega),
\end{equation}
where last term corresponds to the interband contribution:
\begin{equation}\label{meq3}
\delta\epsilon^{(ib)}(\omega) = \sum_{\omega_n > 0}
    \frac{f_n}{\omega_n^2 - \omega^2}.
\end{equation}
The condition in Eq.~(\ref{meq1}) can be written as
\begin{equation}\label{meq4}
\Omega^2 - (\omega^{D}_p)^2 + \Omega^2 \delta\epsilon^{(ib)}(\Omega) = 0.
\end{equation}
The oscillator strength, $f_n$ in Eq.~(\ref{meq3}), is a measure of the effective
number of electrons that participate in the excitation of a given
mode `n', with energy $\omega_n$. The dominant contribution at
frequencies below $\omega^c_p$ comes from the excitation of conduction
electrons. The Drude plasma frequency~\cite{EP62} $\omega^{D}_p =
9.2$ eV, is slightly higher than $\omega^{c}_p$.
The ratio of the two plasma frequencies is usually
expressed in terms of the ``optical'' mass, $m_{opt}$ (recall that
 $m_e=1$ in atomic units):
\begin{equation}\label{meq5}
\left( \frac{\omega^{c}_p}{\omega^{D}_p} \right)^2  = m_{opt} = 0.95
\end{equation}
for Ag and $m_{opt} = 1.35$  for Cu.~\cite{EP62} In the RPA, $m_{opt}$
corresponds to the effective mass of electrons averaged
over the occupied part of the conduction band.

   Eq.~(\ref{meq4}) can be read in the following way. We assume that
the term $\Omega^2 \delta\epsilon^{(ib)}(\Omega)$ plays the role of  a
``self-energy'' for the Drude plasmon. In simple metals, this
term is usually small and varies slowly with $\omega$.
This leads to a small shift of the Drude frequency. The plasma
frequency also acquires a small imaginary part, which accounts
for the fact that the plasmon state is degenerate in energy
with other configurations of the system
involving one or several electrons excited from one band to
another. This idea can be easily realized in a model in which the
plasmon is viewed as a boson coupled to an energy continuum. In an
independent electron picture, the continuum corresponds to the
electron-hole pair excitations.  A simple model
Hamiltonian that can be set up
to describe this situation  is the following :
\begin{eqnarray}\label{meq6}
H  &=& \sum_{\bf q} \; \omega^{D}_p \:
       b^{\dagger}_{\bf q} b_{\bf q}   \nonumber\\
    &+& \sum_{k,{\bf q}} \; \omega_{k,{\bf q}}\:
        m^{\dagger}_{k,{\bf q}} m_{k,{\bf q}}\nonumber \\
    &+& \sum_{{\bf q},k} \;
               g_{k,{\bf q}} \, (b_{\bf q} + b^{\dagger}_{\bf -q})
     (m^{\dagger}_{k,{\bf q}} + m_{-k,{\bf -q}} ).
\end{eqnarray}
In this model, plasmons are presented as a field of oscillators carrying
momentum $\bf q$. Thus, $b^{\dagger}_{\bf q}$ and $b_{\bf q}$ are
plasmon creation and annihilation operators, respectively. To keep
the model as simple as possible, but retaining all the important
physical effects, the energy continuum is also represented
by a set of harmonic oscillators whose quanta are created and
annihilated by the operators $m^{\dagger}_{k,{\bf q}},\, m_{k,{\bf
q}}$. The label $k$ stands for the additional quantum numbers
(spin, relative momentum, band indices...) carried by a given
continuum mode.  The operator $m_{-k,{\bf -q}}$
($m^{\dagger}_{-k,{\bf -q}}$) anhilates (creates) a continuum mode with all
the momentum quantum numbers reversed. The assumed interaction
between plasmons and the energy continuum conserves the momentum
$\bf q$. Terms connecting plasmons and modes with momentum
differing by a finite reciprocal lattice vector are therefore
neglected. This approximation amounts to disregarding crystal
local field corrections, which in general is a good approximation
for metals. Indeed, in the frequency range
in which we are interested, namely $\omega \sim \omega^{D}_p$ , as
it can be seen from the {\it ab initio} calculations in Cu\cite{CRP99} and Ag (see
below), the local field effects at small ${\bf q}$ amount
to a few percent correction.

 The above model, Eq.~(\ref{meq6}), can be solved exactly.~\cite{G81,M90}
Indeed, it is related to the models introduced by Fano~\cite{F61} and
Anderson~\cite{A61} to study the coupling of a discrete state to
a continuum. More precisely, the present model can be thought of as an extension
of Gadzuk's work~\cite{G81} for a  localized vibrational mode
in an electron liquid. Thus, the energies of the plasmonic
modes are given by the solutions of the following equation:
\begin{equation}\label{meq7}
\Omega^2 - (\omega^{D}_p)^2 + 2 \omega^{D}_p M({\bf q},\Omega) = 0,
\end{equation}
where $M({\bf q}, \omega)$ is the plasmon self-energy. It can be
expressed in terms of its imaginary part:
\begin{eqnarray}\label{meq8}
{\rm Im}\; M({\bf q},\omega) = \pi
   \sum_{k} | g_{k,{\bf q}}|^2 \: [&&
      \delta (\omega - \omega_{k,{\bf q}}) \nonumber \\
    &&  - \delta(\omega + \omega_{k,{\bf q}}) ],
\end{eqnarray}
by means of the expression:
\begin{equation}\label{meq9}
M({\bf q},\omega) = \int \; \frac{{\rm d}\omega'}{\pi}\,
  \frac{{\rm Im}\; M({\bf q},\omega')}{\omega' - \omega - i\eta},
\end{equation}
with $\eta \to 0^+$.

Working by analogy with Eq.~(\ref{meq4}), leads us to tentatively
identify
\begin{equation}\label{meq10}
M(|{\bf q}| \to 0,\omega ) =
  \frac{\omega^2}{2 \omega^{D}_p} \delta \epsilon^{(ib)}(\omega).
\end{equation}
If the RPA expression of $\delta\epsilon^{(ib)}(\omega)$ is used,
we can express it in terms of single-electron
orbitals, $|n\: {\bf k}  \rangle$, and energies,
$\varepsilon_{n {\bf k}}$:
\begin{eqnarray}\label{meq11}
 {\rm Im}&& M(|{\bf q}|\to 0, \omega) =
\frac{2 \pi^2}{3 V \omega^D_p}  \sum_{{\bf k} \in BZ} \sum_{m \neq n}
\, f(\varepsilon_{n {\bf k }}) \nonumber \\
&& \times  \left[ 1 - f(\varepsilon_{m {\bf k }} ) \right]
 | \langle m \: {\bf k} | {\bf p} | n\: {\bf k}  \rangle |^2
 \delta(\omega - \varepsilon_{m {\bf k }} + \varepsilon_{n {\bf k }} ),
\end{eqnarray}
where $V$ is the volume of the system and ${\bf p}$ the momentum
operator. We have assumed cubic symmetry and a local pseudopotential;
if the pseudopotential is not local, the
matrix element of $\bf p$ must be replaced by a more complicated
expression. For $\omega = \omega^{D}_p$ this equation matches the
Golden Rule result for the plasmon line width in the $|{\bf q}|=
0$ limit,~\cite{BB94} which confirms the identification made
above in Eq.~(\ref{meq10}). For arbitrary $\omega$, however, it is
better to look at ${\rm Im}\; M({\bf q},\omega)$
as a weighted density of the continuum modes, cf. Eq.~(\ref{meq8}).
Thus, it follows from Eqs.~(\ref{meq7}) and (\ref{meq9})
that if $\omega^{D}_p$ lies above (below) the region where
${\rm Im}\; M({\bf q},\omega)/\omega^D_p$ is large
(provided  ${\rm Im}\, M({\bf q},\omega^D_p)/\omega^D_p$ be small), the
plasma frequency will be blue (red) shifted. This
result is in agreement with the above mentioned work of Wilson.~\cite{W60}
When ${\rm Im}\, M({\bf q},\omega^D_p)/\omega^D_p$  is
large, we will not observe a well-defined plasma resonance
in the spectrum.  This seems to be the case of noble metals due to
the presence of a manifold of occupied $d$-like
bands, which gives rise to a broad band of modes extending
from the interband threshold energy  $\omega_T \approx 3.9$  eV
(for Ag and $2.1$ eV for Cu), to well above the Drude
plasma frequency ($9.2$ eV for Ag and $9.3$ eV for Cu).
If we assumed that the onset of interband
transitions in Ag took place so sharply as to produce a
discontinuity in the weighted density of continuum modes:
\begin{equation}\label{meq12}
{\rm Im}\, M(|{\bf q}|\to 0, \omega) \sim \theta(\omega_T - \omega)
\end{equation}
for $\omega \sim \omega_T$, then ${\rm Re}\, M({\bf q}|\to 0,
\omega)$ would develop a logarithmic behavior $\sim -
\log(\omega - \omega_T)$ for $\omega < \omega_T$. This means that
Eq.~(\ref{meq7}) would have a dampingless solution just below the
threshold, $\omega_T$. This polaronic-like solution~\cite{ES63}
corresponds to an hybrid plasmonic mode, in which the conduction
electrons oscillate coherently without exciting the
electrons in the $d$-like bands. These are just polarized by the
electric field set up by plasmonic mode,
so that the plasma frequency is lowered.
This point of view agrees well with more phenomenological
approaches,~\cite{L93,SR97} which simply assume that the red shift in the plasma frequency
can be accounted for by an effective dielectric function $\epsilon_d$, for
the electrons in the $d$-like bands so that $\omega^{*}_p =
\omega^D_p/\sqrt{\epsilon_d} \approx 3.8$ eV.

 Additionally, there is another complex solution of Eq.~(\ref{meq7}),
for which the conduction electrons  move as an overdamped oscillator,
rapidly decaying into an electron hole-pair. For $|{\bf q}| \to 0$, this
excitation corresponds to an optical interband transition.
Two types of optical  transitions  may occur above the
threshold.~\cite{MP67} Two examples have been indicated by arrows in Fig.~\ref{fig1}, which
represents the calculated first-principles LDA band structure of Ag.
Type A corresponds to the excitation of an electron from one
flat $d$-like band to an unoccupied state in the conduction band or above. The other
type (B), whose threshold occurs at approximately the same energy in the three
noble metals, corresponds to the excitation from an occupied
state in the conduction band (near $L$ point in the BZ) to the flat
part of an unoccupied band or above. This type of transitions
seems to be responsible for the
damping of the plasmon when the $4d$ electrons are
included in the pseudocore, i.e., they are considered to be frozen.
In Sect.\ref{results}, we shall explain this issue more in detail.
These two types of transitions are important because
they often involve flat bands, which may lead to van Hove singularities in the
joint density of states.

 In the case of Ag, the onsets for the transitions of type A and B
nearly overlap in energy. This produces a sharp interband onset and
gives rise to the pronounced peak  that ${\rm Im}\,
\epsilon(\omega)$ exhibits around $4$ eV.
From Eq.~(\ref{meq10}), it follows then that ${\rm Im}\, M(|{\bf q}|\to 0, \omega)$ will also
have a maximum at the same energy. The onset of the interband transitions at
$\omega_T$ is not as sharp as suggested by Eq.~(\ref{meq12}).
However, it still leads to a zero of $\Omega^2 - (\omega^D_p)^2  + 2 \omega^D_p {\rm
Re}\,M(|{\bf q}|\to 0, \Omega)$ for a real $\Omega$, which occurs
below the threshold, where the density of continuum modes is very small.
This point is illustrated by Fig.~\ref{fig2}, which displays
the graphical  solution of Eq.~(\ref{meq7}).
This provides an explanation for the
existence of a subthreshold plasma resonance in Ag. However,  it remains to
be explained why a similar phenomenon is not observed in the
spectra of Cu or Au. The difference stems from the specific details of the
band structure and, in particular, from the value of the threshold energy. In
Fig.~\ref{fig3}, we have plotted the graphical solution of Eq.~(\ref{meq7})in
the case of Cu. The threshold for the excitation of $d$-like bands
occurs this time at a lower energy ($2.1$ eV) because the onset for the
interband transitions of type A and B does not overlap.
Comparison of Figs.~\ref{fig2} and \ref{fig3} shows the
important role played by the value of the threshold energy,
cf. Eq.~(\ref{meq10}), in enhancing the maximum of ${\rm Im}\, M(|{\bf q}|\to 0, \omega)$,
and thus in ${\rm Re}\, M(|{\bf q}|\to 0, \omega)$.
Therefore, for Cu no  zero of Eq.~(\ref{meq7}) close to the real axis exists
below the threshold, where this plasma mode would have a small
density of continuum modes to decay. A similar conclusion can be drawn
for the case of Au.

  To sum up, band structure effects
modify the original Drude plasmon, which
turns out to be no longer a well defined plasmonic excitation. Thus,
it shows up in the Ag spectrum as broad peak with a maximum around $8$ eV.
Moreover, a combination of a high value for $\omega_T$
along with the sudden onset of interband transitions
produces a narrow plasma resonance in Ag. In this respect, our conclusions agree with those of Ehrenreich
and Phillipp.~\cite{EP62} Furthermore, the line width of this ``delicate'' feature in
the energy loss spectrum also depends on the sharpness of the interband onset.
The sharper onset the narrower the resonance, which in any case has very little
oscillator strength.  This makes this plasmonic excitation very
sensitive to changes in the band structure,
impurities or defects. This is consistent
with the experimental data from alloying experiments,~\cite{W72,ML68,IHW70}
which indicate that the plasmon is strongly damped by increasing the concentration
of the other component in the alloy. Finally, it is worth pointing out
that the existence of this resonance also depends on the value of the Drude
plasma frequency. If the value of $\omega^D_p$ is increased,
for example by increasing the electron density in the system, the
resonance could disappear from the spectrum.
This prediction is confirmed by our {\it ab initio} results (shown below)
for the energy loss spectrum of Ag under externally applied
hydrostatic pressure.

\section{Details of the  {\it ab initio} calculations}\label{abinitio}

 By measuring the energy lost by electrons or
X-ray photons in their interaction with matter, one can access the
dynamical structure factor, $S({\bf q},\omega)$ of a system. The
fluctuation-dissipation theorem~\cite{PN66} allows us to relate
$S({\bf q},\omega)$ to the density correlation function:
\begin{equation}\label{aeq1}
\chi({\bf r}, {\bf r'}, t) = -i \theta(t) \;
  \langle \left[ \delta n_H({\bf r}, t), \delta n_H({\bf r}',0)
        \right] \rangle,
\end{equation}
where  $\delta n_H({\bf r}, t) = n_H({\bf r}, t) - \langle n_H({\bf
r}, t) \rangle$, and suffix $H$  means that  operators are in
the Heisenberg picture~\cite{M90}; the brackets mean that the average
is taken over the ground state of the system. Now, time-dependent
density functional theory~\cite{tddft} (TDDFT) provides
us with a method to compute this correlation
function by solving the following integral equation:
\begin{eqnarray}\label{aeq2}
&& \chi_{{\bf G}{\bf G'}}({\bf q},\omega)
   = \chi^0_{{\bf G}{\bf G'}}({\bf q},\omega)
    + \sum_{{\bf G}_1,{\bf G}_2} \,
       \chi^0_{{\bf G}{\bf G}_1}({\bf q},\omega) \nonumber \\
&& \times \left[ v_{{\bf G}_1{\bf G}_2}({\bf q})
 +  K^{xc}_{{\bf G}_1{\bf G}_2}({\bf q},\omega) \right]
     \chi_{{\bf G}_2{\bf G'}}({\bf q},\omega).
\end{eqnarray}
Here, we have exploited crystal symmetry and time translation invariance
by introducing:
\begin{eqnarray}\label{aeq3}
\chi_{{\bf G}{\bf G'}}({\bf q},\omega) = \frac{1}{V}
 \int {\rm d}t \; &&
  {\rm d}^3{\bf r} \; {\rm d}^3{\bf r'} \,
  e^{i \omega t}
  e^{-i({\bf q} + {\bf G})\cdot {\bf r}} \nonumber \\
 && \times e^{i({\bf q} + {\bf G'})\cdot {\bf r'}}\;
  \chi({\bf r}, {\bf r'},t),
\end{eqnarray}
where $\bf q$ belongs to the BZ, and
${\bf G},{\bf G'}$ to the reciprocal lattice.
Similar definitions hold for the other functions that appear in
Eq.~(\ref{aeq2}). For example,
\begin{equation}\label{aeq4}
v_{{\bf G}{\bf G'}}({\bf q}) =
  \frac{4\pi}{|{\bf q} + {\bf G}|^2} \delta_{{\bf G},{\bf G'}}
\end{equation}
corresponds to the Fourier components of the Coulomb interaction between metal
electrons. The function
\begin{eqnarray}\label{aeq5}
\chi^0_{{\bf G}{\bf G'}}({\bf q},\omega) = \frac{1}{V}\:
  \sum_{{\bf k} \in {\rm BZ}} && \sum_{n,m}\,
  \frac{ f(\varepsilon_{n{\bf k}}) - f(\varepsilon_{m{\bf k}+{\bf q}}) }
    {\omega + \varepsilon_{n{\bf k}} - \varepsilon_{m{\bf k}+{\bf q}}
 + i \eta} \nonumber  \\
&& \times \langle n\: {\bf k} | e^{-i({\bf q} + {\bf G})\cdot {\bf r}}
 | m\: {\bf k} + {\bf q} \rangle \nonumber  \\
&& \times
\langle m\: {\bf k} + {\bf q} | e^{i({\bf q} + {\bf G'})\cdot {\bf r}}
  | n\: {\bf k}  \rangle
\end{eqnarray}
with $\eta \to 0^+$, is the density correlation function
for a system of independent particles, with single-particle energies
$\varepsilon_{n{\bf k}}$ and orbitals  $|n\: {\bf k} \rangle$.
These are eigenvalues and eigenstates of a single-particle
Hamiltonian. The potential energy in this Hamiltonian is given by
the Kohn-Sham~\cite{KS65} potential. In our case, the Kohn-Sham potential has been
calculated within LDA~\cite{LDA,KS65}.

  We have used a plane wave basis to expand the Bloch eigenstates.
For this procedure to be computationally efficient, it is
necessary to replace the lattice potential by a superposition of
ionic pseudopotentials having the same scattering properties.
However, as electrons in the outer atomic $4d$ shell play
an important role in the response of
noble metals, they are treated as a part of the valence electrons. The
remaining core electrons are replaced by a frozen pseudocore
using a norm-conserving pseudopotential scheme.~\cite{TM91}
Therefore, electrons included in the pseudocore
will not respond to an external perturbation.
This is not an important drawback as long as the
excitation energy remains below the typical excitation
energies of core electrons $\sim 100$ eV.

 Including $4d$ electrons in the valence produces a pseudopotential
that requires a relatively high number of plane waves to get good
convergence. Here, we  have used a scalar relativistic
pseudopotential generated according to the scheme
of Troullier and Martins.~\cite{TM91} It requires $\approx 1,400$
plane waves per Bloch state up to a cutoff energy of $\approx 1.1$ keV
in order to obtain a well converged ground state energy (less than $0.01$ eV),
and single-electron energies and orbitals
up to $\approx 100$ eV. The lattice constant has been set to
the experimental value of the FCC crystal structure of Ag,
namely $a = 4.09\, {\rm \AA}$. We have also solved the Kohn-Sham equations
with $a = 3.70\, {\rm \AA}$, which corresponds to an
externally applied hydrostatic pressure of $63$ GPa.

 The RPA is equivalent to solving Eq.~(\ref{aeq5})
with  $K^{xc}=0$. Therefore,
all the exchange correlation effects beyond the RPA
are contained in the kernel $K^{xc}_{{\bf G}{\bf G'}}({\bf q},\omega)$.
Consistently with the use of LDA eigenstates and eigenvalues in
the calculation of $\chi^0$, this
function is approximated by the following expression:
\begin{equation}\label{aeq6}
K^{xc}_{{\bf G}{\bf G'}} =  \frac{1}{V}
 \int {\rm d}^3{\bf r} \; e^{- i ({\bf G} - {\bf G'}) \cdot {\bf r} } \;
f_{xc}({\bf r}),
\end{equation}
where
\begin{equation}\label{aeq6b}
f_{xc}({\bf r}) = \left[ \frac{dV_{xc}(n)}{dn}\right]_{n_{gs}({\bf r})}.
\end{equation}
Here $V_{xc}$ is the exchange correlation part of the
Kohn-Sham  LDA potential, and $n_{gs}({\bf r})$ is the ground state
density. This approximation  is usually called {\it
adiabatic} local density approximation (ALDA). However, it is
observed that for metals in the long wavelength limit, ALDA and
RPA yield very similar results. This has been found in simple
metals,~\cite{KE99,QE93} and  also in Cu.~\cite{CRP99}
Thus, when computing the
spectrum, we have used the RPA.

 Instead of solving Eq.~(\ref{aeq5}), we have evaluated the
RPA dielectric matrix:
\begin{equation}\label{aeq7}
\epsilon_{{\bf G}{\bf G'}}({\bf q},\omega) =
\delta_{{\bf G'},{\bf G'}} - v_{{\bf G}{\bf G'}}({\bf q})
  \chi^0_{{\bf G}{\bf G'}}({\bf q},\omega).
\end{equation}
In the RPA, the dynamic structure factor $S({\bf q,\omega})$ is  proportional
to the loss function, ${\rm Im}\,\left[-\epsilon^{-1}_{\bf 0
0}({\bf q},\omega)\right]$, which is the ${\bf G} = {\bf G'} = {\bf 0}$
element of the inverse dielectric matrix.  To calculate the dielectric matrix
we have used a cutoff of $41$ eV in ${\bf G},{\bf G}'$, which amounts to
using $59\times59$ matrices.

 To  evaluate  $\chi^0$ using Eq.~(\ref{aeq5}),
we have cut off the sum over $n$ and $m$ using
bands up to an energy of $\approx 100$ eV. The integration over
the first BZ has been performed using a $20\times 20\times 20$ Monkhorst-Pack
mesh~\cite{MP76} in the irreducible BZ, which allows to deal
effectively with the topology of the Fermi surface of Ag. Indeed,
this is a very dense mesh and the
calculation of $\chi^0$ becomes certainly time and memory consuming. However, if
one is to resolve the very narrow plasma resonance that Ag
presents near the interband excitation threshold, it is then
necessary to work with such a fine mesh.

In order to further reduce the numerical damping
$\eta$ we have
also computed $\chi^0$ on the imaginary frequency axis by
making the replacement $\omega + i \eta \to i \nu$ in Eq.~(\ref{aeq5}).
Analytic continuation to the real axis has been carried out by using a Pad\'e
approximant.~\cite{LC96} This procedure allows us to use a
numerical damping as small as $0.001$ eV. This is to be compared
with the $0.2$ eV value used in the real frequency calculations
also presented in this work. However, we  shall show the
results of boht types of calculations together, to
indicate that one must be very careful with the analytic
continuation procedure using Pad\'e approximants. Thus,
the results obtained using this method tend to be
smoother than what is obtained with a calculation with real
frequencies and larger values of $\eta$. Thus, it usually happens
that some features of the spectrum, as obtained in a real frequency calculation,
are completely missing from the
analytically continued result (see following section). In fact,
although we have used Lentz algorithm~\cite{PTVF92}
to evaluate the corresponding continued fraction, we have been unable to ensure that the procedure
converges by increasing the number of points over the same interval of $\nu$.
When fitting the Pad\'e approximant~\cite{points} to the calculated values of
$\epsilon_{{\bf G},{\bf G}}({\bf q}, i\nu)$, we have
found that the results are unstable with respect to increasing
the number of points. Thus, we have obtained differences of
the order of $10$ \% when the number of points over the same
interval was simply doubled. Finally, we have
decided to use the modulus of the remainder
of the truncated continued fraction
to estimate  the error, and always compare the outcome with a
real frequency calculation prior to consider it as reliable.

\section{Results}\label{results}

 In Fig.~\ref{fig1} we have plotted the band structure of
Ag as calculated within the LDA. Two examples of the main types
of interband transitions (A and B) which couple
to  plasmons at energies $\omega \sim \omega^D_p$ are indicated by arrows.
In the LDA, we obtain a threshold
energy for  type A transitions $\approx 3$ eV, whereas experimentally
it is found $4$ eV.  Since  LDA eigenvalues enter the expression for $\chi^0$,
this turns out to be the reason for underestimating  the
threshold frequency, $\omega_T$, as shown in
Figs.~\ref{fig4} and \ref{fig5}. An underestimate
of this kind was already observed for Cu, where $\omega_T \approx
1.5$ eV in contrast with the experimental value of $2.1$ eV.
The source of such an underestimate may be
either a failure of the LDA or an indication that a frequency dependent and perhaps
nonlocal $K^{xc}$ may be required to fully account for the
experimental value of $\omega_T$. In our opinion, however, it
is related to a failure of the LDA, which may not
be taking into account the strong correlations occurring amongst
the electrons in the flat $d$-like bands of Ag. Therefore, a more detailed
description of correlation for the $d$-like bands seems to be
required.\cite{Jorge}

Fig.~\ref{fig4} presents our results for the energy loss function
for  ${\bf q} = (0.05\, 0\, 0) (2 \pi / a)$ as obtained from the
LDA band structure. Since our interest is focused on the plasmonic excitations,
which occur for small ${\bf q}$, we need not consider exactly the
$|{\bf q}| = 0$ limit. Thus, for the small $\bf q$ considered,
the theory developed in Sect.~\ref{model} still applies.
The continuous line  in this figure corresponds to
${\rm Im}\:\left[-\epsilon^{-1}_{\bf 0 0}({\bf q},\omega)\right]$, whereas the
dashed line corresponds to ${\rm Im}\, \left[-1/\epsilon_{\bf 0 0}({\bf
q},\omega)\right]$, i.e. with and without including crystal
local field effects respectively. As we have already said
in Sect.~\ref{model}, crystal local field effects are at most
a correction $\sim 10\,\%$ in the considered frequency range. Dots are the experimental
data~\cite{P85} obtained by means of optical measurements (i.e., for
$|{\bf q}|  \to 0$), and  should be compared with the dashed line.
As a consequence of  the LDA underestimate of $\omega_T$,
the plasmon peak appears at a lower energy
$\approx 3.4$ eV than observed experimentally~\cite{EP62}
$\approx 3.8$ eV. This peak also looks more damped either
because the high value of the numerical damping
employed in the calculation and  because
of the LDA underestimate of $\omega_T$.

 From Figs.~\ref{fig4} and \ref{fig5} we can conclude that:
(i) The Drude plasmon becomes so broadened  by
its coupling to interband transitions that we can no longer speak of
a well defined excitation. The high energy tail of this broad
feature around $8$ eV is more strongly affected by crystal local
field effects. This seems to  indicate that those decay processes
involve the excitation of  an electron in a more localized state,
probably lying deeper in the $d$-like band. The important role played
by the $d$-like bands in destroying the Drude plasmon can be seen in
Fig.~\ref{fig6}, where the dashed line corresponds to a calculation
of the energy loss function putting the ten $4d$ electrons of atomic
Ag in the frozen pseudocore. In this case, the Landau damping of the
Drude plasmon is only due to transitions of type B, which have a much
smaller spectral weight than those of type A. Furthermore, the
plasma frequency is pushed upwards by this continuum of modes, which
mainly lies below $\omega^D_p$. However, when electrons from the $d$-like
bands are considered as active, the effect goes in the opposite direction
thus shifting the maximum of the Drude peak  from $\omega^D_p$ eV to $\approx 7.5$ eV;
(ii) In agreement with the experiments,
a double hump structure in the region between 15 to 30 eV is found.
This is a band structure effect, which turns out to be a general feature
in the noble metal systems.~\cite{CRP99}
As it can be seen in Fig.~\ref{fig5}, although the
calculation using a Pad\'e
approximant (dot-dashed line) gives a reliable description  of
the low frequency part, it fails to reproduce the double hump structure.
This seems to be  a rather general feature of analytic continuation
by Pad\'e approximants.  Whenever the function to be continued
presents two poles close to each other and having a large imaginary part,
the values for the function on the real axis correspond  to an
average of the two poles, and a single broad peak shows up in
the outcome of the analytical continuation.  (iii) Altogether,
the {\it ab initio} LDA calculations presented here give a
good description of the high energy features (i.e. above $\omega^D_p$)
of the energy loss spectrum of Ag. These include the
Landau damping of the Drude plasmon and
the double hump structure mentioned above. This could be expected
on the basis that the spectral weight
in the high energy region of the spectrum comes,
to a large extent, from the excitation
of electrons to bands well above the Fermi level, which are
expected to be more accurately reproduced by the LDA,
as it is an approximation suitable for an homogeneous
electron system. However, the LDA underestimate
of the threshold energy is responsible for
failing to reproduce the low energy structure.
This issue is further investigated below.

 In Fig.~\ref{fig6}, we show the results for the real
and imaginary parts of $\epsilon_{\bf 0 0}({\bf q}, \omega)$
for ${\bf q} = (0.05 \, 0 \, 0) (2 \pi /a)$. The same remarks
made above apply in this case as well. The experimental
results, full circles for the imaginary part and  open squares for the real part
are from Ref.~\onlinecite{P85}. In this figure we can see more
easily that the onset for the excitation of interband transitions
occurs, not only at a lower energy, but that it is less sharp as well.
As it has been remarked in the previous section, this leads to a
more damped plasmonic mode below the threshold.

  In order to investigate the influence of the interband threshold
energy on the existence of a subthreshold plasma resonance,  we
have performed a series of numerical ``experiments''. The
motivation of these experiments is to understand why Ag presents a
narrow plasma resonance in its energy loss spectrum, while
Cu (and Au) does not, despite the similarity of their band structures.
To modify the threshold energy $\omega_T$ without opening unphysical
gaps in the band structure, we have scaled by a factor $\alpha$
the separation between the Fermi energy ($\varepsilon_F$)
and the LDA eigenvalues of the occupied states:
\begin{equation}\label{aeq8}
\varepsilon'_{n{\bf k}} = \varepsilon_F +
\alpha (\varepsilon_{n{\bf k}} - \varepsilon_F).
\end{equation}
Thus, using this new set of eigenvalues, we have calculated $\chi^0$
and hence obtained the energy loss function. This rather {\it ad hoc}
procedure does not modify the wave functions accordingly. Furthermore, it
introduces a discontinuity of $\nabla_{\bf k} \varepsilon'_{n{\bf k}}$ at the
Fermi surface, so that it might lead unphysical results. However,
as shown in Fig.~\ref{fig7},
this seems not to be the case, at least for $\omega \sim \omega_T$ and
$\bf q$ finite. As predicted in Sect.~\ref{model}, the subthreshold
plasmon peak is enhanced and shifts up in energy
as $\omega_T$ increases. In this figure, continuous lines
correspond to calculations for real frequencies with
$\eta = 0.2$ eV, while dashed lines
correspond to the analytic continuation using a Pad\'e approximant
from the results obtained on the imaginary frequency axis.
Both calculations show the same tendency. On the other hand, when
the threshold energy is further decreased, from its LDA value $\approx 3$ eV,
to $\omega_T \approx 1.5$ eV,  the plasma resonance disappears
and the spectrum becomes more ``Cu-like''.

  Finally, in Fig.~\ref{fig8} we show the  energy loss spectrum of Ag
under an external pressure of 63 GPa. As in previous cases, the
continuous line corresponds to a calculation using real frequencies
while the dashed line is the result of an analytical continuation using
a Pad\'e approximant. Notice that in this case the
lattice constant is $a_{press} = 3.7\, {\rm \AA}$ so that now $q = (0.05\,  0 \, 0)
(2 \pi / a_{press})$, corresponds to a slightly larger momentum
($|{\bf q}_{press}| = 8.5 \times 10^{-2}\, {\rm \AA}^{-1}$ vs
$|{\bf q}| = 7.7 \times 10^{-2}\,{\rm \AA}^{-1}$).

Decreasing the lattice constant changes the width of the bands so that they
disperse more rapidly with $\bf k$ (i.e., the band width becomes larger).
Furthermore, electronic density also increases. This shifts the
maximum previously at $\omega^{max} \approx 7.5$ eV to
$\omega^{max}_{press}\approx 8.8$ eV. Since this maximum
is related to the Drude plasmon, whose {\it bare} frequency scales with
the lattice parameter as $a^{-3/2}$, we could understand the shift
(neglecting the change in the optical mass) as due to the
increase in the electr\'onic density. Indeed, this argument
works reasonablely well as $(a_{press}/a)^{-3/2} = 1.16$, and
$\omega^{max}_{press}/\omega^{max} \approx 1.17$~. Moreover, as the
bands become wider, specially those corresponding to excited states,
more oscillator strength is transferred to higher energies.
This decreases the density of modes at the Drude frequency,
and thus the corresponding peak is now less broad. What is more, in this spectrum
the subthreshold plasmon is not present any more. Again,
this is due to the increased electronic density, which  yields
a higher value for the Drude plasma frequency so that Eq.~(\ref{meq7})
has no longer a solution close to the real axis. The behavior
under pressure predicted here remains to be confirmed by experiments.

\section{Conclusions}\label{conclusions}

  Here, we summarize the main conclusions of the present work.
We have shown that a first principles calculation of the energy loss
spectrum is able to reproduce, within the RPA, the high energy structure
in the spectrum of Ag. However, since the LDA underestimates the
interband energy threshold, the agreement is not as good as for the high
energy region around $4$ eV. For example, the plasma resonance observed
experimentally at $3.8$ eV appears as  a more damped peak at $3.4$ eV.
In our opinion, this underestimate may be related to a failure of the LDA to
account for  electron correlation effects that take place in the flat d-like
bands of Ag. In fact, an {\it ad hoc} scaling of the occupied bands
so that the threshold energy is brought closer to the
experimental value gives rise to a better defined
subthreshold plasmonic excitation. We  have also shown
the important role played by the $d$ symmetry
bands in determining the existence of the plasma resonance at $3.8$ eV, just
below the interband excitation threshold. By using a simple model, it is shown
that the appearance of this excitaci\'on has to do with the sharpness of
the interband onset and the relative values of the threshold energy and the Drude
plasma frequency. Since it contains very little oscillator strength, it is
very sensitive to impurities, defects, and changes in the band
structure. Thus, we have demonstrated by performing numerical experiments that its presence
depends on the value of the threshold energy.
Finally, we have also computed the energy loss spectrum
of Ag under pressure, and found that its features can
be aslo understood in terms of the simple model presented here. We remark that
the techniques used and analysis carried out in this work are not specific
of Ag, and could be readily extended to other metallic compounds with
fully occupied $d$ symmetry bands close to the Fermi level.

\acknowledgements

 We thank N. Lorente,  J.~F. Dobson, and L. Hedin  for fruitful
discussions, and I. Campillo for many useful suggestions.
This work was supported by the Basque {\it Unibertsitate
and Hezkuntza Saila}, JCyL (VA28/99), and {\it Iberdrola, S.A}.
M.A.C. acknowledges a scholarship granted from the Basque Government, and
J.S.D. a grant from The University of Basque Country/{\it Euskal Herriko Unibertsitatea}.
M.A.C. and J.S.D. wish to acknowledge the kind hospitality of the {\it Departamento de
F\'{\i}sica Te\'orica} in Valladolid.


\section*{Abstract}

Separated flows past complex geometries are modelled by discrete
vortex techniques. The flows are assumed to be rotational and inviscid, and a new technique is described to determine the streamfunctions for linear shear profiles. The geometries considered are the snow cornice and the backward-facing step, whose edges allow for the separation of the flow and reattachment downstream of the recirculation regions. A point vortex has been added to the flows in order to constrain the separation points to be located at the edges, while the conformal mappings have been modified in order to smooth the sharp edges and let the separation points be free to oscillate around the points of maximum curvature. Unsteadiness is imposed on the flow by perturbing the vortex location, either by displacing the vortex from equilibrium, or by imposing a random perturbation with zero mean on the vortex in equilibrium. The trajectories of passive scalars continuously released upwind of the separation point and trapped by the recirculating bubble are numerically integrated, and concentration time series are calculated at fixed locations downwind of the reattachment points. This model proves to be capable of reproducing the trapping and intermittent release of scalars, in agreement with the simulation of the flow past a snow cornice performed by a discrete multi-vortex model, as well as with direct numerical simulations of the flow past a backward-facing step. The simulation results indicate that for flows undergoing separation and reattachment the unsteadiness of the recirculating bubble is the main mechanism responsible for the intense large-scale concentration fluctuations downstream.

\section{Introduction}

Transport and dispersion of passive scalars in turbulent boundary layers over complex geometries is an important problem of practical relevance, which has been widely investigated in recent years mainly because of a growing concern about the environmental issue, but also because of its engineering applications such as optimization of combustors or evaluation of mixing properties of flows in chemical reactors.

However, the complexity of flows in which separation and reattachment occur make them less suitable to analytical studies than regular boundary layer flows. Relatively little is known about the behaviour of a plume near hills, or buildings, or, in general, over irregular terrain (Hosker 1984). For a classical problem like the flow over a backward-facing step, there is not an applicable mathematical technique for modeling the evolution of the scalar field. Laboratory experiments generally provide an Eulerian portrait of this flow (e.g. Bradshaw \& Wong 1972), but with no insight into the dispersion process. More recently, Le, Moin \& Kim (1997) performed a direct numerical simulation of a backward-facing step for Reynolds numbers as high as 5100 revealing, among the others, some of the unsteady flow characteristics. They detected the oscillation of the (spanwise-averaged) reattachment location, in agreement with the results of several other works (e.g. Driver, Seegmiller \& Marvin 1983, 1987). The shear layer composed of many small high-intensity vortices, which extends to the reattachment point, rolls up forming a large-scale structure in the recirculation region: the periodic detachment of this large-scale vortex from the step causes the oscillation of the reattachment point.
The main aim of this paper is to model the unsteadiness of reattached flows, and to assess how it affects the dispersion of passive scalars released by a continuous point source.

We will model separated flows by means of point vortex techniques,
assuming the flow as two-dimensional and inviscid. We will consider two different geometries: a backward-facing step, and a ``snow cornice''. Snow cornices are natural devices that control the flow separation on mountain crests by trapping vortices. The modeling technique is not new: Ringleb (1961) studied the steady separated flow past a snow cornice assuming a two dimensional potential flow, which allowed him to solve the problem analytically by means of the classical method of conformal mapping. To represent a snow cornice, Ringleb devised mappings to transform the real axis in the complex $\zeta$-plane onto a line forming a sharp edge in the complex $z$-plane. In this way the complex potential of a uniform flow over the real axis in the $\zeta$-plane was ascribed to the flow on the corresponding region in the $z$-plane. A point vortex (together with its reflected image) was added to the flow to model the recirculating region. The vortex was in equilibrium and satisfied a steady Kutta condition. In a more recent paper, Cortelezzi, Leonard \& Doyle (1994) apply point vortex techniques to model the unsteady separated flow over a semi-infinite plate, and to assess active control techniques.

Because of the incompressibility condition, the Lagrangian motion of a passive scalar is governed by a Hamiltonian system whose conjugate variables are the Cartesian coordinates of the scalar, and whose Hamilton function is the streamfunction. The motion of a point vortex is governed by a Hamiltonian system as well, where the Hamilton function is related to the Green function for the Laplacian of the flow region (Masotti 1931; Lin 1941).

In this respect, the motion of the vortex in the Ringleb model is
integrable since it is governed by a one-degree of freedom autonomous
Hamiltonian system. Incidentally, although we found that Ringleb made
a formal error in determining the law of the motion of a vortex in
the presence of a wall, and therefore failed to derive the correct
equilibrium conditions, his results are qualitatively correct (see Appendix B).

The streamfunction, namely the Hamiltonian of the fluid particles,  depends on the vortex position and therefore is in general time-dependent (non-autonomous). For such a system we expect a non integrable, chaotic, particle motion (Novikov \& Sedov 1979; Aref 1983). However, if the vortex is located at its equilibrium point, the Hamiltonian of the particles is autonomous and the motion is integrable and regular. For the integrable case, the flow field is virtually divided in two main bodies: the fluid entrained by the vortex (closed streamlines encircling the vortex), and the free-stream flow (open streamlines). The two regions do not exchange fluid with each other and are separated by a streamline which leaves the sharp edge of the wall, and reattaches downstream (i.e. a heteroclinic orbit). For instance, in the transformed $\zeta$-plane this flow is a vortex pair, where the real axis is the line of symmetry.

The present model is different from Ringleb's in three main aspects:

(i) We assume a shear flow as the asymptotic upstream condition in order to have a more meaningful flow than the potential one, but still simple enough to allow for an analytical study.

(ii) The system has been made unsteady by perturbing the equilibrium of the vortex.

(iii) Finally, we consider geometries where the sharp corners are
smoothed to blunted edges with large but finite curvature, to satisfy
an unsteady Kutta condition for a vortex of constant circulation not
in equilibrium.

The large curvature edge still causes separation, in that according to the Kutta condition it is possible to define a vortex whose strength produces separation. A vortex slightly displaced from its stable equilibrium moves along a small periodic orbit, and causes the oscillation of the detachment point around the maximum curvature point. Since the edge is smooth, the flow around the edge does not have any singularities and the Kutta condition can be considered fulfilled.

A typical representation of the mean flow vertical profile by power laws is
used in the atmospheric boundary layer in stable conditions. For instance, a linear shear profile has been considered by Hunt, Leibovich \& Richards (1988) to study how the shear affects the flow past hills with low slopes.

On the transformed $\zeta$-plane, the flow reduces to the oscillating vortex pair (OVP). Rom-Kedar, Leonard \& Wiggins (1990) studied the OVP flow extensively, and elucidated the mechanisms of the phase-space Lagrangian chaos. The stable and unstable manifolds of the heteroclinic orbit intersect each other, and drive the fluid to be entrained and detrained by the vortex. The Eulerian counterpart of this phenomenon is the chaotic mixing of the fluid in the region surrounding the separatrix streamline. We are interested in the effects of these phenomena on the dispersion of scalars.

We investigated how the separation and reattachment of an unsteady flow influence the behaviour of effluent plumes released near or inside the recirculating region. We found that the presence of the recirculation bubble dramatically affects the large-scale dynamics of passive scalars, causing important fluctuations of the concentration field downstream.

In the following sections the Ringleb model for steady potential flow is recalled, then a method is proposed to study analytically unsteady rotational flows with constant vorticity. Dynamical systems formalism is used to find equilibrium configurations, to discuss their stability and to recognize and quantify chaotic mixing. Numerical simulations for two different geometries are presented, along with the concentration time series at several sampling points.

\section{Mathematical formulation of the steady potential flow}

In this section we describe a mathematical model of recirculating
flow over complex geometries such as a snow cornice and a backward-facing step. This model, which was proposed by Ringleb (1961), is able to represent a two-dimensional steady separation from the edge. It constitutes the basis for the model of unsteady, recirculating shear flow, which we present in \S 3.

We consider a region bounded by a piecewise analytical curve in the $z$-plane, where $z=x+iy$. The solid boundary extends to infinity along the $x$-axis, and is characterized by a sharp corner. As a consequence of the Riemann mapping theorem (see, e.g., Nehari 1975; Henrici 1974), such a region can be obtained as the conformal image of the half-plane $\eta\ge 0$ in the $\zeta$-plane, where $\zeta=\xi+i\eta$. The infinity in the $z$-plane corresponds to the infinity in the $\zeta$-plane. For instance, the function suggested by Ringleb (1961),
\be
z=\zeta+\frac{\zeta_1^2}{\zeta-\zeta_1}\label{a1}
\ee
where $\zeta_1=\xi_1+i\eta_1$ is a complex constant with $\eta_1<0$, yields to the geometry plotted in figure 1a. The backward-facing step shown in figure 1b can be obtained by the Schwartz-Christoffel mapping:
\be
z=\sqrt{\zeta^2-1}+\log \left(\zeta+\sqrt{\zeta^2-1}\right)\;.\label{a2}
\ee

The model consists of a vortex of constant strength which satisfies a Kutta condition, in equilibrium with a steady free-stream velocity.
The complex potential of the flow is built by superimposition of basic flows: in the $\zeta$-plane it is expressed by the function:
\be
w=q_\infty\zeta+\frac{\gamma}{2\pi i}\log\left(\frac{\zeta-\zeta_0}{\zeta-\zeta_0^\star}\right)\label{a3}
\ee
where $q_\infty$ is the free-stream velocity, and the second term on the right hand side is the complex potential of a point vortex of strength $\gamma$ located at $\zeta_0$ and of its image at the complex conjugate position $\zeta_0^\star$. Equation (\ref{a3}) represents the complex potential of a vortex pair of opposite signs in the $\zeta$-plane (see figure 7).

Since the complex potential is invariant under a conformal mapping $z=z(\zeta)$, the constant $q_\infty$ is equal to the free-stream velocity $Q_\infty$ in the $z$-plane:
\be
Q_\infty=\lim_{z\rightarrow\infty} \frac{dw}{dz}=\lim_{\zeta\rightarrow\infty}\frac{dw}{d\zeta} \left(\frac{dz}{d\zeta}\right)^{-1}=q_\infty\;.
\label{a4}
\ee

The trajectories of a vortex in the physical $z$-plane are the solutions of the Hamiltonian system:
\be
\dot x_0=\frac{\prz H}{\prz y_0};\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \dot y_0=-\frac{\prz H}{\prz x_0}\label{a6}
\ee
where the Hamiltonian $H$ can be derived from the Hamiltonian $H'$ of a vortex in the $\zeta$-plane, according to the Routh rule (e.g. Clements 1973):
\be
H=H'+\frac{\gamma}{4\pi}\log\left|\frac{d z_0}{d \zeta_0}\right|\label{a7}
\ee
where:
\be
H'=q_\infty\eta_0+\frac{\gamma}{4\pi}\log{\eta_0}\;.
\label{a8c}
\ee
Derivations of the Hamiltonian $H'$ and the Routh transformation
rule are given in Appendix A.
\subsection{Steady Kutta condition}
The separation of a flow around a sharp or a large curvature corner is captured by the model if it satisfies the Kutta condition, that is, a relationship between the strength of the vortex and its position in order to eliminate the singularity of the flow. Therefore, the condition of zero velocity at the edge is imposed:
\be
\left(\frac{dw}{d\zeta}\right)_{\zeta=\zeta_c}=0
\label{b1}
\ee
where $\zeta_c$ indicates the location of the corner in the $\zeta$-plane.

For the complex potential $w(\zeta)$  (\ref{a3}), the Kutta condition is satisfied when
\be
\frac{\gamma}{q_\infty}=-\frac{\pi|\zeta_c-\zeta_0|^2}{\eta_0}
\label{b2}
\ee
where we notice that there are infinite possible values for the ratio $\gamma/q_\infty$. For the mapping (\ref{a1}) $\zeta_c=0$, whereas $\zeta_c=-1$ for (\ref{a2}).

On the physical plane, the complex velocity at the corner $z_c$ is
\be
\lim_{z\rightarrow z_c}\frac{dw}{dz}=\lim_{\zeta\rightarrow \zeta_c}\frac{dw}{d\zeta}\frac{d\zeta}{dz}
\label{b3}
\ee
whose value is finite when the corner is a cusp, as in the case of mapping (\ref{a1}), while it is null in the other case.

In order to complete the modeling of the recirculating flow, we are
now concerned with finding the equilibrium conditions for the vortex. We recall that the phase space of the dynamical systems (\ref{a6}) coincides with the flow region, in the sense that the Cartesian coordinates of the vortex ($x_0$, $y_0$) represent the conjugate variables, and the level lines of the Hamiltonian $H$ form the pattern of the possible trajectories of a vortex of assigned strength. Therefore the equilibrium locations ($\dot x_0=0$, $\dot y_0=0$) are the fixed points of this map: the elliptic points correspond to stable equilibrium, the hyperbolic ones to unstable equilibrium. The model of a steady recirculating flow behind a corner is represented by a vortex in stable equilibrium (i.e.  a {\em trapped vortex}), and satisfying the Kutta condition.
The solutions to the equilibrium equations, if any, depend on the geometry of the solid boundary and on the ratio of the vortex intensity to the asymptotic velocity, i.e. $\gamma/q_\infty$. Since the geometries we considered were found to allow vortex capturing, we expect a set of stable solutions.

Therefore, the vortex coordinates ($\xi_0,\eta_0$) on the $\zeta$-plane have to satisfy simultaneously the equilibrium condition and the Kutta condition (\ref{b2}). They are the solutions of the equation:
\be
\frac{|\zeta_c-\zeta_0|^2}{4\eta_0}\left\lbrack\frac{1}{\eta_0}\frac{d\zeta_0}{dz_0}-i\frac{d}{dz}\left(\log
\frac{d\zeta_0}{dz_0}\right)\right\rbrack=1\;.
\label{b4}
\ee
The vortex strength is then obtained by the Kutta condition (\ref{b2}).
\subsection{Comparison with a multiple vortex method}
The accuracy of the above vortex model has been checked by comparing some of the results with those generated by a more sophisticated vortex method. We used a discrete vortex method (Ceschini 1993; Ferlauto 1996) to describe the impulsive start of the flow characterized by vortex shedding. We performed the simulation for the snow-cornice geometry. The method consists of assuming the initial flow around the wall as potential, and satisfying an unsteady Kutta condition by adding a free point vortex of suitable strength in a prescribed location close to the corner. The vortex moves because of, and according to, the wall and the background flow. The motion is then numerically integrated. At given time intervals a new vortex is introduced into the flow to restore the regularity. Figure 2 shows the pattern of the vortices shed throughout the transient, where the vortices strengths are proportional to the symbol sizes. The transient is characterized by a vortex sheet rolling up around the point of highest absolute circulation. This result is consistent with the description of the transient as given by Cortelezzi, Leonard \& Doyle (1994) and Le, Moin \& Kim (1997) for other vortex trapping flows. The vortex sheet increases in size and circulation until eventually reaching a more steady configuration, with small oscillation of its global circulation and centroid.

We compared the asymptotic circulation and the centroid location with
the strength and location of the single point vortex used in our
model (the centroid location is evaluated assuming the strengths as
masses). We plot in figure 3a the total circulation as a function of
time (solid line), as well as the circulation for the single vortex
model (dotted line). In figure 3b we plot the distance of the centroid for
the multiple vortex method from the equilibrium point of a single
vortex. The two models are in good agreement with each other as far
as the intensity of the vorticity field and the location of the
centre of vorticity are concerned, showing that the large-scale characteristics of the flow are well captured by a single vortex model. Of course, the strength of the single vortex model is in its simplicity, which permits the analytical description of a large-scale mechanism of dispersion, as discussed in \S 4.

\section{Shear flow}
In this section we describe a method to model the separation and recirculation of a shear flow over complex geometries.

We assume a rotational flow field, with vorticity $\nabla \times \bq=\omega=const$ and with the following velocity profile in the far field:
\be
\lim_{x\rightarrow\pm\infty}q=u_\infty-\omega y\;.\label{c1}
\ee
The motion of passive scalars is governed by the Hamiltonian dynamical system:
\be
\dot x=\frac{\prz \psi}{\prz y};\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \dot y=-\frac{\prz \psi}{\prz x}\label{a5}
\ee
where the streamfunction $\psi$ represents the Hamilton function and the scalar coordinates ($x,y$) are canonical variables.
The streamfunction $\psi$ is related to the vorticity $\omega$ by the linear Poisson equation:
\be
\nabla^2\psi=-\omega\;.\label{c2}
\ee
To solve the Poisson equation (\ref{c2}) we developed an analytical procedure based on conformal mapping, which is inspired by the method described by Tsien (1943).\\
Let us assume
\be
\psi=\psi_\omega+\psi_p\label{c2a}
\ee
where the streamfunction $\psi_\omega$ is relevant to the shear flow $u=-\omega y$, i.e.
\be
\psi_\omega=-\frac{1}{2}\omega\; y^2\label{c3}
\ee
and the function $\psi_p$ is harmonic in the flow region.\\
Since
\be
\nabla^2\psi_\omega=-\omega\label{c4}
\ee
and
\be
\nabla^2\psi_p=0\label{c5}
\ee
equation (\ref{c2}) is satisfied. However, the problem has been
reduced to the determination of the harmonic streamfunction $\psi_p$
which ensures the fulfillment of the boundary conditions, that is,
the far field velocity equal to $(u_\infty-\omega\,y)$ and the
impermeability of the wall.

Therefore, we seek a complex potential $w_p$ whose imaginary part is
$\psi_p$. We define the flow regions in the $z$-plane as the conformal representation of the half-plane $\Im(\zeta)>0$ according to the mappings:
\be
z=\zeta+\frac{\zeta_1^2}{\zeta+i\delta -\zeta_1}\label{c6}
\ee
and
\be
z=\sqrt{(\zeta+i\delta)^2-1}+\log\left(\zeta+i\delta+\sqrt{(\zeta+i\delta)^2-1}\right)\label{c7}
\ee
which are the modified Ringleb mapping (\ref{a1}) and the Schwartz-Cristoffel mapping (\ref{a2}), respectively, where we introduced the real positive parameter $\delta$ in order to smooth the corners (the smaller $\delta$, the sharper the corner).

The further mapping:
\be
\lambda=\frac{i-\zeta}{i+\zeta}\label{c7a}
\ee
transforms the half-plane $\Im(\zeta)>0$ onto the interior of the
unit circle $|\lambda|=1$ in the $\lambda$-plane. Summing up, the chain mapping $z\rightarrow\zeta\rightarrow\lambda$ maps the physical boundary of the flow region onto the unit circle in the $\lambda$-plane, where $\lambda=-1$ corresponds to infinity in the $z$-plane.

On the physical plane, the shear flow (\ref{c3}) induces the following velocity component $\tilde{v}$ normal to the boundary:
\be
\tilde{v}=\omega\; y_b\; \sin \beta\label{c7b}
\ee
where the subscript $b$ denotes boundary, and $\beta=\arg\left(dz/d\zeta\right)$ is the angle between the boundary and the $x$-axis.

According to equation (\ref{c7b}), the condition of impermeability for the boundary in the $z$-plane is satisfied when the velocity component normal to the boundary, due to the complex potential $w_p$, is equal to $-\tilde{v}$, that is:
\be
-\Im\left(\frac{dw_p}{dz}e^{i\beta}\right)=-\tilde{v}\label{c9}
\ee
which translates on the $\lambda$-plane as
\be
\Re\left(\frac{dw_p}{d\lambda}\lambda\right)_{|\lambda|=1}=-\tilde{v}\left|\frac{dz}{d\lambda}\right|
_{|\lambda|=1}\;.
\label{c10}
\ee

As there are no point vortices (they can be added to the flow later), the complex velocity $dw_p/dz$ cannot be singular inside the flow field. Since the mappings $z\rightarrow\lambda$ are regular inside the unit circle $|\lambda|\le 1$ ($\lambda\ne -1$) in the $\lambda$-plane, then the complex function $dw_p/d\lambda$ must not be singular inside the unit circle as well, except for the point $\lambda=-1$ that corresponds to $z=\infty$. Therefore we can assume:
\be
\frac{dw_p}{d\lambda}=F(\lambda)+\sum_{n=1}^\infty(a_n-ib_n)\lambda^{n-1}\label{c11}
\ee
where $F(\lambda)$ contains the possible singularities located in $\lambda=-1$, that is, in general,
\be
F(\lambda)=\sum_{j=1}^\infty c_j (\lambda+1)^{-j}\;.
\label{c12}
\ee

To determine the $c_j$ coefficients, we recall that the velocity at infinity ($z=\infty$), due to the complex potential $w_p$, is $u_\infty$. Since for the mappings (\ref{c6}) or (\ref{c7}) the derivative at infinity is
\be
\lim_{z\rightarrow\infty}\frac{d\zeta}{dz}=1\label{c13}
\ee
we have
\be
u_\infty=\lim_{\zeta\rightarrow\infty}\frac{dw_p}{d\zeta}=
\lim_{\lambda\rightarrow -1}\left(\frac{dw_p}{d\lambda}\frac{d\lambda}{d\zeta}\right)
=\lim_{\lambda\rightarrow -1}\left[\frac{dw_p}{d\lambda}\;\frac{(\lambda+1)^2}{-2i}\right]
\label{c14}
\ee
which implies that all the coefficients $c_j$ for $j>2$ are null, and that
\be
c_2=-2i\;u_\infty\;.
\label{c15a}
\ee

We determine the coefficient $c_1$ by considering that the mass flow $Q$ across the boundary, due to the shear flow, has to be balanced by the potential flow. We have
\be
Q=\int_{-\infty}^{+\infty}-\omega\; y_b\;\frac{dy_b}{dx_b}\;dx_b=-\frac{\omega}{2}\;\left[ y_b^2\right]_{-\infty}^{+\infty}
\label{c15}
\ee
and also, according to the invariance property of the complex potential,
\be
-Q=-i\;\oint_{|\lambda|=1} \frac{dw_p}{d\lambda}= \pi\;c_1\;.\label{c16}
\ee
The function $F(\lambda)$ is then fully determined:
\be
F(\lambda)=-2i\;u_\infty\;\frac{1}{(\lambda+1)^2}-\frac{Q}{\pi}\;\frac{1}{\lambda+1}
\label{c17}
\ee
According to equation (\ref{c15}), the mass flow $Q$ relevant to the mapping (\ref{c6}) is null, while for the mapping (\ref{c7}) it is
\[Q=-\frac{\omega}{2}(\pi^2+2\pi\;\delta)\;.\]

To calculate the second term on the right hand side of equation (\ref{c11}), we combine equation (\ref{c11}) with equations (\ref{c10}) and (\ref{c17}) evaluated at $|\lambda|=1$, that is $\lambda=\exp(i\vth)$:
\begin{eqnarray}
\sum_{n=1}^\infty(a_n\cos n\vth+b_n\sin n\vth)&=&-\tilde{v}\left|\frac{dz}{d\lambda}\right|
-\Re\left(-2i\;u_\infty\;\frac{\lambda}{(\lambda+1)^2}-\frac{Q}{\pi}\;\frac{\lambda}{\lambda+1}\right)\nonumber\\
 &=& -\tilde{v}\left|\frac{dz}{d\lambda}\right|+\frac{Q}{2\pi}\;.
\label{c18}
\end{eqnarray}
Since the right hand side of (\ref{c18}) is not singular, the coefficients $a_n$ and $b_n$ on the left hand side tend to zero for large $n$. The Fourier series on the left hand side can be suitably truncated and the coefficients $a_n$ and $b_n$ evaluated numerically. The FFT algorithms are efficient tools for this purpose.

The complex potential $w_p$ is finally obtained by integrating equation (\ref{c11}):
\be
w_p=\frac{2i\;u_\infty}{\lambda+1}-\frac{Q}{\pi}\log(\lambda+1)+\sum_{n=1}^\infty\frac{a_n-ib_n}{n}\;\lambda^n\;.
\label{c19}
\ee

In general, for geometries defined by the mappings (\ref{c6}) and
(\ref{c7}), the flow does not satisfy the Kutta condition, i.e. it does not separate at the corner. The rotational flow field does exhibit a recirculating bubble, but the separation point is unphysically located downstream from the high curvature corner. To enforce the Kutta condition, a vortex can be added to the flow, as shown in \S 2 for a potential flow: the same considerations regarding its strength and equilibrium hold for a shear flow.

By adding a point vortex to the flow, the complex potential becomes
\be
w_p=\frac{2i\;u_\infty}{\lambda+1}-\frac{Q}{\pi}\log(\lambda+1)+\sum_{n=1}^\infty\frac{a_n-ib_n}{n}\;\lambda^n
+\frac{\gamma}{2i\;\pi}\log\left(\frac{\zeta-\zeta_0}{\zeta-\zeta^\star_0}\right)\;.
\label{c20}
\ee

The complete streamfunction, relevant to the sum of the shear and the potential component of the flow, is
\be
\psi=-\frac{1}{2}\omega\;y^2+\Im(w_p)\;,\label{c20a}
\ee
which provides the following flow velocity:
\be
\dot z^\star=-\omega\;y+\left\lbrack F(\lambda)+\sum_{n=1}^\infty(a_n-ib_n)\lambda^{n-1}\right\rbrack
\frac{d\lambda}{dz}+\left\lbrack\frac{\gamma}{2i\;\pi}\left(\frac{1}{\zeta-\zeta_0}-\frac{1}{\zeta-\zeta_0^\star}
\right)\right\rbrack\frac{d\zeta}{dz}\;.
\label{c21}
\ee
The vortex velocity can be obtained as
\be
\dot z_0^\star=-\omega\;y_0+\lim_{z\rightarrow z_0}\left(\frac{dw_p}{dz}-\frac{\gamma}{2\pi i}\;\frac{1}{z-z_0}\right)
\label{c22}
\ee
and is governed by a Hamiltonian system as in the case of the
potential flow discussed in \S 2:
\be
\dot x_0=\frac{\prz H}{\prz y_0}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \dot y_0=-\frac{\prz H}{\prz x_0}\label{c23}
\ee
where the Hamilton function $H$ is:
\be
H=-\frac{1}{2}\;\omega\;y^2+\Im\left\lbrack\frac{2i\;u_\infty}{\lambda+1}-\frac{Q}{\pi}\log(\lambda+1)+\sum_{n=1}^\infty\frac{a_n-ib_n}{n}\;\lambda^n\right\rbrack-\frac{\gamma}{4\pi}\left\lbrack\log \eta_0-\log\left|\frac{dz}{d\zeta}\right|\right\rbrack\;.
\label{c24}
\ee

The strength and location of the vortex is then determined according to equations (\ref{c22}) and (\ref{c23}), imposing the flow at rest on the corner $z_c$ and the vortex at its equilibrium location ($\dot x_0=\dot y_0=0$). Figure 4 shows streamlines and velocity vectors for a separated shear flow for a geometry given by the map (\ref{c6}); the flow field is characterized by $\omega/u_\infty=-10$ and $\gamma/u_\infty=-2.243$.

The stability of the equilibrium can be inferred from figure 5, which displays the contour levels of the Hamiltonian of the vortex, namely the possible trajectories of the vortex. The equilibrium point, marked with a square, coincide with an elliptic, i.e. stable, fixed point in the phase space.

\subsection{The unsteady flow}

We examine the Lagrangian transport of passive scalars in the case of unsteady flow. Recalling equations (\ref{a5}) and (\ref{c20a}), the motion is governed by the system of equations:
\be
\dot x=\frac{\prz}{\prz y} \psi(x,y;x_0,y_0);\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \dot y=-\frac{\prz}{\prz x} \psi(x,y;x_0,y_0)\;,\label{g1}
\ee
where the streamfunction $\psi$ depends on both the particle coordinates ($x,y$) and the vortex location $(x_0,y_0)$. For a vortex located at its equilibrium position ($x_0^E,y_0^E$), the one-degree of freedom Hamiltonian system (\ref{g1}) is autonomous, and hence integrable. In this case, the phase portrait coincides with the streamline pattern for the steady flow as shown in figure 4. The phase space is characterized by a streamline connecting two stagnation points, or according to dynamical system jargon, a heteroclinic orbit connecting two hyperbolic fixed points. Figure 6a shows this streamline in  the physical $z$-plane. Such a feature has great relevance because it is capable of triggering chaotic mixing if the system undergoes a time-dependent perturbation, as shown below.

If the vortex is displaced from its equilibrium location ($x_0^E,y_0^E$) by a small quantity $\varepsilon$, its trajectory follows a closed orbit, as shown in figure 5. Its periodic motion obeys equations of the form:
\be
x_0=x_0^E+f(x_0^E, y_0^E,\varepsilon,t),\;\;\;\;\;\;\;\;\;\;\;y_0=y_0^E+g(x_0^E, y_0^E, \varepsilon,t)\;.\label{g2}
\ee
In this case, system (\ref{g1}) is no longer autonomous, because now the streamfunction (\ref{c20a}) depends on time through the time-dependent vortex coordinates $(x_0(t),y_0(t))$. Formally, we can write
\begin{eqnarray}
\psi&=&\psi(x,y;x_0^E,y_0^E,\varepsilon)\nonumber\\
 &=&\psi^E(x,y;x_0^E,y_0^E)+\varepsilon\left(\frac{\prz\psi}{\prz x_0}\frac{\prz f}{\prz\varepsilon}
+\frac{\prz\psi}{\prz y_0}\frac{\prz g}{\prz\varepsilon}\right)+O(\varepsilon^2)\;,\label{g3}
\end{eqnarray}
and the governing equations (\ref{g1}) reduce to the time-dependent perturbation of an integrable Hamiltonian system:
\be
\label{g4}
\begin{array}{lcc}
\dot x&=& \displaystyle\frac{\prz\psi^E}{\prz y}+\varepsilon\displaystyle\frac{\prz}{\prz y} \left(\displaystyle\frac{\prz\psi}{\prz x_0}\displaystyle\frac{\prz f}{\prz\varepsilon}
+\displaystyle\frac{\prz\psi}{\prz y_0}\displaystyle\frac{\prz g}{\prz\varepsilon}\right)+O(\varepsilon^2)\nonumber\\
\dot y&=& -\displaystyle\frac{\prz\psi^E}{\prz x}-\varepsilon\displaystyle\frac{\prz}{\prz x} \left(\displaystyle\frac{\prz\psi}{\prz x_0}\displaystyle\frac{\prz f}{\prz\varepsilon}
+\displaystyle\frac{\prz\psi}{\prz y_0}\displaystyle\frac{\prz g}{\prz\varepsilon}\right)+O(\varepsilon^2)\;.
\end{array}
\ee

As shown in figure 6a, for $\varepsilon=0$ the unstable manifold leaving the hyperbolic fixed point $P^-$ joins smoothly the stable manifold going into the other fixed point $P^+$, and the heteroclinic orbit of the unperturbed system is a single line connecting $P^-$ and $P^+$. However, it is well known that $\varepsilon\ne 0$ dramatically affects the phase portrait in that the time dependence of the system, if not associated to any kind of symmetry, causes non-integrability (Novikov \& Sedov 1979; Aref 1983). The stable and unstable manifolds do not join smoothly but intersect each other transversally forming a paradigmatically chaotic {\em tangle} (e.g. Guckenheimer \& Holmes 1983, Tabor 1989, Wiggins 1992). Since any intersection point maps into another intersection point, the number of transverse intersections is infinite, whereas the distance between intersection points tends to zero as they approach a hyperbolic fixed point.  At the same time, the system preserves the area of the lobes embraced by the arcs of the manifolds clipped between two intersections.
Figure 6b shows a Poincar\'e section of the flow over the snow cornice, enlightening the scenario of the stable and unstable manifold intersecting transversally.

The chaotic character of a system can be deduced whenever evidence of such transverse crossings is provided. For instance, several authors used the Melnikov technique (Melnikov 1963) to detect analytically the transverse crossings between stable and unstable manifolds for (time-periodically) perturbed integrable systems (e.g. Rom-Kedar, Leonard \& Wiggins 1990, Zannetti \& Franzese 1994, Del Castillo-Negrete 1998 and references therein).

\section{Results}

The representation on the transformed $\zeta$-plane of a steady shearless flow over a step or a snow cornice is similar in many respects to a vortex pair flow, except for the differences due to the Routh correction (\ref{a7}), which affects the vortex velocities. The streamline pattern of a vortex pair flow in the frame moving with the vortices is shown in figure 7, where we can see the two hyperbolic fixed points ($P^-,P^+$) connected by heteroclinic orbits. Rom-Kedar, Leonard \& Wiggins (1990) perturbed the vortex pair by an external periodic strain-rate field, and studied extensively the chaotic dynamics of the resulting Oscillating Vortex Pair flow (OVP). They focused on the dynamics of the lobes enclosed by the perturbed heteroclinic orbits. Rom-Kedar, Leonard \& Wiggins (1990), and Wiggins (1992), provided valuable tools for a quantitative evaluation of the mixing process in terms of the amount of fluid involved and the residence time of the fluid particles in the chaotic region.

The unsteady flows over complex geometries that we described in this paper show the same basic mechanisms of mixing as for the OVP flow.
Besides the fluid permanently entrained by the vortex, and the fluid flowing unperturbed downstream, the scenario is enriched by fluid particles temporarily trapped by the vortex and then detrained, before flowing downstream. This phenomenon, known as {\em transient chaos}, was acknowledged and investigated a few years ago by Pentek, Tel \& Toroczkai (1995) in the context of the advection problem of passive tracers in the velocity field of vortex pairs.

We focus on the phenomenon of heteroclinic chaos as a large scale mechanism of intermittent release, causing intense concentration fluctuation of passive tracers. The entangling of the heteroclinic orbits is a consequence of the unsteadiness of the recirculating region and, in particular, of the oscillatory motion of the reattachment location. This oscillation has been observed in large eddy simulations as well as in two and three-dimensional direct numerical simulation of high Reynolds number flows over a backward-facing step (Driver, Seegmiller \& Marvin 1983, 1987; Le, Moin \& Kim 1997 and references therein). The periodic movement of the reattachment location is caused by the formation and detachment from the step of large-scale coherent vortical structures. The discrete vortex method that we applied to the snow cornice reproduces the same phenomena: figure 3a shows the fluctuations of the total circulation during the transient, due to the periodic detachment of clusters of vortices from the recirculating region. The phenomenon is evident during the transient because no perturbations were introduced in this simulation, but a natural oscillation of the centre of vorticity takes place during the transient, as shown in figure 3b. The transient is characterized by the roll-up of the vortex sheet which is shed into the fluid (see figure 2), consistently with other discrete vortex models of reattached flows (Cortelezzi, Leonard \& Doyle 1994) and DNS of the backward-facing step flow (Le, Moin \& Kim 1997).

As shown below, the unsteady single vortex model  is able to display
the periodic detachment of coherent structures, even though these structures are not vortices but fluid particles trapped into the lobes of perturbed heteroclinic orbits.

We present the results of two simulations showing the intermittent release of tracers. The first refers to the snow cornice generated by the modified Ringleb transformation (\ref{c6}). We perturbed the flow by displacing the vortex by $\varepsilon=0.005$ and letting it move freely along its periodic orbit. Passive tracers are released at a constant rate from a point upstream of the edge and close to the wall; their time evolution is described by integrating the governing equations (\ref{g1}) by means of a fourth-order Runge-Kutta algorithm.

The time evolution of the plume is obtained by recording the position of the tracers at fixed time intervals. Also, we record the scalar concentration, averaged over a small control volume, as a function of time.
Figure 8 shows the distribution of scalars in the field and their
concentration in a sampling volume, marked by a shaded rectangular
area, at two different times. The sampling volume is  downstream. The periodic character of the concentration is demonstrated by its time history, plotted in the upper portion of each frame. The wave form of the concentration was found to be dependent on the location of the monitoring point.

A second example of intermittent release is given by the simulation performed for the geometry obtained by the modified Schwartz-Christoffel mapping (\ref{c7}), i.e. the backward-facing step with smoothed corner. The unsteadiness of the flow is simulated by imposing a random perturbation with zero mean on the position of the vortex in equilibrium. Again, the concentration at a point downstream as a function of time in figure 9 illustrates the intermittent formation of lobes with high density of tracers which are responsible for the concentration fluctuations. Because of the random nature of the perturbation, the concentration is not periodic in time.

The strongly evident feature of both the simulations is that the perturbation affects the detachment region slightly but the reattachment region quite strongly: the upstream stagnation point moves imperceptibly while the downstream stagnation point has a patently larger movement. This is consistent with the behaviour of real flows around edges, where viscosity effects constrain the separation to be located at the edge, making the separation region rather insensitive to perturbations. From the analytical point of view the two stagnation points in the $\zeta$-plane have oscillations of comparable amplitude, however the mapping onto the physical $z$-plane alters the displacements. The lengths in the $z$-plane are proportional, to first order, to the absolute value of the mapping derivative $|dz/d\zeta|$, whose minimum value is at the edge (for a sharp convex edge $dz/d\zeta=0$). This analytical property provides a reasonable extension of the above mechanism of scalar intermittent release to the geometries characterized by high curvature edges, as well as by separation and reattachment of the flow.

\section{Conclusions}
This study has examined the effects of the unsteadiness of reattached flows on the large-scale dynamics of dispersion of passive scalars. We considered the two cases of flow past a snow cornice and a backward-facing step, and we used point-vortex models to simulate unsteady separated shear flows. A new technique has been described to determine the streamfunction in a field of constant vorticity.

A rotational flow over the above geometries reproduces the separation, even though the location of the separation point must be corrected by adding a point vortex to the flow, in order to have separation at the point of maximum curvature. The classical Ringleb mapping for snow cornices, and a Schwartz-Christoffel mapping for step flows have been modified by smoothing the corners. Such modification is necessary to introduce the unsteadiness of the recirculating bubble, in that the separation point is allowed to move around the point of maximum curvature, while having the flow still satisfy an unsteady Kutta condition.

The simulation performed for the snow cornice by a more sophisticated multiple-vortex model detected the oscillation of the center of vorticity, thus showing that the oscillatory motion of the reattachment point is associated with the unsteadiness of the recirculating bubble. The same kind of unsteadiness was imposed in the simulation by the simpler single-vortex model, by perturbing the equilibrium position of the vortex. As the multiple-vortex simulation shows, the flow is characterized by a periodic release of clusters of vortices from the recirculating region behind the crest, in agreement with the direct numerical simulations of the backward-facing step of Le, Moin \& Kim (1997). The same mechanism of releasing vorticity holds for passive tracers.

Two different perturbations were applied to the system: a periodic, natural perturbation obtained by displacing the vortex from its equilibrium location, and a random motion with zero mean imposed on the vortex in equilibrium. The scalar concentration averaged over small control volumes was calculated as a function of time. For the case of periodic vortex motion, our model detects the intermittent release of blobs of tracers with the same period as the movement of the reattachment location (i.e. the same period as the vortex motion). For the case of random oscillation of the recirculating bubble, the release of coherent structures of tracers was still observed, although the wave form of the concentration is not longer periodic.

Intermittent dispersion is associated with the phase-space lobe
dynamics (Rom-Kedar, Leonard \& Wiggins 1990, Wiggins 1992), which is
an intrinsically intermittent phenomenon. Since very small
perturbations are sufficient to start the above large-scale
phenomenon, we expect intense concentration fluctuations to appear
downstream whenever a flow undergoes separation and reattachment.

\bigskip

Much of this work was done while P. Franzese was a postdoctoral
fellow at CSIRO Atmospheric Researh at Aspendale, Australia. We express our appreciation to Dr. Mark Hibberd for his thorough comments,
as well as for his continuous support. We are also grateful to Dr.
Mike Borgas for many helpful discussions and for his encouragements.

\section*{Appendix A}
\subsection*{The Hamiltonian of a vortex}
The Green function $G(\xi,\eta;\xi_0,\eta_0)$ for a flow region $\cal R$ confined by a line $\cal C$ that extends to infinity is defined as
\[G(\xi,\eta;\xi_0,\eta_0)=g(\xi,\eta;\xi_0,\eta_0)+\frac{1}{2\pi}\log\sqrt{(\xi-\xi_0)^2+(\eta-\eta_0)^2}\]
where $g(\xi,\eta;\xi_0,\eta_0)$ is such that $G(\xi,\eta;\xi_0,\eta_0)=0$ on $\cal C$, and is harmonic with respect to $(\xi,\eta)$ on $\cal R$.
Moreover, the reciprocity property of Green functions, i.e. $G(\xi,\eta;\xi_0,\eta_0)=G(\xi_0,\eta_0;\xi,\eta)$, leads to the relationships:
\be
\label{a8a}
\begin{array}{lcl}
\displaystyle \frac{\prz}{\prz \xi_0}g(\xi_0,\eta_0;\xi_0,\eta_0)&=&2\displaystyle\lim_{\zeta\rightarrow \zeta_0}\displaystyle\frac{\prz}{\prz \xi}g(\xi,\eta;\xi_0,\eta_0)\\
\displaystyle \frac{\prz}{\prz \eta_0}g(\xi_0,\eta_0;\xi_0,\eta_0)&=&2\displaystyle\lim_{\zeta\rightarrow \zeta_0}\displaystyle\frac{\prz}{\prz \eta}g(\xi,\eta;\xi_0,\eta_0)\;.
\end{array}
\ee
The streamfunction $\psi(\xi,\eta;\xi_0,\eta_0)$ due to a free vortex located in $\zeta=\zeta_0$ must satisfy the boundary condition $\psi=const$ on $\cal C$, and has to be harmonic in the whole region except at $\zeta_0$, where it behaves as $(\gamma/2\pi)\log\sqrt{(\xi-\xi_0)^2+(\eta-\eta_0)^2}$.
Therefore, the streamfunction $\psi$ is related to the Green function by:
\be
\psi(\xi,\eta;\xi_0,\eta_0)=\psi_1(\xi,\eta)+\gamma G(\xi,\eta;\xi_0,\eta_0)\label{a7a}
\ee
where $\psi_1$ is the streamfunction of other flows possibly superimposed on the vortex.

A vortex does not induce velocity on itself. Therefore, it moves as a particle in the flow field from which the vortex induction has been removed. We can calculate the complex conjugate velocity of a vortex as:
\be
\dot \zeta_0^\star=\dot \xi_0-i\;\dot \eta_0=\lim_{\zeta\rightarrow \zeta_0}\left(\frac{dw}{d\zeta}-\frac{\gamma}{2\pi i}\frac{1}{\zeta-\zeta_0}\right)
=\left(\frac{\prz\psi_1}{\prz\eta}+i\frac{\prz\psi_1}{\prz\xi}\right)_{\zeta=\zeta_0}+\gamma \lim_{\zeta\rightarrow \zeta_0}\left(\frac{\prz g}{\prz \eta}+i\frac{\prz g}{\prz\xi}\right)
\label{a8}
\ee
where $w$ is the complex potential. According to equations (\ref{a8a}), the velocity of a vortex can be expressed as:
\be
\dot\zeta_0=\frac{\prz H'}{\prz\eta_0}+i\frac{\prz H'}{\prz\xi_0}
\label{a8aa}
\ee
where the Hamiltonian $H'$ is given by:
\be
H'=\psi_1+\frac{\gamma}{2}g\;.
\label{a8b}
\ee
For instance, in our case the curve $\cal C$ and the flow region $\cal R$ are the real axis and the positive imaginary half-plane respectively, and the flow is the superimposition of a uniform flow and a vortex. Therefore the vortex Hamiltonian $H'$ (\ref{a8b}) becomes:
\be
H'=q_\infty\eta_0+\frac{\gamma}{4\pi}\log{\eta_0}\;.
\label{a8cb}
\ee

\subsection*{The Routh rule}
When the motion of a vortex is conformally mapped onto the physical $z$-plane, its Hamiltonian has to be corrected according to the so-called Routh rule (Routh 1881). The following is a brief account of its derivation.

The complex conjugate velocity of a vortex in the $z$-plane is
\be
\dot z_0^\star=\dot x_0-i\;\dot y_0=\lim_{z\rightarrow z_0}\left(\frac{dw}{dz}-\frac{\gamma}{2\pi i}\;\frac{1}{z-z_0}\right)\;.
\label{a11}
\ee
This can be written as:
\be
\dot z_0^\star=\lim_{\zeta\rightarrow \zeta_0}\left\lbrack\left(\frac{dw}{d\zeta}-\frac{\gamma}{2\pi i}\frac{1}{\zeta-\zeta_0}\right)   \frac{d\zeta}{dz}+\frac{\gamma}{2\pi i}\left(\frac{1}{(\zeta-\zeta_0)\displaystyle\frac{dz}{d\zeta}}-\frac{1}{z-z_0}\right)\right\rbrack\label{a12}
\ee
that is
\be
\dot z_0^\star=\dot\zeta_0^\star\frac{d\zeta_0}{dz_0}-\frac{\gamma}{2\pi i}\frac{\displaystyle\frac{d^2 z_0}{d\zeta^2_0}}{2\left(\displaystyle\frac{d z_0}{d\zeta_0}\right)^2}
=\dot\zeta_0^\star\frac{d\zeta_0}{dz_0}-\frac{\gamma}{4\pi i}\frac{d}{dz_0}\log\left(\frac{dz_0}{d\zeta_0}\right)\;.
\label{a13}
\ee
Recalling equations (\ref{a8aa}), it follows that
\begin{eqnarray}
\dot z_0^\star&=&\frac{\prz H'}{\prz\xi_0}\frac{\prz\xi_0}{\prz y_0}+\frac{\prz H'}{\prz\eta_0}\frac{\prz\eta_0}{\prz y_0}
+i\left(\frac{\prz H'}{\prz\xi_0}\frac{\prz\xi_0}{\prz x_0}+\frac{\prz H'}{\prz\eta_0}\frac{\prz\eta_0}{\prz x_0}\right)+
\frac{\gamma}{4\pi}\left(\frac{\prz}{\prz y_0}\log\left|\frac{dz_0}{d\zeta_0}\right|+i\frac{\prz}{\prz x_0}\log\left|\frac{dz_0}{d\zeta_0}\right| \right)\nonumber\\
 &=&\left(\frac{\prz}{\prz y_0}+i\frac{\prz}{\prz x_0}\right)\left(H'+\frac{\gamma}{4\pi}\log\left|\frac{dz_0}{d\zeta_0}\right|\right)
\label{a14}
\end{eqnarray}
and finally
\be
H=H'+\frac{\gamma}{4\pi}\log\left|\frac{d z_0}{d \zeta_0}\right|\,.\label{a7b}
\ee
The relationship (\ref{a8b}) between a vortex Hamiltonian and a Green function for a bounded simply connected region was first given by Masotti (1931) (but see also Caldonazzo 1931 and Pelosi 1926). Successively, Lin (1941) extended the Masotti theory to multiply connected domains containing more than one vortex. Surprisingly, Ringleb (1961) seems to have been unaware of the previous studies, as he deduced independently the Routh rule. However, he considered that it was only applicable for mappings based on many-valued functions.

\newpage
\begin{center}
\section*{\normalsize REFERENCES}
\end{center}

\hangindent=4truemm\hangafter=1
\noindent{Aref, H. 1983 Integrable, Chaotic, and Turbulent Vortex
Motion in Two Dimensional Flows. {\em Ann. Rev. Fluid Mech.} {\bf 15}, 345-389.}

\hangindent=4truemm\hangafter=1
\noindent{Bradshaw, P. \& Wong, F. W. F. 1972 The reattachment and relaxation of a turbulent shear flow. {\em J. Fluid Mech.} {\bf 52}, 113-135.}

\hangindent=4truemm\hangafter=1
\noindent{Caldonazzo, B. 1931 Sui moti liquidi piani con un vortice libero. {\em Rend. Circ. Matem. Palermo} {\bf LV}, 369-394.}

\hangindent=4truemm\hangafter=1
\noindent{del Castillo-Negrete, D. 1998 Asymmetric transport and non-Gaussian statistics of passive scalars in vortices in shear. {\em Phys. of Fluids} {\bf 10}, 576-594.}

\hangindent=4truemm\hangafter=1
\noindent{Ceschini, E. 1993 {\em Interferenza vortice-parete}. Laurea thesis, DIASP, Politecnico di Torino.}

\hangindent=4truemm\hangafter=1
\noindent{Clements, R. R. 1973 An inviscid model of two-dimensional vortex shedding. {\em J. Fluid Mech.} {\bf 57}, 321-336.}

\hangindent=4truemm\hangafter=1
\noindent{Cortelezzi, L., Leonard, A. \& Doyle, J. C. 1994 An example of active circulation control of the unsteady separated flow past a semi-infinite plate. {\em J. Fluid Mech.} {\bf 260}, 127-154.}

\hangindent=4truemm\hangafter=1
\noindent{Driver, D. M., Seegmiller, H. L. \& Marvin, J. 1983 Unsteady behaviour of a reattaching shear layer. {\em AIAA Paper} 83/1712.}

\hangindent=4truemm\hangafter=1
\noindent{Driver, D. M., Seegmiller, H. L. \& Marvin, J. 1987 Time-dependent behaviour of a reattaching shear layer. {\em AIAA J.} {\bf 25}, 914-919.}

\hangindent=4truemm\hangafter=1
\noindent{Ferlauto, M. 1996 {\em Studio e controllo di strutture vorticose di parete.} Ph.D. thesis, DIASP, Politecnico di Milano.}

\hangindent=4truemm\hangafter=1
\noindent{Guckenheimer, J. \& Holmes P. 1983 {\em Non-linear oscillations, dynamical systems and bifurcation of vector fields}. Springer-Verlag.}

\hangindent=4truemm\hangafter=1
\noindent{Henrici, P. 1974 {\em Applied and Computational Complex Analysis}, John Wiley \& Sons.}

\hangindent=4truemm\hangafter=1
\noindent{Hosker, R. P. 1984 Flow and diffusion near obstacles. {\em Atmospheric science and power production}, U.S. Dept. of Energy, 241-326.}

\hangindent=4truemm\hangafter=1
\noindent{Hunt, J. C. R., Leibovich, S. \& Richards, K. J. 1988 Turbulent shear flows over low hills. {\em Q. J. R. Meteorol. Soc.} {\bf 114}, 1435-1470.}

\hangindent=4truemm\hangafter=1
\noindent{Le, H., Moin, P. \& Kim, J. 1997 Direct numerical simulation of turbulent flow over a backward-facing step. {\em J. Fluid Mech.} {\bf 330}, 349-374.}

\hangindent=4truemm\hangafter=1
\noindent{Lin, C. C. 1941 On the Motion of Vortices in Two
Dimensions - I. Existence of the Kirchhoff-Routh Function, {\em Proc.
N.A.S.} {\bf 27}, 570-575.}

\hangindent=4truemm\hangafter=1
\noindent{Masotti, A. 1931 Sulla funzione preliminare di Green per un'area piana. {\em Atti Pont. Accad. Sci. Nuovi Lincei} {\bf 84}, 209-216. Sul moto di un vortice rettilineo. {\em Ibidem}, 235-245. Sopra una proprieta' energetica di un vortice rettilineo. {\em Ibidem}, 464-467. Sopra una relazione di reciprocita' nella idrodinamica. {\em Ibidem}, 468-473. Sulle azioni dinamiche dovute ad un vortice rettilineo. {\em Ibidem}, 623-631.}

\hangindent=4truemm\hangafter=1
\noindent{Melnikov, V. K. 1963 On the stability of the center for time periodic perturbations. {\em Trans. Moscow Math. Soc.} {\bf 12}, no. 1.}

\hangindent=4truemm\hangafter=1
\noindent{Nehari, Z. 1975 {\em Conformal Mapping}, Dover Publications Inc.}

\hangindent=4truemm\hangafter=1
\noindent{Novikov, E. A. \& Sedov Y. B. 1979 Stochastization
of Vortices. {\em JETP Lett.} {\bf 29}, no. 12, 677-679.}

\hangindent=4truemm\hangafter=1
\noindent{Pelosi, L. 1926 Un'applicazione idrodinamica della funzione di Green. {\em Atti Reale Accad. Sci. Torino} {\bf LXI}, 569-583.}

\hangindent=4truemm\hangafter=1
\noindent{Pentek, A., Tel, T. \& Toroczkai, Z. 1995 Chaotic advection in the velocity field of leapfrogging vortex pairs. {\em J. Phys. A Math. Gen.} {\bf 28}, 2191-2216.}

\hangindent=4truemm\hangafter=1
\noindent{Ringleb, F. O. 1961 Separation control by trapped vortices, in {\em Boundary Layer and Flow Control}, edited by G. V. Lachmann, Vol. 1, Pergamon Press, 265-294.}

\hangindent=4truemm\hangafter=1
\noindent{Rom-Kedar, V., Leonard, A. \& Wiggins, S. 1990 An analytical study of transport, mixing and chaos in an unsteady vortical flow. {\em J. Fluid Mech.} {\bf 214}, 347-394.}

\hangindent=4truemm\hangafter=1
\noindent{Routh, E. J. 1881 Some applications of conjugate functions. {\em Proc. Lond. Math. Soc.} {\bf 12}, 83.}

\hangindent=4truemm\hangafter=1
\noindent{Tabor, M. 1989 {\em Chaos and integrability in non linear dynamics}, John Wiley \& Sons.}

\hangindent=4truemm\hangafter=1
\noindent{Tsien, H. S. 1943 Symmetrical Joukowsky Airfoils in shear flow. {\em Q. Appl. Math.} {\bf 1}, 130-148.}

\hangindent=4truemm\hangafter=1
\noindent{Wiggins, S. 1992 {\em Chaotic Transport in Dynamical Systems}. Springer-Verlag.}

\hangindent=4truemm\hangafter=1
\noindent{Zannetti, L. \& Franzese, P. 1994 The non-integrability of the restricted problem of two vortices in closed domains. {\em Physica D} {\bf 76}, 99-109.}

\clearpage

\begin{center}
\section*{\normalsize \bf FIGURE CAPTIONS}
\end{center}

Fig. 1: (a) The snow cornice obtained by the Ringleb mapping; (b) the backward-facing step obtained by the Schwartz-Christoffel mapping.

Fig. 2: Simulation of the flow past a snow cornice by a multiple vortex method. The solid circles represent the vortices shed throughout the transient. The vortices strengths are proportional to the symbol sizes.

Fig. 3: (a) Total circulation for the multiple vortex method (solid
line) and circulation for the single vortex model (dotted line), as a function of time; (b) distance of the centroid for the multiple vortex method from the equilibrium point of a single vortex.

Fig. 4: Streamlines and velocity vectors for a separated shear flow for the blunt snow cornice.

Fig. 5: Pattern of the possible trajectories of the vortex. The equilibrium point is marked with a square.

Fig. 6: (a) The separatrix streamline connecting two fixed points ($P^-, P^+$) in the unperturbed case; (b) as a consequence of the perturbation, the streamline splits in two branches which intersect each other: unstable manifold (bold line) and stable manifold (thin line).

Fig. 7: Streamline pattern of a vortex pair flow in the frame of reference moving with the vortices.

Fig. 8: Distribution of scalars in the field behind a snow cornice and concentration time series sampled in the volume marked by the shaded rectangular area downstream, at two different times. The concentration time series is plotted in the upper portion of each frame. The system is periodically perturbed by the vortex autonomously oscillating around its equilibrium location.

Fig. 9: Distribution of scalars in the field behind a backward facing
step with smoothed corner and concentration time series sampled in the volume marked by the shaded rectangular area downstream, at four different times. The concentration time series is plotted in the upper portion of each frame. A random perturbation is imposed on the vortex equilibrium location.

\end{document}