\documentclass{article}
\usepackage{amsmath}

\begin{document}
\subsection{Acoustic Wave Equation}

The equation of acoustic wave motion for a fluid, that in the absence of sound,
has uniform density $\rho$, pressure $P$, and temperature $T$ and is everywhere
at rest will be derived in this section.  We assume that the medium is non-viscous
and has zero heat conductivity.  
$$
\nabla^2p-\frac{1}{c^2} {\partial^2 p \over \partial t^2} = 0
\label{eqn:lineareuler}
$$
where $c$ is the phase velocity and is given by the relation $c= \sqrt{\beta/\rho_0}$, where $\beta$ is the
bulk modulus.

\subsection{Derivation of the Acoustic Wave Equation}

Let\footnote{This follows the derivation in Landau and Lifshitz, \textit{Fluid Mechanics}.}
the pressure $p$ and the density $\rho$ be expressed as small changes
from a constant pressure $p_0$ and density $\rho_0$
$$
p = p_0 + p', \qquad \rho = \rho_0 + \rho'
$$
where $p'$ ane $\rho'$ are the variations and $p'\ll p_0$, $\rho'\ll\rho_0$.
Substituting these expressions into the equation of continuity
$$
{\partial \rho\over\partial t} + \nabla(\rho \mathbf{v}) = 0,
$$

Start a new paragraph just to test indentation.  So that one obtains
$$
{\partial \rho'\over\partial t} + \rho_0\nabla(\mathbf{v}) + \nabla(\rho'\mathbf{v}) = 0,
$$
since $\rho'$, $p'$ and $\mathbf{v}$ are all of first order, the last term can
be neglected to get
$$
{\partial \rho'\over\partial t} + \rho_0\nabla(\mathbf{v}) = 0,
$$
This is the first equation.

Euler's equation is
$$
{\partial\mathbf{v}\over\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} = -{1\over\rho} \nabla p
$$
when one substitutes for the density and pressure
$$
{\partial\mathbf{v}\over\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} = -{1\over\rho_0} \nabla p'
$$
and using the fact that $\rho_0\gg\rho'$.  Keeping only terms of first order
simplifies this equation to
$$
{\partial\mathbf{v}\over\partial t} + {1\over\rho_0} \nabla p'=0
$$
This is the second equation.

The third equation comes from the fact that a sound wave in an ideal
fluid is adiabatic.  This means that the propagation does not change
the entropy of the system.  Thus
$$
p' = \left({\partial p \over \partial \rho_0}\right)_s  \rho'
$$
where the $s$ denotes constant entropy.  Substituting into the continuity equation yields
$$
{\partial p'\over\partial t} + \rho_0\left({\partial p\over \partial \rho_0}\right)_s\nabla\mathbf{v} = 0,
$$

To combine the two equations into one, then introduce a velocity potential $\psi$, which
is related to the pressure by
$$
p' = -\rho {\partial \psi\over\partial t}
$$
The velocity potential is related to the velocity by
$$
\mathbf{v}(\mathbf{r},t)=\nabla \psi(\mathbf{r},t)
$$
so we obtain
$$
-\rho {\partial^2 \psi\over\partial^2 t} + \rho_0\left({\partial p\over \partial \rho_0}\right)_s\nabla^2 \psi = 0,
$$
Now if we define the thermal compressiblity
$$
\beta = \left({\partial p\over \partial \rho_0}\right)_s
$$
then
$$
{\partial^2 \psi\over\partial^2 t} = {\rho_0\over \beta} \nabla^2 \psi = 0,
$$
The speed of sound is $c=\sqrt{\beta/\rho_0}$.

\subsection{Velocity Potential}

Since the curl of a gradient vanishes, equation~\ref{eqn:lineareuler}
shows that the particle velocity is irrotational, so that it can be
expressed as the gradient of a scalar function, $\mathbf{v}=\nabla \phi$,
where $\phi$ is the velocity potential.  Substituting the velocity
potential into the Euler equation, ~\ref{eqn:lineareuler}, shows that
the time derivative of the velocity potential is related to the acoustic
pressure by the density
$$
p=-\rho_0 {\partial \phi\over \partial t}
$$
Thus, the velocity potential is also a solution to the wave equation. 
The velocity potential can be interpreted physically as the absence or
shear, turbulence, or rotational flow, since the particle velocity is
the gradient of the velocity potential, and thus its curl must vanish. 
This is true only with the idealized derivation shown here, though this
linear model does extremely well for acoustic propagation for non--shock
waves.  The velocity potential is useful in photoacoustics
in that it can be easily obtained by integrating the pressure over time,
giving a positive waveform that can be used, for example, in
backprojection image reconstruction.\footnote{Viator's thesis}

\subsection{Stress Confinement}
The acoustic characteristic time is
$$
t_c = d/c
$$
where $d$ is the characteristic distance in the problem and $c$ is the speed of sound.
The characteristic distance is typically the smallest dimension of interest
in the problem.  For example, a sphere with radius $r$ in an unbounded medium
will have $d=2r$.  However, if that same sphere is irradiated with light and the
absorption coefficient for the sphere is sufficiently high $1/\mu_a\gg r$ then
the light will be deposited as thin cap on the sphere.  In this case the characteristic
dimension is the thickness of that cap and so $d = 1/\mu_a$.  Of course, if the 
absorption coefficient is lower $r < 1/\mu_a$ then the characteristic dimension 
becomes the radius once again.  

It is also convenient to define a dimensionless laser pulse duration in terms of
the characteristic time as
$$
\tau_p = t_p/t_c
$$


\subsection{General Solution of the Acoustic Wave Equation}
Since in soft tissues the optical contrast, which determines the optoacoustic sources, is generally 
much higher than the small variations of acoustic impedance, 
an acoustically homogeneous medium is assumed. 

The propagation of acoustic waves is based on the linear wave equation
$$
\nabla^2 P(\mathbf{r},t) - {1\over c^2} {\partial^2 P(\mathbf{r},t) \over \partial t^2} = 0
$$
where $P$ denotes the pressure and $c$ the speed of sound in the medium.
This equation must be supplied with apropriate boundary conditions,
which we will neglect for the moment. For pulsed optoacoustic wave
generation and for finite element formulations it
is, however, more appropriate to work with the acoustic velocity
potential $\psi$ and not with directly with the acoustic pressure $P$. These quantities and
the particle velocity $\mathbf{v}$ in the acoustic medium are related by
$$
P = \rho {\partial \psi(\mathbf{r},t)\over\partial t} \qquad \mathrm{and}\qquad \mathbf{v}(\mathbf{r},t)=\nabla \psi(\mathbf{r},t)
$$
where $\rho$ is the density.  It is easily observed, that in the linear case the acoustic velocity
potential satisfies the same equation as the acoustic pressure
$$
\nabla^2\psi(\mathbf{r},t) - {1\over c^2} {\partial^2 \psi(\mathbf{r},t) \over \partial t^2} = 0
$$
When there is a acoustic source then the equation becomes
$$
\nabla^2\psi(\mathbf{r},t) - {1\over c^2} {\partial^2 \psi(\mathbf{r},t) \over \partial t^2} = {\beta\over \rho c_p} S(\mathbf{r},t)
$$
where $c_p$ is the specific heat
capacity at constant pressure, $\beta$ is the thermal expansion
coefficient, and $S(\mathbf{r},t)$ is the heat generated per unit volume
and time.  

A further assumption 
is that heating is instantaneous and is described by a delta 
function at time $t = 0$. It is advantageous to use 
the solution for the velocity potential, which is related to 
the acoustic pressure by
$$
P(\mathbf{r},t)= - \rho {\partial \psi(\mathbf{r},t)\over \partial t}
$$
where $\psi$ is the velocity potential, $\rho$ the density 
and $P$ the acoustic pressure. A solution based on a retarded Green's 
function is given by the integral.
$$
\psi (\mathbf{r}_i ,t)=-{\beta \over 4\pi\rho c_p } {1\over t}
\iint\limits_{\left| \mathbf{r}'-\mathbf{r}_i \right| =c t} W(r')\,dS
$$
where $W$ the volumetric 
energy density caused by absorption of the laser pulse. The positions
$\mathbf{r}'$ and $\mathbf{r}_i$ 
are points in the source volume and the position of the i$^{th}$ 
detector in the array, respectively. The surface element $dS$ allows 
integration over the surface of a sphere with radius $c t$ from 
the detector point. A solution for a laser pulse with finite 
duration is obtained by convolving $\psi$ with the temporal 
pulse profile $g(t)$. 
$$
\psi (t)=\psi_\delta (t) \otimes g(t)
$$
The symbol $\otimes$ denotes convolution and $\psi_\delta$ 
is the solution of Eq. (2). A Gaussian function 
with a $1/e$ width corresponding to the pulse duration is a
convenient representation for $g(t)$
$$
g(t)=e^{-(2t/t_p )^2}
$$
To calculate the velocity potential 
for an arbitrary energy distribution $W$, a finite sum over small 
sources distributed over a grid with spacings $\Delta x$, $\Delta y$,
and $\Delta z$ replaced the integral in Eq.(2). To avoid discretization 
noise, these sources had to overlap. A method to create overlapping 
sources is to first assume point sources and to perform the summation. 
The noisy result of this summation is convolved with a Gaussian 
function (using Eq.3) that has a $1/e$ duration equal to the average 
propagation time of the acoustic wave through a volume element
$$
t_{prop} ={\Delta x+\Delta y+\Delta z \over 3 c} 
$$
Ideally, this time should be identical to the laser pulse duration. 
Due to computer storage limitations we had to choose a larger 
voxel size and $t_{prop}$ replaced $t_p$ in Eq. (3). For the sake of 
simplicity, the operation described in (3) that transforms information 
about the absorbed energy distribution W into velocity potential 
signals $\psi$ is symbolically written as
$$
\psi =A W
$$
where $A$ denotes the optoacoustic ``transform''. For the back projection 
first a quantity $s_i(t)$ is derived from Eq.(2) that is given 
by the surface integral of $W$ over a sphere around the i$^{th}$ detector,
$$
s_i (t)\equiv \iint\limits_{\left| r' -r_i \right| =c t} W(r' ) dS=
-{4\pi  \rho  C_p \over \beta }   t  \psi (r_i ,t)
$$
$s_i$ has some similarities with the 3-D Radon transform with 
its value at a given time being determined by a surface integral 
over $W$. Where the Radon transform uses planar projections in 
different directions, the operation described in Eq.(7) uses 
integration over the surfaces of spheres and the direction of 
the projection is given by the detector position, $r_i$ 

\subsection{Acoustic Source}
\subsection{Point Source Solution}
The Green's function for the Wave Equation

Let\footnote{Morse and Feshbach, \textit{Methods of Theoretical Physics}}
$G$ be a solution of
$$
\nabla^2G(\mathbf{r},t; \mathbf{r}_0,t_0)-{1\over c^2}{\partial^2\over\partial t^2}G(\mathbf{r},t; \mathbf{r}_0,t_0)= -4\pi\delta(\mathbf{r}-\mathbf{r}_0)\delta(t-t_0)
$$
satisfying homogeneous boundary conditions on surface S and obeing the causality requirement that $G$ and 
${\partial G\over\partial t}=0$ everywhere for $t<t_0$.  The reciprocity relation is
then
$$
G(\mathbf{r},t\vert\mathbf{r}_0,t_0) = G(\mathbf{r},-t\vert\mathbf{r}_0,-t_0)
$$
The closed form solution fot he Green's function fo the infinite domian is
$$
G(\mathbf{r},t\vert\mathbf{r}_0,t_0) = {1\over \vert\mathbf{r}-\mathbf{r}_0\vert}\delta(\vert\mathbf{r}-\mathbf{r}_0\vert/c - (t-t_0))
$$

\subsection{Solid Sphere Solution}
\subsection{Disk Solution}
\subsection{Boundary Conditions}

\section{Preliminaries}
\subsection{Euler's Equation}
For an inviscid fluid, one that has no viscosity, a fluid particle
accelerates according to Newton's second law,
$\mathbf{a}=d\mathbf{f}/dm$, where $\mathbf{a}$ is the acceleration,
$d\mathbf{f}$ is the force, and $dm$ is the mass of the fluid particle. 
The force can be written
$$
d\mathbf{f}=-\nabla P dV
$$
The acceleration is
$$
\mathbf{a}={\partial \mathbf{v}\over \partial t} + u_x{\partial \mathbf{v}\over \partial x} + u_y\frac{\partial
\mathbf{v}}{\partial y} + u_z{\partial \mathbf{v}\over \partial z}
$$
With the operator defined as
$$
(\mathbf{v} \cdot \nabla) = u_x{\partial\over \partial x} + u_y{\partial\over \partial y} +
u_z{\partial\over \partial z}
$$
then the acceleration is
$$
\mathbf{a}={\partial \mathbf{v}\over \partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v}
$$
With the substitution of $dm=\rho dV$ into $d\mathbf{f}=\mathbf{a} dm$ we have the Euler equation
$$
-\nabla P = \rho \left [ {\partial \mathbf{v}\over \partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} \right ]
$$

\subsection{Equation of Continuity}
The equation of continuity relates the net flow through
a volume element with the rate of mass increase and is expressed
$$
{\partial \mathbf{v}\over\partial t} + u_x{\partial \mathbf{v}\over \partial t} + \nabla \cdot(\rho \mathbf{v}) =0
\label{eqn:continuity}
$$
where $\mathbf{v}$ is the particle velocity.  Using the fact that $\rho
= \rho_0 (1+s)$ and $\rho_0$ is constant for both spatial and time
derivatives, then the  equation of continuity~\ref{eqn:continuity}
becomes linear
$$
{\partial s\over \partial t} + \nabla \cdot \mathbf{v} =0
\label{eqn:linearcontinuity}
$$

\subsection{Equation of State}
The adiabatic equation of state is 
$$  
p={P\over P_0} =\rho_0 \left ({\partial P \over\partial \rho} \right)_{\rho_0, s}
\label{eqn:eosadiabatic}
$$
where $p$ is the acoustic pressure and $s={\rho-\rho_0\over \rho_0}$ is the condensation.  


\end{document}