\section{Introduction}
\label{intro} Simulations of penetrative convection above a
stellar-type convection zone have been performed in two-dimensions
by Hurlburt et al. (1986), and in three-dimensions by Singh et al.
(1994, 2001) and Robinson et al. (2004). In the studies by Hurlburt
et al. and Singh et al., several models were computed to study the
overshoot of convective motions from unstable into the upper stable
layer in a non-rotating configuration. It was found that the motions
from the convective region penetrate a significant fraction of
pressure scale height into the stable layer above. Singh et al.
(1994) also found that the penetration distance above the convection
zone ($\Delta_u$) scales as $\Delta_u \sim (F_b/\rho_{ctop})^{1/3}$,
where $F_b$ is the input flux and $\rho_{ctop}$ is the density at
the top of the convection zone. Robinson et al. (2004) performed
three-dimensional simulations of the upper radiation-convection
transition layer with more realistic modeling for three subgiant
stars. For the present Sun they found the overshoot to be $0.5 H_p$,
for the $11.3$-Gyr subgiant it was $0.75 H_p$, while for a
$11.6$-Gyr subgiant it was approximately $1 H_p$, $H_p$ being the
local pressure scale height.

Simulations of the effects of rotation on convective penetration or
overshooting have been studied by Brummel et al. (2002), Ziegler \&
R\"{u}diger (2003), Browning et al. (2004), K\"{a}pyl\"{a} et al.
(2004), and Pal et al. (2007).
 Brummel et al. (2002) and Pal et al.
(2007) examined the behaviour of penetrative convection below a
convection zone under the influence of rotation by means of
three-dimensional simulations. They found that with an increase in
rotational velocity, the downward penetration decreased. A similar
behaviour was observed when the stability of the lower stable layer
was increased in a rotating configuration. Furthermore, the relative
stability parameter $S$ showed an $S^{-1/4}$ dependence on the
penetration distance implying the existence of a thermal adjustment
region in the lower stable layer rather than a nearly adiabatic
penetration zone. Ziegler \& R\"{u}diger (2003) and K\"{a}pyl\"{a}
et al. (2004), using their 3D MHD codes also found that, as a
general feature, the overshooting at the bottom decreases as a
function of increasing rotation at a given latitude.

Browning et al. (2004) performed three-dimensional simulations of
core convection within A-type stars at a range of rotation rates.
They found that as convective motions enter the stable
stratification of the radiative envelope, increasing the rotation
rate enhances the overshooting.

In this paper, we perform large eddy simulations of turbulent
convection in a {\it{f}}-plane configuration to study the effect of
rotation on the penetration distance above the convection zone. A
total of 17  two-layer (lower unstable -- upper stable) models have
been set up to examine the effect of rotation rate
({\boldmath$\rm\Omega$}) and colatitude ($\theta$) on the
penetration distance. In some models, we systematically vary the
stability of the upper stable layer to study its effect on the
penetration distance. In one model, we increase the horizontal
resolution to see if the penetration distance is altered.

In the next section, we give the essential ingredients of the
simulations and the parameters of the computed models. In Section 3,
we provide the results and their discussion. Important conclusions
of the study are listed in Section 4.

\section{Parameters of Computed Models}

The general behaviour of convective transport in a stellar-type
setting has been studied using Large Eddy Simulation (LES) approach
by several groups (Chan \& Sofia 1986, Hossain \& Mullan 1991, 1993,
Muthsam et al. 1995, Singh \& Chan 1993, Singh et al. 1994, 1995,
1996, 1998a,b, 2001 Saikia et al. 2000, Chan 2001, Pal et al. 2007).
LES, being less demanding on speed and memory as it can have a
coarser grid, allows large scale flows to be  modeled explicitly
while the smaller scales are modeled by some sort of sub-grid scale
formulation (Smagorinsky 1963).

Pal et al. (2007) numerically solved the Navi\'{e}r-Stokes equations
and incorporated the rotational effects by considering an
{\it{f}}-plane configuration. A plane parallel layer of perfect gas
was considered in a rectangular box which can be viewed as a small
portion of a spherical shell. We use a configuration that is similar
to that of Pal et al. (2007) except that the computational domain
has two layers (lower unstable -- upper stable) rather than their
three layer sandwich (stable-unstable-stable) configuration. The
choice is justified as we are interested in studying the behaviour
of penetrative convection above a convection zone rather than below
as was the aim of Pal et al. (2007). The spherical shell rotates
around the polar axis from west to east  and the angular velocity
vector {\boldmath$\rm\Omega$} points toward the north pole. We use a
right handed cartesian coordinate system in which $X$ and $Y$ denote
the horizontal direction and $Z$ denotes the upward vertical
direction. The gravity vector is denoted by {\boldmath$\rm g$} and
is along the negative $Z$ direction.  The angle between
{\boldmath$\rm\Omega$} and the $Z$ axis is denoted by $\theta$ and
the rotation vector {\boldmath$\rm\Omega$} lies in the $XZ$ plane.
The tilted rotation vector is kept constant implying a uniform
angular velocity of the sphere.

The rectangular computational domain has an aspect ratio of 1.5 and
a mesh of $35 \times 35 \times 96$ points. The top and bottom
boundaries are kept impenetrable and stress-free while the side
boundaries are made periodic. The domain is divided into two layers
with the convectively unstable layer positioned below the
convectively stable one. A constant flux {\boldmath$\rm{F}_b$} is
fed from the bottom and all the thermodynamic variables are
expressed in units which set the total depth and initial density,
pressure and temperature at the top to unity.

We use the Large Eddy Simulations (LES) approach to solve the
Navi\'{e}r-Stokes equation for an ideal gas having the ratio of
specific heats ($\gamma$) as $5/3$. The equations of the problem are
given in Pal et al. (2007). We use the Smagorinsky coefficient of
viscosity to represent the sub-grid-scale eddy viscosity in the
form:

\begin{equation}
\mu=\rho(c_{\mu}\Delta)^2(2\mathbf\sigma:\sigma)^{1/2},
\end{equation}
where $c_{\mu}$ is the Deardorff coefficient (Deardorff 1971),
$\sigma$ is the strain rate tensor and the colon sign inside the
bracket denotes tensor contraction and
$\rm\Delta^2=(\rm\Delta_x\rm\Delta_y)^{1/2}\rm\Delta_z$. The SGS
turbulent diffusivity is computed from this viscosity by assuming a
constant Prandtl number $\rm{Pr}=1/3$

\begin{table}
\caption{Physical Parameters for Models R1 to R17}\label{table:1}
\centering
\begin{tabular}{cccccrrrcc}
\hline\hline &     &       & \multicolumn{2}{c}{Polytropic Indices}&
&
        & &\multicolumn{2}{c}{Layer Thickness(PSH)} \\
\cline{4-5} \cline{9-10}\\
Models & Grid & $\rm{F_b}$ & Top & Bottom & $P_b$ & $T_b$ &$\rho_b$ & Top & Bottom \\
\hline
R1-R13 & $35 \times 35 \times 96 $ & 0.125 & 2.0 & 1.5 & 4655 & 19.5 & 238 & 6.04 & 2.40 \\
R14 & $35 \times 35 \times 96$ & 0.125 & 3.0 & 1.5 & 19767 & 17.8 & 1106 & 7.12 & 2.77 \\
R15 & $35 \times 35 \times 96$ & 0.125 & 4.0 & 1.5 & 65092 & 16.8 & 3862 & 8.01 & 3.07 \\
R16 & $35 \times 35 \times 96$ & 0.125 & 5.0 & 1.5 & 178290 & 16.1 & 11024 & 8.76 & 3.33 \\
R17 & $46 \times 46 \times 96$ & 0.125 & 2.0 & 1.5 & 4655 & 19.5 & 238 & 6.04 & 2.40 \\
\hline
\end{tabular}
\end{table}

The governing fluid equations are solved using an implicit scheme
(Chan \& Wolf 1982) and the time step is characterized by the CFL
number

\begin{equation}
N_{CFL}=C_s\Delta{t}/\Delta_{min},
\end{equation}
where $\rm\Delta_{min}$ is the minimum grid size in any direction
and $C_s={\gamma}^{1/2}$ is the dimensionless sound speed at the
top. The boundary conditions imposed are :
\begin{eqnarray}
{\mathbf {F_b}} = constant = 0.125 \,\,\,(at \,\,the\,\,
bottom),\nonumber
\end{eqnarray}
\begin{eqnarray}
T=T_t=constant \,\,\,(at \,\,the \,\,top),\nonumber
\end{eqnarray}
\begin{eqnarray}
v_z=0\,\,;\,\,\,\frac{\partial{v_x}}{\partial{z}}=\frac{\partial{v_y}}{\partial{z}}=0.\nonumber
\end{eqnarray}

\begin{figure}
\centering{\includegraphics[width=8cm]{fig1new.eps}} \caption
{Profiles of the time- and horizontally-averaged energy fluxes for
(a) case R1 $(\mathbf{\Omega}=0,\, \theta=0\,^\circ)$, (b) case R7
($\mathbf{\Omega}=0.5,\, \theta=45\,^\circ$). \label{fig.1}}
\end{figure}

We have computed $17$ models and the physical parameters for the
models are listed in Table 1. The numerical parameters of the
simulations are given in Table 2. Columns (6) - (9) in Table 3 list
the values of Coriolis number ($Co$), Rossby number ($Ro$),
effective Reynolds number ($Re$), and Taylor number ($Ta$) for the
simulated cases. Following Chan (2001), they are defined as
\begin{eqnarray}
Co={\mathbf\Omega} d /\langle v^{\prime\prime} \rangle, \nonumber
\end{eqnarray}
\begin{eqnarray}
Ro=1/Co,                     \nonumber
\end{eqnarray}
\begin{eqnarray}
Re=\langle v^{\prime\prime} \rangle d/ \langle
\bar{\mu}/{\bar{\rho}} \rangle,      \nonumber
\end{eqnarray}
\begin{eqnarray}
Ta=(2{\mathbf\Omega}d^2/ \langle \bar{\mu}/\bar{\rho} \rangle)^2,
\nonumber
\end{eqnarray}
with $v^{\prime\prime} \equiv ({{v_x}^{\prime\prime}}^2+
{{v_y}^{\prime\prime}}^2+{{v_z}^{\prime\prime}}^2)^{1/2}$ where
$v_x^{\prime\prime}$ denotes the root mean square (rms) fluctuation
of $v_x$ etc. The extent of domain is denoted by $d$.

Models R1 to R13 are computed to examine the effect of rate of
rotation ({\boldmath$\rm\Omega$}) and the angle ($\theta$) between
the rotation vector and the vertical axis on the extent of
penetration below the convectively stable layer. The polytropic
index of the upper stable layer is taken to be 2.0 in all these
models (R1-R13). In models R2 to R5, $\theta$ is varied from
$22.5\,^\circ$ to $90\,^\circ$ while {\boldmath$\rm\Omega$} is fixed
at 0.25. In the second set of models R0, R10, R11 and R13,
{\boldmath$\rm\Omega$} is varied from 0 to 1 in steps of 0.25,
respectively, keeping $\theta$ fixed at $0\,^\circ$. In models R3,
R7, and R12, {\boldmath$\rm\Omega$} is 0.25, 0.5 and 1.0,
respectively, while $\theta$ is fixed at $45\,^\circ$. In the last
set of models R14 to R16, polytropic index of the upper stable layer
is varied from 3 to 5 while all the other parameters are the same as
model R7. Another model R17, having a grid of $46 \times 46 \times
96$ points and all other parameters same as R7
({\boldmath$\rm\Omega$} $= 0.5,\, \theta = 45\,^\circ$) is computed
to examine the effect of increase in the number of grid points
(resolution) in the horizontal plane.

As described , four models, namely, R7, R14, R15 and R16 have been
set-up to examine the effect of stability of the upper stable layer
on the penetration height of these rotating configurations. Each of
these four models have different polytropic indices in the upper
stable layer, namely, 2, 3, 4 and 5 for the cases R7, R14, R15 and
R16, respectively. For all these four models, the values of
{\boldmath$\rm\Omega$} and $\theta$ are kept constant at 0.5 and
$45\,^\circ$ respectively. Similar to the relative stability
parameter $S$ for the lower stable layer (Hurlburt et al. 1994,
Singh et al. 1995, Pal et al. 2007), we can define a relative
stability parameter for the upper stable layer for these four cases
as:

\begin{equation}
S=\frac{m_{i}-m_a}{m_{2}-m_a}\,\,;\,\, m_i=2,3,4,5,
\end{equation}
where $m_{i}$ denote the polytropic indices of the upper stable
layer for various cases, $m_{2}$ is the polytropic index of the
upper stable layer for our case R7 and $m_a=1/(\gamma-1)$ is the
adiabatic index. Here, we have taken $\gamma=c_p/c_v$ as 5/3. Thus,
for our reference case R7 the relative stability parameter $S$ is
equal to unity and for cases R14, R15 and R16 it comes out to be 3,
5 and 7, respectively.

\begin{table}
\caption{Numerical Parameters for the 17 models}
\label{table:2}
\centering
\begin{tabular}{ccccccc}
\hline\hline
Models & $g$ & $C_{\mu}$ & $\rm\Delta t$ & $t$ & $N_{CFL}$ & $\rm\Delta_{min}$ \\
\hline
R1-R15,R17 &  50 & 0.2 & 0.00054795 & 1698 &  0.1 & 0.007074 \\
R16        & 50  & 0.2 & 0.00027397 & 1698 & 0.05 & 0.007074 \\
\hline
\end{tabular}
\end{table}

\begin{table}
\caption{Dynamical Parameters of the Computed Models}
\label{table:3}
\centering
\begin{tabular}{ccccccccccc}
\hline\hline
Model&{\boldmath$\rm\Omega$}  & $\theta/\pi$ & $\theta$&
$\langle{v}^{\prime\prime}\rangle$
& $Co$ & $Ro$ & $Re$ & $Ta$ & $\rm\Delta_u$ & $\rm\Delta_p$ (PSH)\\
\hline
R1 & 0    & 0   & $0^o$    & 0.074 & $0.000$& $\infty$ & 1088.883 &  0.00E+0 & 0.055& 0.384 \\
R2 & 0.25 & 1/8 & $22.5^o$ & 0.059 & 4.231  & 0.236 & 894.394 &  5.728E+7 &0.057& 0.399  \\
R3 & 0.25 & 1/4 & $45^o$   & 0.064 & 3.914  & 0.255 & 947.021 &  5.502E+7 &0.057& 0.399  \\
R4 & 0.25 & 3/8 & $67.5^o$ & 0.070 & 3.579  & 0.279 & 1022.810 &  5.360E+7 &0.051& 0.354 \\
R5 & 0.25 & 1/2 & $90^o$   & 0.066 & 3.776  & 0.265 & 957.403 &  5.227E+7 &0.031& 0.210  \\
R6 & 0.50 & 1/4 & $22.5^o$ & 0.056 & 8.980  & 0.111 & 840.540 &  2.279E+8 &0.060& 0.422 \\
R7 & 0.50 & 1/4 & $45^o$   & 0.059 & 8.432  & 0.119 & 879.430 &  2.200E+8 &0.055& 0.385 \\
R8 & 0.50 & 1/4 & $67.5^o$ & 0.067 & 7.411  & 0.135 & 975.684 &  2.091E+8 &0.052& 0.362 \\
R9 & 0.50 & 1/2 & $90^o$   & 0.064 & 7.795  & 0.128 & 875.589 &  1.864E+8 &0.022& 0.147 \\
R10& 0.25 & 0   & $0^o$    & 0.058 & 4.278  & 0.234 & 889.006 &  5.787E+7 &0.058& 0.407 \\
R11& 0.50 & 0   & $0^o$    & 0.054 & 9.283  & 0.108 & 822.673 &  2.331E+8 & 0.063& 0.445 \\
R12& 1.0  & 1/4 & $45^o$   & 0.069 & 14.474 & 0.069 & 1011.623 &  8.576E+8 &0.052& 0.363 \\
R13& 1.0  & 0   & $0^o$    & 0.050 & 20.004 & 0.050 & 768.311 &  9.448E+8 &0.063& 0.446 \\
R14& 0.50 & 1/4 & $45^o$   & 0.036 & 13.714 & 0.073 & 876.052 &  5.774E+8 & 0.013& 0.111 \\
R15& 0.50 & 1/4 & $45^o$   & 0.026 & 19.013 & 0.053 & 895.473 &  1.160E+9 &0.006& 0.061 \\
R16& 0.50 & 1/4 & $45^o$   & 0.019 & 25.469 & 0.039 & 799.152 &   1.659E+8 &0.005& 0.058\\
R17& 0.50 & 1/4 & $45^o$   & 0.062 & 7.971  & 0.126 & 1107.676 & 3.118E+8 & 0.049& 0.340 \\
\hline
\end{tabular}
\end{table}

\begin{figure}
\centering{\includegraphics[width=8cm]{fig2new.eps}}
\caption{Distribution of time- and horizontally averaged kinetic
energy fluxes with height for models (a) R1, R10, R11 and R13. The
corresponding values of angular rotational velocity
({\boldmath$\rm\Omega$}) are 0, 0.25, 0.50 and 1.0, respectively.
For these four cases the angle ($\theta$) between the rotation
vector and the vertical axis is $0\,^{\circ}$, (b) R10, R2, R3, R4
and R5. All these models have {\boldmath$\rm\Omega$} fixed at 0.25
while $\theta$ changes from $0\,^{\circ}$ to $90\,^{\circ}$ in steps
of $22.5\,^{\circ}$, (c) R7, R14, R15 and R16. The polytropic
indices in the upper stable layer for these four models are,
respectively, 2, 3, 4 and 5. The corresponding relative stability
parameters are 1, 3, 5 and 7, respectively. All the four models have
{$\mathbf{\Omega}=0.5$}  and $\theta=45\,^{\circ}$. \label{fig.2}}
\end{figure}

\begin{figure}
\centering{\includegraphics[width=10.9cm]{fig3new.eps}} \caption{(a)
Horizontally averaged kinetic energy flux as a function of time.
Total of 1,00,000 profiles were collected out of which 500 (every
200th) are plotted for case R1.  The height is denoted in the grid
points starting from one at the bottom (b) Penetration height
($\rm\Delta_u$) with time for model R1 for 1,00,000 time steps
(computed at every 10th time step). \label{fig3}}
\end{figure}

\section{Results and Discussion}
All 17 models computed in this study have a two layer configuration
in which the stable region is placed above the convectively unstable
region. Model R1 is non-rotating while models R2 to R17 are all
rotating, differentiated by a set of values of rate of rotation
({\boldmath$\rm\Omega$}), the angle between the $Z$-axis and the
rotation vector ($\theta$) and the polytropic index (or stability)
of the upper stable layer. Tables 1 and 2 list the physical
parameters of all the models.

The interface of the unstable-stable layer is located at a height of
0.6 from the bottom. This implies that the thickness of the lower
stable layer is 0.6 or 60\% of the total domain of computation while
the upper stable layer has a thickness of 0.4 corresponding to 40\%
of the total domain. The total domain of computation contains about
8.44 pressure scale heights (PSH) for models R1-R13 and R17 with the
upper stable layer containing 6.04 PSH while the lower convective
layer contains 2.40 PSH. For models R14, R15 and R16 the total
domain has around 9.89, 11.08 and 12.09 PSH, respectively with the
upper stable layer containing around 7.12, 8.01 and 8.76 PSH,
respectively.

It may be noticed that although the upper stable layer occupies only
40\% of the domain, it contains more than 70\% of the PSH. This is
because the length of the PSH decreases as we go from bottom to the
top. The density contrast ($\rho_{bottom}/\rho_{top}$) of models
R1-R13 is 238.0 while the temperature contrast
($T_{bottom}/T_{top}$) is 19.5 (cf. Table 1). For model R16 the
density contrast is the maximum (11024).

After the fluid has thermally relaxed, it is further evolved for
another $1,00,000$ time steps and the time and horizontal averages
are taken for quantities of interest. Figure 1 shows the
distribution of various energy fluxes with height for a number of
models. The calculation of the extent of the penetration into the
upper stable layer has been based on time- and horizontally averaged
kinetic flux ($\rm{F_k}$). Such a choice is obvious as the kinetic
flux is directly related with the motions and the profile is also
convenient for estimation of the extent of penetration (Hurlburt et
al. 1986, Hurlburt et al. 1994, Singh et al. 1994,1995, Saikia et
al. 2000, Pal et al. 2007). We illustrate this point by plotting the
distribution of kinetic energy flux for several sets of models in
Fig. 2. The kinetic energy flux is negative in most of the
convectively unstable region. Near the interface in the unstable
region it increases and becomes positive and then falls to zero in
the upper stable region. Ideally, the extent of penetration should
correspond to the first zero of the kinetic energy flux in the
stable layer above the unstable-stable interface. However, we define
the penetration height ($\rm\Delta_u$) to be the distance from the
interface of the unstable-upper stable layer to where $\rm{F_k}$ has
fallen to a value of $0.0001$ in the upper stable layer. The extent
of penetration distance ($\rm\Delta_u$) and penetration distance in
PSH ($\rm\Delta_p$) for all the 17 models are given in Table 3.

In Fig. 3(a), we show the horizontally averaged kinetic energy flux
as a function of depth and time for model R1.  We have plotted $500$
profiles of $\rm{F_k}$ which  corresponds to every 200th time step.
In Fig. 3(b), we show a time series of penetration distance
$\rm\Delta_u$ calculated by using the criterion outlined above. The
mean penetration is  indicated by the horizontal dashed line.

We now describe the effect of varying various parameters, e.g.,
{\boldmath$\rm\Omega$}$,\, \theta$, $S$, and the horizontal
resolution on the penetration distance above the convectively
unstable region.

\begin{figure}
\centering{\includegraphics[width=10.9cm]{fig4new.eps}}
\caption{Dependence of penetration height in PSH ($\rm\Delta_p$)
with (a) {\boldmath$\rm\Omega$} varying from 0, 0.25, 0.50 and 1.0
for the models R1, R10, R11 and R13, (b) colatitude ($\theta$)
varying from $0\,^{\circ}$  to $90\,^{\circ}$ in steps of
$22.5\,^{\circ}$ for models R10, R2, R3, R4 and R5 having
{$\mathbf{\Omega}=0.25$}. In this as well as in all the subsequent
figures, error bars denote the standard deviation in $\rm\Delta_p$
computed for $1,00,000$ time steps after the fluid has thermally
relaxed. \label{fig.4}}
\end{figure}

\subsection{Dependence of penetration height on \boldmath$\rm\Omega$}

Two sets of models R1, R10,  R11 and R13 and R3, R7 and R12 have
been computed to examine the effect of rate of rotation on the
penetration distance. In models R1, R10,  R11 and R13, the rotation
rate ({\boldmath$\rm\Omega$}) is systematically increased from 0 to
1.0 in steps of 0.25 (cf. Table 3). In all these four models, the
rotation vector ($\theta$) is kept at $0\,^\circ$, implying that the
rotation vector coincides with the $Z$-axis or the vertical
direction. Figure 4(a) shows the dependence of penetration distance
on {\boldmath$\rm\Omega$} for these four models. In model R1, which
is a non-rotating case ({\boldmath$\rm\Omega$}$\,=0$), the
penetration distance $\rm\Delta_u=0.055$ or 0.384 PSH. As the
rotation rate increases to 0.25 in model R10, the distance
$\rm\Delta_u$ increases to 0.058 which corresponds to 0.407 PSH. For
models R11 and R13 having {\boldmath$\rm\Omega$}$\,=0.5$ and 1.0
respectively, the penetration distance increases further to 0.445
PSH and 0.446 PSH respectively.

In the second set of models R3, R7 and R12, while the rotation rates
are, respectively, 0.25, 0.5 and 1.0, $\theta$ is changed to
$45\,^\circ$. We find that unlike in the first set of models (with
$\theta = 0\,^\circ$), the penetration distance decreases with
increase in the rotation rate. For example, the upward penetration
decreases from 0.399 PSH (Model R7, {\boldmath$\rm\Omega$} $= 0.25$)
to 0.363 PSH (Model R12, {\boldmath$\rm\Omega$} $= 1.0$). The
buoyancy-breaking is accentuated as rotation rate is increased when
$\theta$ is $45\,^\circ$.

The penetrative behaviour is markedly different from the case of
downward penetration (Pal et al. 2007), in which the penetration
distance ($\Delta_d$) decreased with increase in the rotation rate,
for any fixed value of colatitude $\theta$.

\begin{figure}
\centering{\includegraphics[width=8cm]{fig5new.eps}}
\caption{Variation of penetration height with Rossby number for (a)
three models R10, R11 and R13 having $\mathbf{\Omega}=0.25,\,0.50$
and $1.0$, respectively, and $\theta=0\,^\circ$, and (b) three
models R3, R7 and R12 having $\mathbf{\Omega}=0.25,\,0.50$ and
$1.0$, respectively, and $\theta=45\,^\circ$ The dashed lines
represent scaling laws: (a) $Ro^{-0.059}$ and (b) $Ro^{0.072}$.
\label{fig.5}}
\end{figure}

\begin{figure}
\centering{\includegraphics[width=8cm]{fig6new.eps}} \caption{Plot
of $(\rm{F_b}/\rho_{ctop})^{1/3}/\Delta_p$ for different $\Omega$
for (a) cases R10, R11 and R13 having $\theta=0\,^\circ$ and (b)
cases R3, R7, R12 having $\theta=45\,^\circ$. The dashed line shows
the scaling relationship (a) $1/\Omega^{-0.065}$ and (b)
$1/\Omega^{0.068}$ for the two sets of models. \label{fig.6}}
\end{figure}

For the two sets of models with $\theta=0\,^\circ$ and
$\theta=45\,^\circ$, we examine the behaviour of pentration height
with Rossby number. Figure 5(a) shows the variation of penetration
height in PSH ($\rm\Delta_p$) with Rossby number ($Ro$) for three
cases R10, R11 and R13 having $\theta=0\,^\circ$ and
$\mathbf{\Omega}=0.25,\,0.50$ and $1.0$, respectively. A smaller
Rossby number implies a larger rotation rate and the penetration
distance increases with decrease in the Rossby number from $0.407$
PSH for model R10 to $0.446$ PSH for model R13. For this set of
models with $\theta=0\,^\circ$, we see a relationship $\rm\Delta_p
\sim 0.38\, Ro^{-0.059}$.

Fig. 5(b) shows the variation of $\rm\Delta_p$ with $Ro$ for the
second set of models R3, R7 and R12 having $\theta=45\,^\circ$ and
$\rm\Omega\,=0.25,\,0.50$ and $1.0$, respectively. For this set, the
penetration distance decreases with the decrease in Rossby number.
Model R3 with  $Ro = 0.255$ has $\Delta_p = 0.399$ PSH, while R12
with $Ro = 0.069$ has $\Delta_p = 0.363$ PSH. This set of models
shows a scaling relation of $\rm\Delta_p \sim 0.44\, Ro^{0.072}$.

\subsection{Dependence of penetration height on colatitude $\theta$}

A set of five models R10, R2, R3, R4 and R5 has been examined to
study the effect of co-latitude ($\theta$) on the penetration
distance. All the models have {\boldmath$\rm\Omega$}$\,=0.25$ while
$\theta$ is changed from $0\,^\circ$ (Model R10) to $90\,^\circ$
(Model R5) in steps of $22.5\,^\circ$. Figure 4(b) shows the plot of
$\rm\Delta_p$ against $\theta$ for these five models. At
$\theta=0\,^\circ$, the penetration distance is maximum $0.407$ PSH.
It decreases slightly for successive values of $\theta$, and for
$\theta=90\,^\circ$ (case R5) has the minimum value of $0.210$ PSH.

Singh et al. (1994) computed many non-rotating cases having
different input fluxes and found a scaling relationship between the
upward penetration and the input flux ($F_b$) of the form $\Delta_p
\sim\, (\rm{F_b}/\rho_{ctop})^{1/3}$, where $\rho_{ctop}$ is the
density at the top of the convection zone. We attempt to find a
similar relationship which includes the effects of rotation.

In Fig. 6(a), we have plotted
($\rm{F_b}/\rho_{ctop})^{1/3}/\Delta_p$ for
{\boldmath$\rm\Omega$}$\,=0.25$, $0.5$, and $1.0$ corresponding to
models R10, R11, and R13, respectively. We find a scaling
relationship $\Delta_p \sim\,
(\rm{F_b}/\rho_{ctop})^{1/3}/0.283(\Omega)^{-0.065}$ for the three
cases having $\theta=0\,^\circ $. For cases R3, R7, R12 having
$\theta=45\,^\circ$ (Fig. 6(b)), we have a relationship of the form
$\Delta_p \sim\, (\rm{F_b}/\rho_{ctop})^{1/3}/0.35(\Omega)^{0.068}$.

For a set of models with a particular value of $\theta$, we have

\begin{eqnarray}
\Delta_p = c_1Ro^{x_1} \,\,\,\,\,\ \rm{and} \,\,\,\,\,\, \Delta_p =
\frac{(F_b/\rho_{ctop})^{1/3}}{c_2\Omega^{x_2}},
\end{eqnarray}
where $c_1$ and $c_2$ are constants and $x_1$ and $x_2$ are the
scaling powers for $Ro$ and $1/\Omega$.

One can use relations (4) to compute back the power laws of Fig. 6
by writing

\begin{eqnarray}
\Omega^{x_2} = \frac{(F_b/\rho_{ctop})^{1/3}}{c_1c_2Ro^{x_1}},
\end{eqnarray}
\begin{eqnarray}
x_2 =
\frac{\ln\frac{(F_b/\rho_{ctop})^{1/3}}{c_1c_2Ro^{x_1}}}{\ln\Omega}.
\end{eqnarray}

Substituting the value of $Ro,\, c_1,\,c_2,\,x_1, \rm{and}\,\,
\Omega\,$ in Eqn.(6), we find that for cases with
$\theta=0\,^{\circ}$, $x_2$ to be $-0.063$ and for cases with
$\theta=45\,^{\circ}$ it is $0.062$ which are close to the values
$-0.065$ and $0.068$, respectively, as shown in Fig. 6.

\begin{figure}
\centering{\includegraphics[width=8cm]{fig7new.eps}}
\caption{Instantaneous vertical velocity at a horizontal surface at
unstable-upper stable layer boundary located at a height of 0.6 from
the bottom for two cases R7 ($\mathbf{\Omega}=0.5,\,
\theta=45\,^\circ$) and R16 ($\mathbf{\Omega}=0.5,\,
\theta=45\,^\circ$). Polytropic indices of the upper stable layer
are 2 (Model R7) and 5 (Model R16). \label{fig.7}}
\end{figure}

\subsection{Dependence of penetration height on relative stability parameter $S$}

We have constructed four models with different stability parameter
($S$) defined in Eqn. (3) by varying the polytropic index of the
upper stable layer. We have $S = 1,\,3,\,5$ and $7$ for cases R7,
R14, R15 and R16, respectively. In Fig. 2(c), we have plotted the
distribution of the kinetic energy flux with height for these four
models. The values of the angular velocity ({\boldmath$\rm\Omega$})
and angle between the rotation vector and vertical axis ($\theta$)
are fixed at $0.5$ and $45\,^\circ$, respectively, for these models.
The corresponding penetration heights $\rm\Delta_u$ and
$\rm\Delta_p$ (in PSH) are given in the last two columns of Table 3.
Figure 2(c) clearly shows the differences in the kinetic flux
profile with $S$ especially in the region near the stable-unstable
layer interface. We find a decrease in penetration distance from
0.385 PSH in model R7 ($S=1$) to 0.058 PSH in model R16 ($S=7$).
Figure 7 shows the instantaneous vertical velocities at the
horizontal interface of the unstable-stable interface for models R7
and R16. The vertical velocities for model R16 with stiffer stable
layer (S=7) are much lower compared to that of model R7 (S=1).

Fig. 8 shows the variation of penetration distance with stability
factor $S$ for models R7, R14, R15, \& R16. These models show a
scaling relationship $\Delta_p \sim S^{-1}$ . This means that the
convection is efficient enough to establish a nearly adiabatic
stratification in the upper stable layer and the penetration is
nearly adiabatic (Hurlburt et al. 1994, Zahn 2002). Brummel et al.
(2002) and Pal et al. (2007) found a scaling relationship of
$\Delta_p \sim S^{-1/4}$ for downward penetration for rotating
convection zones associated with a thermal adjustment region at the
bottom.

\begin{figure}
\centering{\includegraphics[width=8cm]{fig8new.eps}}
\caption{Variation of penetration distance $\rm\Delta_p$ with
stability parameter $S$ for models R7, R14, R15 and R16. Here,
$\mathbf{\Omega}=0.5$ and $\theta=45\,^{\circ}$. The dashed line
represents scaling law, $S^{-1}$. \label{fig.8}}
\end{figure}
\begin{table}
\caption{Parameters of models computed to examine the effect of
resolution on the upward penetration}\label{table:4} \centering
\begin{tabular}{ccccccc}
\hline\hline Models & Grid &{\boldmath$\rm\Omega$}  & $\theta/\pi$ &
$\theta$ &$t$ & $\Delta_p$ (PSH)\\
\hline
R1A & $64 \times 64 \times 128 $& 0    & 0   & $0^o$    &1056 & 0.303\\
R7A & $64 \times 64 \times 128$ &0.50 & 1/4 & $45^o$   & 907 & 0.333\\
R7B & $96 \times 96 \times 96$ &0.50 & 1/4 & $45^o$   & 550 & 0.268
\\ \hline
\end{tabular}
\end{table}

\subsection{Dependence of penetration height on resolution}

Table 4 lists additional models computed to see the effect of
resolution on the upward penetration. Model R1A is same as model R1
except that it has a total of $524,288\,(64 \times 64 \times 128)$
grid points instead of $117,600\,(35 \times 35 \times 96)$ grid
points, having better horizontal as well as vertical resolution.
Penetration height for non-rotating model R1 is $0.384$ PSH (cf.
Table 3), while for the non-rotating model R1A with increased
resolution, it is $0.303$ PSH.

Two rotating models R7 and R7A examine the effect of change in both
horizontal and vertical resolution on the penetration distance. We
notice a decrease in penetration depth from $0.385$ PSH in model R7
($35 \times 35 \times 96$) to $0.333$ PSH in model R7A ($64 \times
64 \times 128$). It may be noticed that while models R1 and R7 had
almost same penetration distance, the decrease in $\Delta_p$ with
increased resolution is smaller for model R7A which has a higher
rotation rate compared to model R1A which is non-rotating. It seems
that the effect of increase in resolution on penetration is
different for models with different rotation rates.

Three models R7, R17, and R7B examine the effect of horizontal
resolution on the penetration distance above the convection zone.
While all the models have 96 grid points in the vertical, they
differ in their horizontal resolution. Model R7 has $35 \times 35$
grid points in the horizontal direction, model R17 has $46 \times
46$ and model R7B has $96 \times 96$ grid points. The penetration
distance decreases with the increase in horizontal resolution from
$0.385$ PSH (Model R7) to $0.340$ PSH (Model R17) to $0.268$ PSH
(Model R7B).

\section{Conclusions}

We have computed $17$ models and presented the results of our
three-dimensional numerical simulations of turbulent compressible
convection penetrating into a radiative envelope under the influence
of rotation. We find that the penetration distances for these models
lie in the range $0.058\,\rm{PSH} \le \rm\Delta_p \le 0.446\,
\rm{PSH}$.

Recently, Woo \& Demarque (2001) put an empirical constraint on the
convective core overshoot for intermediate-to-low mass stars by
using Roxburgh's integral constraint. They found that the proper
limit of core overshoot for these stars would be $15\%$ of the core
radius. In Table 5, we have shown this calculation for our 17
models. We have computed the penetration distance in PSH as a
percentage of the total size of the convection zone. The size of the
convection zone is taken from the last column of Table 1 while the
penetration height ($\Delta_p$) has been taken from the last column
of Table 3. As can be seen from Table 5, We find an upper limit on
the penetration into the upper radiative layer to be around
$18.5\%$.

\begin{table}
\caption{Penetration as a percentage of the size of the Convection
Zone }\label{table:5} \centering
\begin{tabular}{cccccccccccccc}
\hline\hline Model& R1 & R2 & R3 & R4 & R5 & R6 & R7 & R8 & R9 & R10
& R11 & R12 &R13  \\
 $\%\Delta_p$ & 16.0
& 16.6 & 16.6 & 14.8 & 8.8 & 17.6 & 16.0 & 15.0 & 6.0 & 17.0 & 18.5
& 15.1 & 18.6 \\
Model & R14 & R15 & R16 & R17 & & & & & & & & &  \\
$\%\Delta_p$ & 4.0 & 2.0 & 1.7 & 14.2 & & & & & & & & &   \\
\hline
\end{tabular}
\end{table}

For rotation about a vertical axis ($\theta=0\,^\circ$), the
penetration into the upper stable region increases as the angular
rotational velocity ({\boldmath$\rm\Omega$}) increases or the Rossby
number ($Ro$) decreases.  However, this trend is reversed for
rotation around an inclined axis. When the angle
($\theta=45\,^\circ$), the penetration distance into the radiative
envelope  decreases with increasing {\boldmath$\rm\Omega$} owing to
horizontal mixing. To see the effect of change of angle of
inclination on penetration, angle $\theta$ is systematically varied
from $0\,^\circ$ to $90\,^\circ$ in steps of $22.5\,^\circ$ for five
models with a fixed angular velocity {\boldmath$\rm\Omega$}$=0.25$.
We again find that the penetration distance decreases as the
colatitude $\theta$ is increased. We do not see this behaviour
changing even when the resolution is increased.

We also find that the penetration distance above the convection zone
obeys a scaling relation of the form $\Delta_p \sim S^{-1}$
appropriate for nearly adiabatic penetration even in the presence of
rotation. The present simulations need to  be extended to include
more realistic input physics and a higher resolution to enable us to
get a better insight into the dynamics of rotating convection near
the convective-radiative interface. Since the motivation of the
study was to look into the effect of varying rotation rates on the
penetration distance, we needed to compute several models. Due to
the restrictions in speed and memory, we were unable to increase the
resolution of all the models. We hope to achieve this in future
studies.

\begin{acknowledgements}
      PSP is grateful to University Grants Commission, India for a
      Senior Research Fellowship. HPS was supported by grants from
      Indian Space Research Organization and Deutscher Akademisher Austauschdienst
      (DAAD), Germany. KLC was partly supported by
      grants from Hong Kong Research Grant Council. Authors acknowledge
the use of High Performance Computing facilities at IUCAA, Pune,
India.
\end{acknowledgements}


\section{Structural parameters of the bulges and disks}

The structural parameters of the bulges and disks of a
magnitude-limited sample of 148 S0-Sb galaxies were derived from 2MASS
$J-$band images. We used GASP2D, our new code for two-dimensional
photometric decomposition (M\'endez-Abreu et al. 2008).
The S\'ersic bulge and exponential disk were assumed to be
characterized by elliptical and concentric isophotes with constant
(but possibly different) ellipticity and position angle.
We found that the surface-brightness radial profiles of larger bulges
are more centrally peaked than those of the smaller bulges. Larger
bulges have a lower effective surface brightness and higher absolute
luminosities. They reside in larger disks, as revealed by the
correlation between central velocity dispersion and disk scalelenght.
This reveals a strong coupling between bulges and disks. Larger disks
have a lower central surface brightness.

\section{The equatorial intrinsic ellipticity of bulges}

The determination of the intrinsic shape of bulges is an ill-posed
problem, that can be tackled only in a statistical way. A new approach
to derive the PDF of intrinsic equatorial ellipticities using only
photometric data was used (M\'endez-Abreu et al. 2008). Fig.~1 shows
the PDF obtained from our galaxy sample. The significant decrease of
probability for $E<0.07$ (or equivalently $B/A<0.93$) suggests that
the shape of bulges is elliptical rather than circular.
The average ellipticity $\langle E \rangle = 0.16$ ($\langle B/A
\rangle = 0.85$) is in agreement with previous findings by Bertola et
al. (1991) and Fathi \& Peletier (2003). The
$\rm{PDF}(E)$ is not related to morphology, light concentration, or
luminosity of bulges.

\begin{figure}
  \plotone{mendezabreu_a1.ps}
  \caption{Left panel: \rm{PDF} of our observable $T$ which is
    calculated from the twist angle ($\delta$) between the bulge and
    disk, bulge ellipticity ($\epsilon$) and galaxy inclination
    ($\theta$). Right panel: \rm{PDF} of $E$ obtained from the
    inversion of $T$. $E$ is related to the intrinsic axis ratio by
    $E$=$(1-B/A)/(1+B/A)$ where $A$ and $B$ are the semi-axis of the
    bulge ellipsoid. The probability is normalized over 10 bins
    geometrically distributed, the width of the first bin is 0.03 and
    the width ratio is 1.25. Error bars correspond to a Poisson
    statistics.}
\end{figure}

\section{Discussion and conclusions}

The photometric and kinematic parameters of the bulges of our sample
galaxies follow the same fundamental and photometric planes of
elliptical galaxies, supporting the idea that both are formed in the
same way. However, the correlations between the bulge and disk
structural parameters are usually interpreted as an indication of
secular evolution processes. Therefore, the above relations are not
enough to distinguish between bulges formed by early dissipative
collapse, merging or secular evolution. All these scenarios could be
tested against the PDF of the intrinsic equatorial ellipticity of
bulges.


\section{Introduction}

It is now possible to make timing observations of millisecond pulsars to a precision of $\sim 100$\,ns.  Such exquisite timing precision allows the pulsar astrometric, spin and orbital parameters to be determined with great accuracy.  This has recently allowed us to limit the rate of change of the gravitational constant, $|\dot{G}/G| < 2.1\times10^{-12}$\,yr$^{-1}$, and to obtain the most accurate astrophysical distance estimate outside of the Solar System (Verbiest et al. 2007, submitted to ApJ).  Multi-frequency pulsar observations have allowed us to study the Solar corona (You et al. 2007, in press) and density fluctuations in the interstellar plasma (You et al. 2007b)\nocite{yhc+07}. One of the most exciting applications of such data-sets is to search for the signatures of gravitational waves (GWs) passing over the Earth.    This is the main goal of the Parkes Pulsar Timing Array (PPTA) project\footnote{http://www.atnf.csiro.au/research/pulsar/ppta}, which aims to observe 20 millisecond pulsars with a timing precision close to 100\,ns over more than five years.

\begin{figure}
\begin{center}\includegraphics[width=8cm,angle=-90]{sens.ps}\end{center}
\caption{Characteristic strain sensitivity for existing and proposed GW detectors as a function of GW frequency. Predicted signal levels from various astrophysical sources are shown.}\label{fg:sens}
\end{figure}

The existence of a GW signal can be identified by analysis of pulsar timing residuals.  In brief, the pulsar timing model (see, e.g. Edwards, Hobbs \& Manchester 2006)\nocite{ehm06} allows the observed arrival times of pulses from a radio pulsar to be compared with a model of the pulsar's astrometric, orbital and rotational parameters.  The differences between the actual arrival times and the predicted arrival times are the ``pulsar timing residuals".  Any features observed in the timing residuals indicate the presence of unmodelled effects which may include calibration errors, (additional) orbital companions, spin-down irregularities, or GWs. The induced pulsar timing residuals due to a GW signal were first calculated by Sazhin (1978)\nocite{saz78} and Detweiler (1979)\nocite{det79}.  They showed that a GW signal causes a fluctuation in the observed pulse frequency $\delta \nu/\nu$ which affects pulsar timing residuals at time $t$ from the initial observation as
\begin{equation}
R(t) = -\int_0^t\frac{\delta \nu(t)}{\nu}dt.
\end{equation}
For a given pulsar and GW source, this effect is only dependent upon the characteristic strain at the pulsar and at the Earth. The GW strains evaluated at the positions of multiple pulsars will be uncorrelated, whereas the component at the Earth will lead to a correlated signal in the timing residuals of all pulsars. As pulsars are typically observed every few weeks for many years and standard pulsar timing techniques absorb any low-frequency GWs by fitting for the pulsar's spin-down\footnote{As the intrinsic pulsar spin-period and rate of spin-down is unknown, a quadratic curve is always fitted and subsequently removed from the pulsar timing residuals.}, pulsar timing experiments are sensitive to GW signals in the ultra low frequency ($f \sim 10^{-9}$\,Hz) band. This technique is therefore complementary to other GW detection methods such as the Laser Interferometer Space Antenna (LISA) and ground based interferometer systems (such as LIGO and VIRGO), which are sensitive to higher frequency GWs ($f \sim 10^{-2.5}$ and $f\sim 10^{2.5}$\,Hz respectively). In Figure~\ref{fg:sens} we plot the sensitivity of the pulsar timing experiments, LISA and LIGO and the expected GW sources.  For the pulsar timing experiments, maximum sensitivity is reached at $f = 1/T$ where $T$ is the data span. As described above, the pulsar experiments are not sensitive to GWs with $f < 1/T$ and, for white pulsar timing residuals, the upper bound on the characteristic strain increases as $f^{3/2}$ for $f > 1/T$.

In this paper we review the potential sources of GWs (\S\ref{sec:sources}) before describing how upper bounds on the amplitude of any existing background can be calculated and highlighting the astrophysical implications of these bounds. In \S\ref{sec:detect} we describe how a stochastic GW background may be detected.  Finally, we highlight some practical issues that need to be addressed in achieving our sensitivity goal (\S\ref{sec:prac}) and describe possible future experiments (\S\ref{sec:future}).

\section{Potential sources of gravitational waves}\label{sec:sources}

Theoretical models predict that pulsar timing experiments will be sensitivity to burst GW sources, individual sources and to a stochastic background.

\subsection{Single sources}

Recently Sudou et al. (2003)\nocite{simt03} reported the possible discovery of a binary supermassive black hole system in the radio galaxy 3C66B.  As shown by Jenet et al. (2004)\nocite{jllw04}, the GW radiation from this postulated system would have produced clearly detectable fluctuations in the timing residuals of existing pulsar data-sets.  The non-detection of the expected GW signal (see Figure~\ref{fg:simul}) ruled out the postulated system with 95\% confidence. However, detecting (or limiting the existence of) individual supermassive binary black hole systems is necessary for testing theories of hierarchical galaxy formation. Order-of-magnitude calculations demonstrate that a system with a chirp mass of $10^{9}$M$_\odot$ and an orbital period of 10\,yr, at a distance of 20\,Mpc will induce a sinusoidal signal in the timing residuals with amplitude $\sim 100$\,ns at a frequency of twice the orbital frequency.

\begin{figure}
\begin{center}\includegraphics[width=5cm,angle=-90]{3c66b_3.ps}\end{center}
\caption{Simulation of the induced timing residuals for PSR~B1855+09 due to the postulated binary black hole system in the radio galaxy 3C66B.}\label{fg:simul}
\end{figure}

\subsection{Burst sources}

Sources of burst GW emission that may be detectable with current and proposed pulsar timing experiments include 1) the formation of supermassive black holes (Thorne \& Braginskii 1976)\nocite{tb76}, which lead to a day-long burst of radiation, 2) highly eccentric supermassive black hole binaries (Enoki \& Nagashima 2007)\nocite{en07}, 3) close encounters of massive objects (Kocsis et al. 2006)\nocite{kgm06} and 4) cosmic string cusps (Damour \& Vilenkin 2001)\nocite{dv01}.  The sensitivity of a pulsar timing array experiment to burst GW emission depends upon the number of pulsars in the array, the timing precision, number of observations, the pulsar positions on the sky and the coordinates of the GW burst source.  With our final data-sets we expect to achieve an angular resolution of approximately 30$^\circ$ (Lommen, private communication). By sharing data with timing array projects being carried out in the Northern hemisphere we expect to improve both the angular resolution of our detector and our sensitivity to burst sources in the Northern Hemisphere.

\subsection{Stochastic background}

A stochastic background of GWs due to binary supermassive black holes (Jaffe \& Backer 2003\nocite{jb03}, Wyithe \& Loeb 2003\nocite{wl03}, Enoki et al. 2004\nocite{eins04}, Sesana et al. 2004\nocite{shmv04}), cosmic strings or relic GWs from the big bang (Maggiore 2000\nocite{mag00}) are all potentially detectable.  In most models, the GW strain spectrum,  $h_c(f)$, can be represented by a power-law in the GW frequency, $f$,
\begin{equation}
h_c(f) = A\left(\frac{f}{\rm yr^{-1}}\right)^\alpha
\end{equation}
where the spectral exponent $\alpha = -2/3$, $-1$ and $-7/6$ for likely GW backgrounds due to coalescing black hole binaries, cosmic strings and relic GWs respectively (Figure~\ref{fg:a_upper}).  The energy density of the background per unit logarithmic frequency interval can be written as
\begin{equation}
\Omega_{\rm GW}(f) = \frac{2}{3}\frac{\pi^2}{H_0^2}f^2h_c(f)^2
\end{equation}
where $H_0$ is the Hubble constant.

\begin{figure}
\begin{center}\includegraphics[width=5cm,angle=-90]{a_upper.ps}\end{center}
\caption{The solid line indicates our current limits on GW backgrounds of specified spectral exponent.  The dashed line underneath indicates the expected limits within 5 years.  The near vertical lines give predicted background amplitudes.  The ``star"-symbol indicates the earlier limit obtained by Kaspi et al. (1994).}\label{fg:a_upper}
\end{figure}

Jenet et al. (2006)\nocite{jhs+06} introduced a method to place an upper bound on $A$ for a given $\alpha$ using actual pulsar timing residuals. In contrast with earlier techniques (e.g.  Kaspi et al. 1994\nocite{ktr+94}) this method used the \textsc{tempo2} software package (Hobbs, Edwards \& Manchester 2006)\nocite{hem06} to account for all the fitting procedures undertaken during the pulsar timing process.  Using observations from the PPTA project along with archival data from the Arecibo telescope, Jenet et al. (2006)\nocite{jhs+06} obtained an upper limit on the energy density per unit logarithmic frequency interval of $\Omega_{\rm GW}[1/(8{\rm yr})]h^2 \leq 2\times10^{-8}$ for $\alpha=-2/3$, which corresponds to $A \leq10^{-14}$ (see Figure~\ref{fg:a_upper} and the ``Current limit'' in Figure~\ref{fg:sens}).  This is the most stringent limit to date on the existence of a GW background in the nano-Hertz frequency band and has astrophysical and cosmological implications.  For instance, this limit was used to 1) constrain the merger rate of massive black hole systems at high redshift, 2) rule out some relationships between the black hole mass and galactic halo mass, 3) constrain the rate of expansion in the inflationary era and 4) provide an upper bound on the dimensionless tension of cosmic strings.

We have recently shown that the Jenet et al. (2006)\nocite{jhs+06} technique described above suffers from two limitations.  First, it does not give the lowest possible upper bound and secondly, it relies on the observed spectrum being white.  Unfortunately, even for millisecond pulsars, it is unusual for the timing residuals not to be affected by unexplained phenomena, such as pulsar ``timing noise" (see e.g. Hobbs, Lyne \& Kramer 2006)\nocite{hlk06}.  A new method has now been developed and will be described in a subsequent paper.  Our new method is applicable to any pulsar timing data-set and will lead to significantly improved bounds on the GW background amplitude.

\section{Detecting a background}\label{sec:detect}

Jenet et al. (2005)\nocite{jhlm05} developed a technique for making a definitive detection of an isotropic, stochastic, background of GWs by searching for correlated signals between pulsar data-sets.  For each pair of pulsars, the zero-lag correlation between the respective timing residuals can be calculated.  The expected correlation as a function of angle between the pulsars for timing residuals dominated by a GW background was first calculated by Hellings \& Downs (1983)\nocite{hd83} and is shown in Figure~\ref{fg:corr}\footnote{This correlation curve assumes the general theory of relativity. We are calculating the equivalent curves for other theories of gravity.  These curves and their implications will be discussed elsewhere.}.  Jenet et al. (2005)\nocite{jhlm05} showed, for expected amplitudes of a GW background, that this signal could be unambiguously detected if 20 or more pulsars were observed over a period of 5 years each with an rms timing residual of ~100---500\,ns.

\begin{figure}
\begin{center}\includegraphics[width=7cm,angle=-90]{hdCurve3.ps}\end{center}
\caption{The expected correlation in the timing residuals of pairs of pulsars as a function of angular separation for an isotropic GW background.}\label{fg:corr}
\end{figure}

We have recently generalised the initial work of Jenet et al. (2005) to produce routines that can be applied to the observed pulsar timing residuals (with the actual sampling, data lengths etc.).  This new technique also increases the significance of any possible detection, by optimally pre-whitening the data-sets. An example of the improvement possible by pre-whitening the timing residuals is shown in Figure~\ref{fg:sig_a} where the significance of detection is plotted versus the GW background amplitude for the PPTA goal of timing 20 pulsars with 100\,ns timing precision over a period of five years. Without any pre-filtering of the data we will only achieve a maximum detection significance of $S \sim 3.5$ (corresponding to the sensitivity labelled as the ``PPTA limit'' in Figure~\ref{fg:sens}).  Using optimal whitening schemes this sensitivity can be significantly enhanced. Initial results suggest that if we can reduce our current timing precision for all pulsars by a factor of $\sim 2$ we will be sensitive, within five years, to the maximum amplitude predicted by models of an isotropic stochastic GW background caused by coalescing black-holes in the centres of galaxies\footnote{Note that this limit is similar to the current limits described in \S2.3 that were obtained using a few long data-sets.  We emphasise that, even though limits can be placed using a few data-sets, detecting a GW background requires observations of $\sim$20\,pulsars.}.

\begin{figure}
\begin{center}\includegraphics[width=10cm]{sig_a.eps}\end{center}
\caption{The sensitivity of the completed PPTA project to a GW background of specified amplitude. The three curves represent no pre-filtering of the data (lowest), using an optimal low-pass filtering technique (middle) and optimally pre-whitening the data-sets (highest).}\label{fg:sig_a}
\end{figure}

\section{Practical issues}\label{sec:prac}

Pulsar timing is affected by the stability of terrestrial clocks, ephemeris errors and the pulsars themselves.  The goal of the PPTA project, described above,  to obtain data-sets for 20 pulsars with approximately two-weekly sampling with rms timing residuals of $\sigma \sim 100$\,ns over five years, is challenging.  To date we have a few pulsars where we can obtain such precision. For instance, van Straten et al. (2001)\nocite{vbb+01} obtained $\sigma = 130$\,ns over 40 months of observing PSR~J0437$-$4715. For the majority of our pulsars $\sigma \sim 1\mu s$.  We have many possibilities for decreasing these timing residuals, including:

\begin{itemize}
\item{\emph{correcting the data-sets for dispersion measure variations}: You et al. (2007) showed that for PSR~J1939+2134, $\sigma$ decreased from 290\,ns to 190\,ns after correction.}
\item{\emph{improving instrumentation}: We have recently developed new digital filterbank systems for the Parkes telescope.  The recently commissioned system provides 1\,GHz of bandwidth at an observing frequency of 3\,GHz.  We are also developing a new wide-bandwidth coherent dedispersion system.}
\item{\emph{new processing techniques}:  We are exploring new methods to calibrate our data-sets which should lead to at least a factor of 2 improvement in our rms residuals.  We have recently produced a new pulsar timing package, \textsc{tempo2} (Hobbs, Edwards \& Manchester 2006)\nocite{hem06}, that  is accurate (for known physics) at the 1\,ns level.}
\end{itemize}

\section{Future possibilities}\label{sec:future}

The sensitivity of a pulsar timing array to a stochastic background is proportional to the average timing precision, the square root of the number of observations and to the number of pulsars in the array.
This paper has concentrated on the Parkes pulsar timing array project.  However, significant improvements will be made to the timing array sensitivity to GW sources by combining our data-sets with Northern Hemisphere pulsar observations.  The North American pulsar timing array project (NanoGrav) has been observing a sample of pulsars using the Arecibo and GreenBank observatories.  These observations provide more than 20\,years of timing PSRs B1855+09 and B1937+21 and were used to provide the first stringent limits on the stochastic gravitational wave background (e.g. Kaspi et al. 1994) and to rule out the postulated binary black hole system in the radio galaxy 3C66B (Jenet et al. 2004).  The European Pulsar Timing Array (EPTA) team have recently started to obtain data on a large sample of pulsars using the four major European telescopes at Jodrell Bank, Effelsberg, Nancay and Westerbork.  We are currently developing techniques to combine observations from different observatories and will soon produce new limits on the existence of a GW background using all available data.

Figure~\ref{fg:sens} shows the sensitivity that should be achievable using the Square Kilometer Array (SKA) telescope hoped to be built by the year 2020.  This sensitivity curve assumes an SKA timing array project that could observe 100 millisecond pulsars with rms timing residuals around 50\,ns for a period of 10\,yr. Provided intrinsic pulsar timing noise does not dominate the effect of the stochastic GW background then this limit corresponds to a detection limit at 3nHz of $\Omega_{gw} \sim 10^{-13}$.  Pre-whitening procedures have the potential to greatly decrease this limit.

\section{Conclusion}

It is possible that gravitational waves will be detected within the next decade by world-wide pulsar timing array projects.  As a by-product of these investigations stringent checks will also be placed on terrestrial time standards and the solar system ephemeris.  Regular dual-frequency observations of multiple pulsars will also provide valuable information about the interstellar medium.  Using a pulsar array as a gravitational wave detector is complimentary to other searches that are attempting to detect much higher-frequency gravitational waves.

\subsection{Acknowledgments}

This work is undertaken as part of the Parkes Pulsar Timing Array project, which is a collaboration between the ATNF, Swinburne University and the University of Texas, Brownsville.  The author wishes to than R. Manchester, W. Coles and F. Jenet for helping to improve this paper. The Parkes radio telescope is part of the Australia Telescope, which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO.

\section*{References}

\bibliographystyle{aipproc}


\section{Introduction}
\indent

A black hole is an object with an extreme gravitational field
whose escape velocity exceeds the speed of light. Classical treatments of static,
spherically symmetric black holes generally use the Schwarzschild solution to
Einstein's equations of general relativity.  For astrophysical black holes which
likely undergo periods of accretion and evaporation,
it is more realistic to consider dynamic rather than static behavior.

We have examined the geometry of an example dynamic black hole that
evaporates at a steady rate, as measured by a distant observer.
A flat Minkowski space-time results once the mass
of the black hole vanishes\cite{JLOct07}.

The metric defining that space-time considered is given by
\begin{equation}
ds^2 = -\left (1-{R_M (ct) \over r} \right ) (dct) ^2 +
2 \sqrt{{R_M (ct) \over r}} \, dct \, \, dr + dr^2 + r^2 \, d \omega ^2
\label{metric}
\end{equation}
where $d\omega^2 \equiv d\theta^2 + sin^2 \theta \, d\phi^2$,
and it is assumed\cite{JLMay07} that $\ddot{R}_M =0$.
In this equation, $R_M \equiv 2 G_N M(ct) / c^2$ is a dynamic
form of the Schwarzschild radius that we refer to as the
\emph{radial mass scale}.  The metric in Eq. \ref{metric} is
seen to have anomalous behavior at $r=R_M (ct)$.

Photons are given by null geodesics. The  radial null surfaces ($ds^2 =0$)
for the metric specified in Eq. \ref{metric} satisfy
\begin{equation}
{d r_\gamma \over dct} = - \sqrt{R_M \over r_\gamma} \pm 1,
\label{photons}
\end{equation}
where the +/- sign describes outgoing/ingoing photons.
The radial coordinate corresponding to the horizon satisfies
\begin{equation}
{d R_H \over dct} = - \sqrt{R_M \over R_H} + 1 \quad , \quad
R_H = {R_M \over \left (1 - \dot{R}_H  \right ) ^2} < R_M \quad ,
\label{horizon}
\end{equation}
and is given by
the particular null surface proportional to the radial mass scale
$R_H=R_M / \zeta_H$, where
\begin{equation}
\dot{R}_M =
\zeta_H \, \left(1-\sqrt{\zeta_H} \right ) .
\label{zetaH}
\end{equation}
The ability to relate the horizon scale to the radial mass scale
in terms of a proportionality constant $\zeta_H$
allows this algebraic solution for a given rate of
evaporation.

We examine the radial motion of photons emitted near the black hole horizon.
It is shown that some photons emitted inside
the radial mass scale (or ``shrinking Schwarzschild radius")
can escape falling into the horizon, and that this behavior is clearly demonstrated
through the use of the Penrose diagram.

\setcounter{equation}{0}
\section{Behavior of Near-Horizon Photons
\label{mainsection}}

\subsection{Conformal coordinates for Penrose diagram}
\indent

The form of the conformal time and radial coordinates
that are used for the construction of the Penrose diagram
have been developed in a companion paper\cite{JLOct07}:
\begin{equation}
\begin{array}{c}
ct_* = {r \over 2} \left (
 exp \left [ \int ^ {R_M (ct) \over r} { \left (
1 + \sqrt{\zeta '}  \right ) d \zeta '  \over
\left \{ \zeta' \left ( 1 + \sqrt{\zeta'} \right )
+ \dot{R}_M \right \} }
\right ]  -
 exp \left [ \int ^ {R_M (ct) \over r} { \left (
1 - \sqrt{\zeta '}  \right ) d \zeta '  \over
\left \{ \zeta' \left ( 1 - \sqrt{\zeta'} \right )
- \dot{R}_M \right \} }
\right ] \right ) \\ \\
r_* = {r \over 2} \left (
 exp \left [ \int ^ {R_M (ct) \over r} { \left (
1 + \sqrt{\zeta '}  \right ) d \zeta '  \over
\left \{ \zeta' \left ( 1 + \sqrt{\zeta'} \right )
+ \dot{R}_M \right \} }
\right ]  +
 exp \left [ \int ^ {R_M (ct) \over r} { \left (
1 - \sqrt{\zeta '}  \right ) d \zeta '  \over
\left \{ \zeta' \left ( 1 - \sqrt{\zeta'} \right )
- \dot{R}_M \right \} }
\right ] \right )  .
\end{array}
\label{conformal}
\end{equation}
This equation relates the space-time coordinates of an
asymptotic observer $(ct, r)$ with the conformal coordinates
$(ct_* , r_*)$.

\subsection{Outgoing photon trajectories}
\indent

We'll first examine the trajectories of outgoing photons corresponding to
the + sign in Eqn \ref{photons}.
All photons are emitted at $ct = -10$ units from various
initial radial coordinates.  The calculated trajectories
are shown in space-time plots in Fig. \ref{Outgo}.
\setfigure{Outgo}{Trajectories of outgoing photons emitted
(a) inside $R_H$, (b) between $R_H$ and $R_M$, and
(c) at $R_M$.}
In this figure, the dashed line closest to the singularity
represents the horizon of the black hole $R_H (ct)$,
while the outer dashed line represents the radial
mass scale $R_M (ct)$ (the coordinate anomaly).
The photon labeled (a) starts inside the black hole horizon.
Although it is an outgoing photon,
the radial coordinate of the photon $r_\gamma$ decreases monotonically
until it reaches the singularity at $r=0$.

The photon labeled (b) starts between the horizon $R_H$ and the
radial mass scale $R_M$.
Initially, the photon falls in toward the singularity ($r =0$),
tracking (but remaining just outside) the horizon.
This outgoing photon falls inward due to the strong gravity near the black hole.
In the late stages of evaporation,
gravity weakens enough to allow the photon to move away from the horizon.
As $r_\gamma$ crosses
$R_M$, the photon is momentarily stationary in the $(ct,r)$ coordinates.
The slope of the photon
approaches 45$^\textnormal{\tiny{o}}$ at spatial infinity as it
escapes the influence of the black hole.

The photon labeled (c) starts on the coordinate anomaly, $R_M$ .
The photon is initially stationary at $R_M$,
and at no point does $r_\gamma$ decrease.
The coordinate anomaly shrinks away from the photon as the black hole evaporates
and gravity weakens, and the photon subsequently escapes.

The trajectories of these outgoing photons
are illustrated on the Penrose diagram in Fig. \ref{OutPen}.
\setfigure{OutPen}{Penrose diagram of an evaporating black
hole with outgoing photon trajectories over plotted.  Red curves
(running vertically in right hand region) represent curves of
constant $r$.  The blue-green curves represent curves of constant $ct$.}
In this Penrose diagram, the space-time structure is represented using
functions of the conformal coordinates Eq. \ref{conformal}
developed in the companion paper \cite{JLOct07}.
On the Penrose diagram, the path of photons are immediately apparent,
since all light-like surfaces lie at 45$^\textnormal{\tiny{o}}$ angles.

\subsection{Ingoing photon trajectories}
\indent

We will next examine the trajectories of ingoing photons.
For an evaporating black hole, there exists a (ingoing) \emph{light-like
boundary} that separates the regions of space-time for
which any object can communicate with the singularity.
This shrinking boundary, with radial coordinate
$R_B$, is shown as an additional dashed
line outside of the horizon and radial mass scale in the
space-time plot of Fig. \ref{Ingo}.
\setfigure{Ingo}{Trajectories of ingoing photons emitted
(d) at $R_M$, (e) between $R_M$ and $R_B$,
and (f) outside $R_B$.}
As with the outgoing photons, all ingoing photons are
emitted at ct=-10 units.
The photon labeled (d) starts on the coordinate anomaly,
crosses the horizon, and hits the singularity.
This photon reaches the singularity within an extremely short time interval after emission.
The photon labeled (e) starts between the coordinate anomaly
$R_M$ and the shrinking light-like boundary $R_B$.
The photon crosses the coordinate anomaly $R_M$ and the black hole event horizon
$R_H$  (the outgoing causal horizon),
and hits the singularity.
The photon labeled (f) starts outside the shrinking
light-like boundary $R_B$, and never crosses the black hole horizon.
These trajectories are plotted on the Penrose diagram in Fig. \ref{InPen}.
\setfigure{InPen}{Penrose diagram with ingoing photon trajectories over plotted.}

\setcounter{equation}{0}
\section{Conclusions}
\indent

We have calculated the paths of outgoing and ingoing photons
originating near an evaporating black hole.  Several behaviors
are found to be of interest.  It is found that outgoing photons
just outside the horizon initially approach the singularity, until the
mass of the black hole decreases to a point that allows the
photons to alter direction and begin to move away from the
singularity.

When plotted on the Penrose diagram for this space-time,
the photon trajectories confirm the expected conformal
coordinate representation.  Consistent with expectations, outgoing photons are always
stationary at $r=R_M$; it is the shrinking mass scale $R_M$
that moves away from those photons.  In addition, because the black
hole eventually evaporates away, there is a class of
photons that can never reach the singularity.

Analogous results have been obtained for an accreting black
hole, which will be presented in a subsequent paper.

\begin{center}
\textbf{Acknowledgments}
\end{center}
BAB would like to acknowledge the support of the NASA Administrator's
Fellowship Program.


\section{Introduction}

The interest in in-medium properties of hadrons has been growing over the
last decade because of a possible connection with broken symmetries of QCD
and their partial restoration inside nuclear matter. Nearly two decades
ago there were sophisticated theoretical predictions that masses of vector
mesons should generally decrease in medium as a function of density due to
a partial restoration of chiral symmetry \cite{Brown-Rho,HL}.
Specifically, there existed well-worked out predictions that the scalar
$2\pi$ strength should decrease in medium and that the vector meson masses
should drop \cite{Kuni,Klingl}. All of these calculations were performed
for idealized situations (infinite cold nuclear matter at rest) and little
attention was paid to the actual observability of these predicted changes.
At the same time experiments (CERES, TAPS) seemed to show the predicted
behavior. The CERES results indicated a significant broadening of the
$\rho$ meson in medium \cite{Ceres}, whereas the TAPS results on the
$2\pi$ strength exhibited the predicted lowering of the $\sigma$ strength
inside nuclei \cite{TAPS2pi}. Most recently, the CB/TAPS@ELSA experiment
has also obtained an indication for a lowering of the $\omega$ meson mass
in nuclei \cite{Trnka}. Other interesting data in this context have been
obtained by groups at JLAB \cite{g7}, KEK \cite{KEK} and CERN \cite{NA60}.
While experiments with relativistic or ultrarelativistic heavy ions have
the advantage that high densities can be produced, they have the
disadvantage that any observed signal naturally is a time-integral over
the whole collision history. On the other hand, experiments with
microscopic probes on cold nuclei have the advantage that they proceed
closer to equilibrium, even though the densities probed are much lower.

We have focussed our work on in-medium properties of hadrons on two steps.
First, we have performed state-of-the-art calculations of vector meson
spectral functions in cold nuclear matter. Second, we have followed
closely the CB/TAPS@ELSA experiment and have performed various feasibility
studies and analyses of this experiment searching for in-medium changes of
the $\omega$ meson in medium. In a third step we have also analyzed the
lowering of scalar strength in nuclei, predicted as a signature for chiral
symmetry restoration, and indeed observed by the TAPS experiment.

\section{In-Medium Properties of Vector Mesons}

\subsection{QCD sum rules and hadronic models}

While the predicted drop of the chiral condensate with density
\cite{Weise} is a rather robust result  the connection between chiral
condensates on one hand and hadronic spectral functions on the other
cannot simply be inferred from looking at the dependence of the condensate
on density and temperature. Instead, the connection is much more indirect;
only an integral over the spectral function can be linked via the QCD sum
rule to the condensate behavior in medium. In an abbreviated form this
connection is given by
\begin{eqnarray}
R^{{\rm OPE}}(Q^2) & = & {\tilde c_1 \over Q^2} + \tilde c_2 -{Q^2 \over
\pi} \int\limits_0^\infty \!\! ds \, {{\rm Im}R^{{\rm HAD}}(s) \over
(s+Q^2)s}
  \label{eq:opehadr}
\end{eqnarray}
with $Q^2:= -q^2 \gg 0$ and some subtraction constants $\tilde c_i$. Here
$R^{\rm OPE}$ represents a Wilson's operator expansion of the
current-current correlator in terms of quark and gluon degrees of freedom
in the space-like region. It is an expansion in terms of powers of $1/Q^2$
and contains the condensates as expansion parameters. On the other hand,
$R^{{\rm HAD}}(s)$ in (\ref{eq:opehadr}) is the same object for time-like
momenta, represented by a parametrization in terms of hadronic variables.
The dispersion integral connects time- and space-like momenta. Eq.
(\ref{eq:opehadr}) connects the hadronic with the quark world. It allows
to determine parameters in a hadronic parametrization of $R^{{\rm
HAD}}(s)$ by comparing the lhs of this equation with its rhs. Invoking
vector meson dominance it is easy to see that for vector mesons Im$R^{{\rm
HAD}}(s)$ in (\ref{eq:opehadr}) is just the spectral function of the meson
under study.

The operator product expansion of $R^{\rm OPE}$ on the lhs involves quark-
and gluon condensates \cite{Leupold,LeupoldMosel}; of these only the
two-quark condensates are reasonably well known, whereas our knowledge
about already the four-quark condensates is rather sketchy.

Using the measured, known vacuum spectral function for $R^{\rm HAD}$
allows one to obtain information about these condensates appearing on the
lhs of (\ref{eq:opehadr}). On the other hand, modelling the
density-dependence of the condensates yields information on the change of
the hadronic spectral function when the hadron is embedded in nuclear
matter. Since the spectral function appears under an integral the
information obtained is not direct. Therefore, as Leupold et al.\ have
shown \cite{Leupold,LeupoldMosel}, the QCDSR provides important
constraints for the hadronic spectral functions in medium, but it does not
fix them. Kaempfer et al have turned this argument around by pointing out
that measuring an in-medium spectral function of the $\omega$ meson could
help to determine the density dependence of the higher-order condensates
\cite{Kaempfer}. Leupold has also argued that the four-quark condensates
are interesting, because a special one of them can be unambiguously linked
to a difference of spectral functions in the axial and vector channel
\cite{Leupold4q}.

Thus models are needed for the hadronic interactions. The quantitatively
reliable ones can at present be based only on 'classical' hadrons and
their interactions. In lowest order in the density the mass and width of
an interacting hadron in nuclear matter at zero temperature and vector
density $\rho_v$ are given by (for a meson, for example)
\begin{eqnarray}     \label{trho}
{m^*}^2 & =  & m^2 - 4 \pi \mathcal{R}e f_{m N}(q_0,\theta = 0)\, \rho_v
\nonumber \\
m^* \Gamma^* & = & m \Gamma^0 -  4 \pi \mathcal{I}m f_{mN}(q_0,\theta =
0)\, \rho_v ~.
\end{eqnarray}
Here $f_{mN}(q_0,\theta = 0)$ is the forward scattering amplitude for a
meson with energy $q_0$ on a nucleon. The width $\Gamma^0$ denotes the
free decay width of the particle. For the imaginary part this is nothing
other than the classical relation $\Gamma^* - \Gamma^0 = v \sigma \rho_v$
for the collision width, where $\sigma$ is the total cross section. This
can easily be seen by using the optical theorem. The calculation of
hadronic spectral functions in medium thus reduces -- in lowest order in
the density -- to a calculation of hadron-nucleon scattering amplitudes.

Unfortunately it is not a-priori known up to which densities the
low-density expansion (\ref{trho}) is useful. Post et al. \cite{Post1}
have recently investigated this question in a coupled-channel calculation
of selfenergies. Their analysis comprises pions, $\eta$-mesons and
$\rho$-mesons as well as all baryon resonances with a sizeable coupling to
any of these mesons. The authors of \cite{Post1} find that already for
densities less than $0.5 \rho_0$ the linear scaling of the selfenergies
inherent in (\ref{trho}) is badly violated for the $\rho$ and the $\pi$
mesons, whereas it is a reasonable approximation for the $\eta$ meson.
This may serve as a warning sign for many in-medium calculations that use
the low-density approximation without any further proof of its
reliability. Reasons for this deviation from linearity are Fermi-motion,
Pauli-blocking, selfconsistency and short-range correlations. For
different mesons different sources of the discrepancy prevail: for the
$\rho$ and $\eta$ mesons the selfconsistency acts against the low-density
theorem by inducing a strong broadening for the $D_{13}(1520)$ and a
slightly repulsive mass shift for the $S_{11}(1535)$ nucleon resonances to
which the $\rho$ and the $\eta$ meson, respectively, couple. The
investigation of in-medium properties of mesons, for example, thus
involves at the same time the study of in-medium properties of nucleon
resonances and is therefore a coupled-channel problem.

\begin{figure}[htb]
\centerline{\epsfig{file=rhospect.eps,height=14cm}} \caption{Spectral
function of the $\rho$ meson in nuclear matter at density $\rho_0$ for
various momenta indicated in the figure. The left column shows the
transverse spectral function, the right column that of longitudinally
polarized $\rho$ mesons. The thin dotted line in each figure is the vacuum
spectral function, the other curves give the effect of selfconsistency and
short-range correlations (from \protect\cite{Post1}).} \label{rhospect}
\end{figure}
The strong interplay between the $\rho$ meson and the D13(1520)-nucleon
hole excitation leads to a dominant lower hump in the $\rho$ spectral
function also in this relativistic and selfconsistent calculation; it
confirms our earlier result obtained in a more simplified approach.
Whereas the longitudinal component of the $\rho$ meson only broadens
somewhat, the transverse component shows a major distortion which evolves
as a function of the $\rho$ momentum (see Fig. \ref{rhospect}) and
reflects the coupling to the D13(1520) resonance. It is obvious that such
a rich structure cannot be obtained from QCD sumrules that involve only an
integral over the spectral function. At the same time, the D13(1520)
resonance broadens considerably due to the opening of phase-space for
$\rho$-decay. For the $\eta$ meson the optical potential resulting from
our model is rather attractive whereas the in-medium modifications of the
S11(1535) are found to be quite small.

\subsection{Coupled channel calculation of hadronic selfenergies}
Besides the range of validity of the low-density approximation the
meson-nucleon cross sections are often not directly measurable and can be
inferred only from unitarity arguments involving many coupled channels.
Unitarity thus provides important contraints for reactions such as $\omega
N \to 2\pi N$ that are experimentally not accessible, but are not only
needed for but contribute significantly to a calculation of the $\omega$
selfenergy. Bearing this in mind we have recently performed a calculation
of the selfenergy of the $\omega$ meson in medium.
\begin{figure}[htb]
\centerline{\epsfig{file=K-omspect.eps,height=6cm}}
 \caption{Spectral function of the $\omega$ meson in nuclear matter at rest,
at densities $0$, $\rho_0$ and 2 $\rho_0$ (from
\protect\cite{Muehlich:2006nn}).} \label{omspectr}
\end{figure}
This calculation is based on a unitary coupled channel analysis of all
existing $\pi N$ and $\gamma N$ data up to an invariant mass of 2 GeV
\cite{Muehlich:2006nn}. The coupled channel character of this calculation
is of utmost importance here because it is the only way to include
experimental constraints on the $2\pi$ decay channel that was found to be
dominant in \cite{Klingl}. This analysis and thus also the selfenergy of
the $\omega$ meson extracted from the $\omega N$ scattering amplitude
gives a broadening of about 60 MeV at $\rho_0$ for an $\omega$ at rest and
a small upward shift of the peak mass. In addition, due to a nonzero
coupling of the $\omega$ to the S11(1535) resonance the $\omega$ spectral
function exhibits a small peak at a mass of around 550 MeV. This
calculation gives, for the first time, the $\omega$ selfenergy also for
nonzero momenta (which corresponds to the experimental situation). The
result of this calculation is shown in Fig. \ref{omspectr}.

For vanishing $\omega$ momentum this result qualitatively agrees with that
of \cite{Lutz} in that it yields a lower mass structure in the spectral
function at an invariant mass of 500 - 600 MeV and only a very small shift
of the main peak; the latter is in contrast to the results of
\cite{Klingl} (for a recent discussion of the results obtained in
\cite{Klingl} see \cite{Eichst}) and of \cite{Thomas}. The latter was
based on a relativistic mean-field model and does not contain any
dispersive effects. Fig.\ \ref{omspectr} also shows clearly that the
$\omega$ meson in medium acquires a significant broadening, caused by
collisions with the surrounding nucleons.

\subsection{Data analysis with GiBUU}
Our second line of approach to the problem of in-medium selfenergies has
concentrated on an analysis of the recent CB/TAPS@ELSA data \cite{Trnka}
Since the experiment looks for the channel $\gamma + A \to A^* + \omega
\to A^* + \gamma + \pi^0$ it is mandatory to control the effects of final
state interactions on the $\pi^0$ in a quantitative way. The only method
available for this is that of coupled channel semiclassical transport
calculations which -- as we had shown earlier in extensive work -- can
give a consistent description of many experimental phenomena, both in
heavy-ion \cite{Lari} as well as in nucleon-, pion- \cite{Buss}, photon-
\cite{Buss-Leitner} and neutrino-induced reactions \cite{Leitner}. For any
reaction on nuclei with hadrons in the final state a state-of-the-art
transport calculation of the final state interactions is an indispensible
part of the theory. We have, therefore, spent significant effort on
developing a new code, dubbed 'GiBUU', for the transport calculations.
This code is written in object-oriented FORTRAN 95/2003 and incorporates
all the experience we have gained with earlier numerical implementations
at Giessen over the last 20 years \cite{GiBUU}. A major effort has gone
into an analysis of the $\omega$ photo-production experiment at
CB/TAPS@ELSA.

For $\omega$ mesons the CB/TAPS@ELSA collaboration has measured the
nuclear transparency \cite{Metag,Kot}, see Fig.\ \ref{TA}. This
transparency gives directly the imaginary part of the meson's selfenergy
in medium; using a low-density approximation one can then extract the
inelastic $\omega N$ cross section. The figure shows an analysis of these
data for various values of the collisional width, assumed at a momentum of
1.1 GeV, the average momentum of the $\omega$ mesons in this experiment.
Again, care has been taken to ensure consistency between the collisional
width and the inelastic $\omega N$ cross sections. It is seen from this
figure that a width of about 130 - 150 MeV describes the data. This width
is larger by about a factor of 3 than that previously assumed. A fit to
the data can actually also determine the momentum-dependence of this cross
section \cite{Kot}.
\begin{figure}[htb]
\centerline{\epsfig{file=TA_omega.eps,height=8cm}} \caption{Experimentally
determined transparency ratio (squares) as a function of massnumber in
comparison with GiBUU calculations assuming various values for the
inleastic width at an $\omega$ momentum of 1.1 GeV (from
\protect\cite{Metag,Kot}).} \label{TA}
\end{figure}
With this method we have also analysed results on the experimental
determination of the nuclear transparency ratio for $\phi$ mesons
\cite{Muhlich:phi}, measured by a group at SPRING8. Also in this case an
unexpectedly large inelastic cross section for $\phi N$ interactions was
extracted. We have found that indeed cross sections about a factor 3
larger than those theoretically expected are needed to explain the
mentioned data, in line with a simple Glauber analysis by the SPRING8
group.

Our simulations give a full event analysis and thus allow to calculate
also background contributions on the same footing as the actual signal.
They also allow insight into the effects of rescattering of the pions
produced in the decay of the $\omega$ meson and have suggested a method to
suppress the rescattered pion background that has actually been adopted by
the experimental group. The result of this analysis is that the $\omega$
spectral function published by the CB/TAPS@ELSA collaboration can be
explained only if a lowering of the $\omega$ meson mass in medium by about
16 \% is assumed together with the appropriate collisional broadening.

A problem in this context is that the experiment does not determine the
spectral function of the $\omega$ meson itself. Instead, we have noted
that the result of the experimental analysis is the product of the primary
production cross section with the spectral function and the partial decay
probability into the channel under study ($\pi^0 \gamma$ here). If the
first and the latter depend strongly on the invariant mass itself, as it
does for the CB/TAPS@ELSA experiment, then significant distortions of the
spectral function may arise. This is a topic that is presently under
intensive study \cite{MosEr}.

\section{In-medium Properties of Scalar Mesons}

Since the restoration of chiral symmetry should lead to a convergence of
scalar-isoscalar and scalar-isovector strength at high densities it was
predicted that the decay of the scalar $\sigma$ meson into two uncharged
pions should become increasingly difficult and the scalar strength should
become narrower at the $2\pi$ threshold \cite{Kuni}. The TAPS
collaboration had indeed initially seen this effect in the $2\pi^0$ data
as predicted, whereas a comparison measurement in the $\pi^0 \pi^{\pm}$
channel did not show such an effect \cite{TAPS2pi}. Various explanations
for these findings have been advanced by the Valencia group \cite{Oset2pi}
and by a group in Lyon \cite{Chanfray} in terms of $\pi - \pi$
correlations or chiral symmetry restoration in nuclei.

\begin{figure}[htb]
\centerline{\epsfig{file=2piCa.eps,height=8cm}} \caption{Data of the TAPS
collaboration (Bloch et al.) for 2$\pi^0$ photoproduction on a $^{40}Ca$
target for two different photon energies. The solid curve gives the result
of a GiBUU calculation \protect\cite{Buss1}, the dashed curves in the
semicharged 2$\pi$ channel are normalized to the data (see text). Data
from \protect\cite{Bloch}.} \label{2pi}
\end{figure}

None of these calculations, however, did look into the simplest possible
explanation of the observed effects in terms of mundane pion rescattering.
We have, therefore, performed such calculations \cite{Buss} using the
GiBUU method which is ideally suited for this task. These calculations,
which did not contain any effects connected with $\pi \pi$ correlations,
did reproduce the observed effect for the $\pi^0 \pi^0$ channel, but they
also predicted a similar effect in the semi-charged channel where it had
not been seen experimentally. However, a more recent analysis with higher
statistics by Bloch et al. \cite{Bloch}  yielded a result for both charge
channels that is in perfect agreement with our calculations (see Fig.
\ref{2pi}; the dashed lines in this figure are normalized in height to the
data, this normalization reflects uncertainties in the elementary cross
sections). In particular the yield in the semi-charged channel is strongly
influenced by a coupled-channel effect, the charge transfer in $\pi N$
interactions; Glauber based absorption models miss this contribution. This
illustrates that a very sophisticated treatment of final state
interactions is absolutely mandatory when looking for more 'exotic'
effects in nuclei. We conclude that any analysis of the $2\pi^0$ data with
respect to a lowering of the scalar strength in nuclei has to take the
pion rescattering effects into account. Present day's data are all
consistent with simple rescattering.

\section{Summary}
The main message we have learned from the studies reported here is that it
is important to calculate not only in-medium properties under idealized
conditions (static, uniform matter in equilibrium), but to also explore
the influence of these properties on actual observables. The spectral
function itself, which contains the information on in-medium selfenergies,
in particular in-medium masses and widths, is not directly observable.
Instead, both the creation of the studied hadron as well as its decay
influence the observables as much as the spectral function itself and thus
have to be under good control. The same holds for the final state
interactions on hadronic decay products. Here a state-of-the-art treatment
of final state interactions is mandatory. There is now general agreement
on the amount of collisional broadening of vector mesons in medium, but
the verification of an actual mass-shift still requires more work, both
theoretical and experimental. Until then the quest for chiral symmetry
restoration in nuclei remains a challenging topic.

\section{Acknowledgement}
This talk is mainly based on results obtained by Oliver Buss and Pascal
Muehlich as part of their doctoral theses. I thank them, and the entire
GiBUU group, for many inspiring discussions.