\def\w{\omega}
\def\xiv{\mbox{\boldmath$\xi$}}
\INT
Passive seismic imaging is an example of wavefield interferometric
imaging.  In this case, the goal is the production of subsurface
structural images by recording the ambient noisefield of the Earth
with surface arrays of seismometers or geophones.  The images produced
with this technique are directly analogous to those produced with
conventional reflection seismic data.  Interferometric imaging is not
new science.  This statement is true for two reasons.  
\par
First, the idea of imaging the subsurface by correlating the 
recordings from subsurface sources was first  
introduced by \longcite{GEO33-02-02640269}.  That work 
provided a one-dimensional proof that the auto-correlation of time
series collected on the surface of the Earth can produce the
equivalent to a zero-offset time section.  Subsequently, 
\longcite{Zhang.sep.60.85.art}, used planewave decomposition 
to prove the result in 3D over a homogeneous medium.  \longcite{fink}
developed the Green's function of a heterogeneous medium
with acoustic waves via correlation and validated the theory with an
ultrasonic experiment.  \longcite{kees04}, through one-way
reciprocity, 
proved that cross-correlating traces of the transmission response of
an arbitrary medium synthesizes the complete reflection response,
in the form of shot-gathers, collected in a conventional active source
experiment.  In fact, in recent years, there has been a sudden
increase in interest around interferometry that has generated a wealth
of literature \cite{speced}.  
\section{Theory}
\subsection{Natural phase encoding}
The most intuitive explanation of how correlation of subsurface
sources returns the conventional active seismic data volume is in
terms of a phase encoding experiment.  Phase encoding shot records for
simultaneous processing was discussed in
\longcite{GEO65-02-04260436}.  The the motivation for phase encoding
is to manipulate the independant shot gathers in a predictable manner
such that the entire data volume can be processed (in this case
migrated) simultaneously.  To that end, each shot gather is convolved
with a unique function that will minimize cross-talk
interference between the shots during a summation over the shot axes
in a 3D seismic survey.  Random time series are quite effective
encoding functions.  Planewave migration strategies can also be viewed
as a phase encoding migration operation.  Consider a set of encoding
functions $a({\bf x};{\bf x}_s,\w)$ for each shot ${\bf x}_s$ with
which shot records are convolved.  Normally, source functions for
active seismic data are bandlimited impulses.  After convolving the
impulses with encoding functions the source wavefields
$\hat{S}$ become  
\begin{equation}
\hat{S}_{z=0}({\bf x},\w)=\sum_{{\bf x}_s} 
a({\bf x};{\bf x}_s,\w)I_{z=0}({\bf x}_s,\w)\;.
\end{equation}
The five dimensional active seismic data, $R$, is similarly convolved with
encoding functions and summed to reduce its dimensionality
\begin{equation}\label{encode}
\hat{R}_{z=0}({\bf x},\w)=\sum_{{\bf x}_s} 
a({\bf x};{\bf x}_s,\w)R_{z=0}({\bf x};{\bf x}_s,\w)\;.
\end{equation}
\par
Special cases for the encoding function have important meaning.  If
$a({\bf x};{\bf x}_s,\w)\equiv I_{z=0}({\bf x}={\bf x}_s,\w)$, the
resulting volume is that obtained by resorting to midpoint-offset
coordinates.  For seismic interferometry, we assume unknown and
arbitrarily complicated encoding functions arrive from the
subsurface.  Therefore a transmission record $T({\bf x},\w)$ recorded
by an array at the surface is exactly
$\hat{R}_{z=0}$.  To recover $R_{z=0}({\bf x};{\bf x}_s,\w)$, we need
to invert equation~\ref{encode}.  Since the adjoint of convolution is
correlation,  correlating every trace with every other trace would
return a scaled version of the active data $R$ except for the presenc
of the sum in equation~\ref{encode}.  Unfortunately, the correlation
of a sum is not equal to the sum of correlations.  Therefore, the
summation over ${\bf x}_s$ will introduce cross-talk artifacts that
are convolutions of the encoding functions and the shot records one is
trying to recover.  The simple way to avoid this problem is to
parameterize the transmission wavefields as a function of individual
subsurface (encoding) functions, $\xiv$, and move the summation
outside of the correlation.  Thus, 
\begin{equation}\label{uncode}
\tilde{R}_{z=0}({\bf x};{\bf x}_s,\w)=\sum_\xiv T({\bf x};\xiv,\w)T^*({\bf
  x}_s,\xiv,\w)=\sum_\xiv|a({\bf x},\xiv,\w)|^2R_{z=0}({\bf x};{\bf
  x}_s,\w)\;. 
\end{equation}
Above, $\tilde{R}$ includes a scalar amplitude term that
can be removed by normalizing by the amplitude spectrum, though is
kinematically equivalent to the active source data assuming that the
sum over $\xiv$ is dense.  If too few subsurface sources are captured,
the resulting volume will not reproduce the kinematics of the active
survey.  In terms of the phase encoding analogy above, this would be
equivalent to decimating the shot axis of the active survery which
will clearly lead to illumination problems.
\par
Figure~\ref{fig:sketch}, simplified for clarity,
shows the basic kinematics exploited in processing transmission
data.  The figure includes two recording stations 
capturing an approximately planar wavefront emerging from a two-layer
subsurface.  Panel ({\tt a}) shows the raypaths associated with the direct
arrival and one reflected both at the free surface and a subsurface
interface.  The second travel path (labeled reflection ray) has the 
familiar kinematics of the reflection seismic experiment if a source
were excited at receiver one.  The transmission
wavefield is depicted in panel ({\tt b}).  Wavelet polarity is
appropriate for direct arrivals and reflection.  The three main
features of transmission wavefields can be appreciated here.  First, the
exact timing of the arrival is unknown.  Second, the phase 
and duration of transmitted energy are unknown and likely
complicated.  Third, if the incident wavetrain coda is long,
arrivals in the transmission record can interfere.
\par
%%%%\sideplot{figure1}{width=3.3in, height=2.4in}{ ({\tt a}) Approximately
\sideplot{sketch}{width=3.3in,height=2.4in}{({\tt a}) Approximately
  planar arrival with rays showing important propagation paths for
  passive imaging.  ({\tt b}) Idealized traces from a transmission
  wavefield. ({\tt c}) Shot-gather (reflection wavefield) synthesized
  using trace \textit{r1} as the source.  Many details, explained in
  the text, are neglected for simplicity.} 
Choosing trace \textit{r1} as the comparison trace, panel ({\tt c})
depicts the correlation spikes associated with the arrivals in the
data panel ({\tt b}), where $\otimes$ is correlation.  A solid line
with linear move-out is super-imposed across the correlated traces
corresponding to the direct arrival recorded at each receiver
location.  The dashed line on panel ({\tt c}) has hyperbolic moveout.
However, no correlation peak exists on the \textit{r1}$\otimes$
\textit{r1} trace under the hyperbola.  Not drawn, the second arrival
on \textit{r2} will have a counterpart on \textit{r1} from a ray
reaching the free surface farther to the left of the model.  In fact,
the correlations produced from a single planewave will produce another
planewave. 
\par
However, each planar reflection is moved to the lag-time associated with
a two-way trip from the surface to the reflector.  Correlation 
removes the wait time for the initial arrival and maintains the 
time differences between the direct arrival and reflections.  Summing the 
correlations from a full suite of planewaves builds hyperbolic events 
through constructive and destructive interference.  Analyzing seismic 
data in terms of planewave constituents is a commonly invoked tool in 
seismic processing.  Summing the correlations from incident planewaves 
is a planewave superposition process that builds hyperbolas.
\par
Passive seismic imaging is predicated on raypaths bouncing every
which-way from every direction.  Cartoons depicting the transmission 
experiment always leave something out that
causes an inconsistency that needs more raypaths and receivers to
explain.  Unfortunately the trend continues nearly forever.  Such a
complication arises with the inclusion of a second 
reflector.  The two reflection rays correlate with each other with a
positive coefficient.  The two travel paths share the time through the
shallow layer, so they correlate at a lag (or time difference), equal
to the two-way travel time through the deep layer.  This correlation
is not a problem however.  
\par
An upcoming interbed multiple event arrives after the direct arrival
delayed by the two-way travel-time through the deeper layer.  This event
has the opposite polarity compared to the direct arrival after once
changing its propagation direction from $\uparrow$ to $\downarrow$.  The 
delay of its arrival at receiver \textit{r2} compared to the direct
arrival at receiver \textit{r1} is the same lag as the correlation
between the reflectors and has opposite sign.  Therefore, the internal
multiple cancels potential artifacts of the correlation.  This shows
the importance of modeling transmission data with a two-way 
extrapolator.  Without all possible multiples, correlation artifacts
will quickly overwhelm the Earth structure within the correlated output. 
\par
Cross-correlation of each trace with every other trace handles the
three main difficulties of transmission data: timing, waveform, and
interference.  First, the output of the correlation is in lag units,
that when multiplied by the time sampling interval, provide the
time delays between like events on different traces.  The zero lag of the
correlation is the zero time for the synthesized 
shot-gathers.  Second, each trace records the character and duration
of the incident energy as it is reflected at the surface.  This
becomes the source wavelet analogous to a recorded vibrator 
sweep.  Third, overlapping wavelets are separated by correlation. 
\par
To calculate the Fourier transform of the reflection response of a
loss-less subsurface, $R({\bf x}_r,{\bf x}_s, \omega)$,
\longcite{kees04}, proved
\begin{equation}\label{rigor}
2\Re [R({\bf x}_r,{\bf x}_s,\omega)]=
\delta({\bf x}_s-{\bf x}_r)-\int_{\partial\! D_m} 
T  ({\bf x}_r,\xiv,\omega)
T^*({\bf x}_s,\xiv,\omega)\;
\mbox{d}^2\xiv \;,
\end{equation}
where $^*$ represents conjugation.  The vector ${\bf x}$ will
correspond herein to horizontal coordinates, where subscripts $r$ and
$s$ indicate different station  locations from a transmission
wavefield.  After correlation $r$ and $s$ acquire the meaning of receiver
and source locations, respectively, associated with an active survey.  The 
RHS represents summing correlations of windows of passive data around
the  occurrence of individual sources from locations $\xiv$.  The
symbol $\partial \! D_m$ represents the domain boundary that surrounds
the subsurface region of interest on which the sources are located.  The 
transmission wavefields $T({\bf x}_r,\xiv,\omega)$ contain the arrival and
reverberations due to only one subsurface source.  To synthesize the
reflection experiment exactly, impulsive sources should completely
surround the volume of the subsurface one is trying to image.  Alternatively, 
many impulses can be substituted with a full suite of planewaves
emerging from all angles and azimuths as in the kinematic explanation
above.
\par
Figure~\ref{fig:domain} is a cartoon that follows the derivation of
the formula given above \cite{kees05}.  Beginning with Gauss' theorem,
it is known that the net flux through a domain boundary must be equal to
the net energy production within the volume of the domain.  
\begin{equation}
\int_D pp
\end{equation}
The physical
experiment representing this is depicted in the upper sketch of the
figure.  Two sources are pictured within the domain, and we
acknowledge that receivers are needed everywhere around the boundary
to make measurements.  Two arbitrary raypaths are drawn out of the
nearly infinite set.  The second sketch switches the roles of sources
and receivers via source-receiver reciprocity.  Of the many
derivations of seismic interferometry produced in the literature, the
lack of constraints on the geometry of the measurements makes this
result very general.  The final sketch shows the special case when the
domain is truncated by a free surface.  The raypath recorded at the
left receiver from the top source is now recovered by a reflection
from the free surface due to the lowest source.  An unknown time delay
has been introduced however.  
\par
For any receiver location in the above cartoon, the physical
experiment shown would simply quantify the size of the energy sink
needed to remove all the acoustic energy within the domain due to the
sources surrounding it.  However, the energy conservation relations
can be rewritten in terms of Green's functions instead of particle
velocities.  Then the volume integral becomes the
difference between the Green's functions at two locations, and 







