%\documentclass{aastex}

%% preprint produces a one-column, single-spaced document:

%\documentclass[12pt,preprint]{aastex}
\documentclass[preprint]{aastex}

%% preprint2 produces a double-column, single-spaced document:

% \documentclass[preprint2]{aastex}

%%%%%%\documentstyle[psfig,graphics,amssymb]{elsart}


\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\lapprox}{\lower 3pt\hbox{$\buildrel < \over \sim\;$}}
\newcommand{\gapprox}{\lower 3pt\hbox{$\buildrel > \over \sim\;$}}
\newcommand{\alfatil}{\~{\alpha}}
\def\etal{{\it et al.}}
\def\lae{\mathrel{<\kern-1.0em\lower0.9ex\hbox{$\sim$}}}
\def\gae{\mathrel{>\kern-1.0em\lower0.9ex\hbox{$\sim$}}}
\def\deg{^{\circ}}
\def\sec{^{\prime\prime}}
\def\kms{km~s$^{-1}$}
\def\mone{$^{-1}$}
\def\mtwo{$^{-2}$}
\def\mthree{$^{-3}$}
\newcommand{\Msun}{$M_{\sun}$}
\newcommand{\Zsun}{$Z_{\sun}$}


\shorttitle{Evolution of Global Properties of Powerful Radio Sources I}
%%\shortauthors{Carvalho et al.}


\begin{document}

\title{Evolution of Global Properties of Powerful Radio Sources I -
Hydrodynamical Simulations in a Constant Density Atmosphere and
Comparison with
Self-similar Models}

\author{Joel C. Carvalho\altaffilmark{1} and Christopher P. O'Dea}
\affil{Space Telescope Science Institute, 3700 San Martin Dr.,
Baltimore, MD 21218}

\altaffiltext{1}{Departamento de Fisica, UFRN, C.P. 1661, CEP 59072-970,
Natal,
RN, Brazil}
%\author{Draft: \today }

\begin{abstract}

We have carried out 2-D axisymmetric numerical simulations of light,
supersonic jets
propagating into constant density atmospheres over a broad range of
parameter space. We
examine the evolution of the global properties of the sources as a
function of source
size, for a range of Mach number, density contrast,  and jet power. We
also compare our
results with the expectations of current analytical self-similar models.
The derived
global parameters are not sensitive to small changes in input
parameters. The material
content of the relativistic jet does effect its propagation and the
overall structure and
energetics (see Paper II) of the source. The jet in Case 1 (pressure
dominated by
relativistic electrons) propagates faster, the cocoon has an smaller
aspect ratio, the
cocoon maintains higher pressure, and the bow shock expands at higher
lateral velocity
than Case 2 jets (pressure dominated by relativistic electrons and
protons). However,
these differences generally become less important as the jet Lorentz
factor increases. We
have confirmed the suggestion of Wilson \& Falle (1985) that the
location of the first
shock in the jet depends linearly on the jet radius and Mach number. We
show that the
result is largely independent of density contrast and that over a wide
range in Mach
number the relation can be used to make a crude estimate of the jet Mach
number. We find
that the pressure in the cocoon and the bow-shocked region varies both
along and across
the source axis. The degree to which the cocoon comes into pressure
balance with the
ambient medium depends on the jet Mach number, with intermediate $M$
jets having cocoons
which are in pressure balance with the ambient medium over much of the
length of the
source; while high and low $M$ jets have cocoons which are over
pressured or under
pressured, respectively. The bow shock lateral expansion decelerates
quickly behind the
jet head, but the velocity remains at least somewhat supersonic. On the
other hand, the
cocoon lateral expansion velocity quickly drops behind the jet head and
becomes subsonic
with respect to the post bow shock gas. The strong time dependence of
the cocoon and bow
shock pressure and lateral expansion speed is expected to influence the
properties of the
emission line nebula associated with the radio source. The slow
expansion of the cocoon
and bow shock over most of their length will limit the maximum velocity
to which ambient
clouds can be accelerated. The behavior of the time evolution of the
source size ($z_h
\propto t^m$) is in general more complex than the predictions of the
self-similar models.
The time dependence of the source size exhibits up to  3 phases whose
relative length
depends on the jet Mach number. The radius of the jet head is a weak
function of density
ratio and a stronger function of Mach number, with higher Mach number
jets tending to have
smaller jet heads. The advance speed of the jet head depends on the head
radius in
agreement with the expectations of ram pressure balance. Due to the
dependence of the
growth rate of the jet head with Mach number, the advance speed of the
head decreases
slowly for low Mach number and remains roughly constant for high Mach
number. The volume
of the region enclosed by the bow shock increases with time at a rate
which is close, but
not identical to the predictions of the constant velocity self-similar
models. The bow
shock expansion rate seems to increase with time for small Mach number
but is relatively
stable at high Mach number. The cocoon volume generally increases more
slowly  than the
volume of the region enclosed by the bow shock. In low Mach number jets,
the dissipation
of energy through internal shocks is able to disrupt the jets and turn
off the supply of
energy to the cocoon. At this stage, the bow shock propagates as a sound
wave and the
morphology of the inferred radio emitting material becomes faint and
diffuse. Thus, we
suggest that low Mach number jets may be able to transition between FR
II and FR I type
radio morphologies.
\end{abstract}

\keywords{galaxies: jets -- galaxies: intergalactic medium --
hydrodynamics -- radio continuum: galaxies}

%******************************************************************************

%******************************************************************************

\section{INTRODUCTION}

Powerful radio galaxies play a critical role in our understanding of
both galaxy evolution
and the phenomenon of activity in galactic nuclei. They provide probes
of the intracluster
medium, particularly in cooling  flow clusters where they reveal the
existence of large
scale ordered magnetic  fields in the ICM (Taylor \etal\ 1994). They
interact with gas in
the ISM of their  host galaxies producing the alignment effect and
therefore can be used
as probes of the content and structure of the ISM (McCarthy 1993). The
higher redshift
sources act as beacons which reveal the location of the earliest sites
of the formation of
giant ellipticals in the universe (e.g., Spinrad \etal\ 1995; Stern \&
Spinrad 1999).
Powerful radio galaxies act
as probes of their environments as a function as epoch and may allow
cosmological
parameters to be estimated (Daly 1994; Guerra, Daly \& Wan 2000). Their
luminosities,
morphologies, and lifetimes constrain the processes through which
angular  momentum is
extracted from gas in the central accretion disk around a super-massive
Black Hole and
jets are formed (e.g., Blandford 1994).  Yet we know  little about how
the radio galaxies
are born and how they subsequently evolve.

The standard view of the evolution of strong double radio sources is
that a pair of
supersonic jets are ejected from the central object in the parent galaxy
and tunnel
through the ambient gas. The interaction of the jet with the medium
develops an expanding
double shock that forms the lobe structure observed in the radio maps.
There is much
recent interest in understanding how the radio luminosity and radio
source size vary over
the lifetime of the radio source, and especially in constraining this
evolution through
use of the power-diameter (P-D) diagram (Baldwin 1982; Fanti et al.
1995; Readhead et al.  1996; Kaiser, Dennett-Thorpe \& Alexander 1997; 
O'Dea \& Baum 1997; Blundell, Rawlings \& Willott 1999; De Vries \etal\
2000; Alexander 2000). 
The main current limitations are the
simplicity of the
analytical models which are used, and selection effects in the current
radio source samples.

In the majority of the models for radio source expansion, it is assumed
that the cavity
opened in the interstellar medium by the jet has a uniform pressure
(e.g., Scheuer 1974;
Carvalho 1985; Begelman \& Cioffi 1989; Begelman 1996). However,  we
observe that the
lobes of powerful FR II radio sources are much brighter near the end of
the jets (i.e.,
where the jet is depositing its momentum) than near the parent galaxy.
Thus, it is clear
that the lobe properties are non-uniform and the implications of this
for the radio source
propagation and evolution need to be understood. Most existing
analytical models assume
self-similarity. Thus, we also address in detail the issue of to what
extent and under
what circumstances the jet propagation and lobe expansion are
self-similar.

In this paper (the first of the series) we report on  results of set  of
numerical
simulations of hydrodynamic supersonic jets covering a very broad range
in parameter space
which we have carried out in order to improve our understanding of  the
global evolution
of radio sources. Our goal is to describe the evolution of global
properties of the radio
source as a function of source size over a broad range of Mach number
and density
contrast. Here we focus on jets propagating into a constant density
atmosphere. In
Carvalho \& O'Dea (2001a, hereafter Paper II) we expand our analysis to
include jets
propagating into atmospheres of declining density and pressure and
discuss the energetics
of sources in both constant and declining density atmospheres. In
Carvalho \& O'Dea
(2001b, hereafter Paper III)  we will apply the knowledge gained from
these simulations to
a more detailed model for radio galaxy evolution.


In Section $\S$~\ref{sec-sim}, we discuss the relation between
relativistic and
non-relativistic simulations and calculate the parameters of the
classical equivalent jets
that will be simulated and describe the numerical technique. In Section
$\S$~\ref{sec-num}
we discuss differences in the results of 3D and 2D simulations with
special emphasis on intrinsic limitations to the results that can be
derived
from axisymmetric 2-dimensional simulations. And we present the results
of a series of
numerical hydrodynamical simulations performed with
the purpose of studying the main morphological and dynamical properties
of the radio
sources. We describe  the structure and dynamics of the jet, cocoon, and
bow shock.  In
Section $\S$~\ref{sec-similar} we review the three basic classes of
self-similar
analytical models with which we will compare our numerical results. In
Section
$\S$~\ref{sec-global} we compare our results on the jet head propagation
and the evolution
of the source volume with the self-similar models.  In Section
$\S$~\ref{sec-jetpower} we
discuss the behavior of very low power jets, the generation of sound
waves by the radio
source, and a possible transition from FR II to FR I type morphology. In
Section
$\S$~\ref{case1}, we discuss the dependence of our results on the
relativistic particle
content of the jets. In Section $\S$~\ref{RunC}, we discuss the
robustness of our results
against small variations in input parameters. In Section
$\S$~\ref{sec-summary} we
summarize our conclusions.


%******************************************************************************

%******************************************************************************

\section{SIMULATIONS OF ASTROPHYSICAL JETS}
\label{sec-sim}

During the past two decades, numerical simulations have arguably
increased
our insight into the behavior and physics of astrophysical jets
including
the structure, dynamics and stability of both non-relativistic (Rayburn
1977;
Norman et al. 1982; Lind et al.  1989;
Clarke, Norman \& Burns 1989; Hardee \& Norman 1990; K\"ossl, M\"uller
\&
Hillebrandt 1990; Cox, Gull \& Scheuer 1991; Falle 1991; Clarke \& Burns

1991; Cioffi \& Blondin 1992; Loken et al. 1992; Hardee et al. 1992;
Ferrari, Massaglia \& Bodo 1996) and relativistic jets (Wilson 1987;
Bowman, Leahy \& Komissarov 1996; Mioduszewski, Hughes \& Duncan
1997; Marti et al. 1997; Komissarov \& Falle 1997, 1998; Rosen et al.
1999),
with and without magnetic field.
In addition, some authors have analyzed the synchrotron emission and
produced radio maps based on the numerical simulations, which enable a
more realistic study of the morphology and of effects such as Doppler
beaming
(e.g., Mioduszewski, Hughes \& Duncan 1997).

%******************************************************************************

%******************************************************************************

\subsection{Relativistic versus Non-relativistic Jet Simulation}

There are several pieces of evidence for relativistic bulk motion in the
jets of
powerful
radio sources including (1) apparent superluminal motion (Zensus 1997),
(2) strong
side-to-side
asymmetry (Bridle \& Perley 1984), and (3) the Laing-Garrington effect
(Laing 1988;
Garrington \etal\ 1988).
However, most of the  hydrodynamic simulations performed to study the
jet
phenomenon,
including the ones in this work, are non-relativistic. This raises the
issue of how
to
compare the non-relativistic simulations with relativistic astrophysical
jets.

In the last few years a number of numerical studies of relativistic jets
have been
carried
out (Duncan \& Hughes 1994; Gomez et al 1995; Marti et al 1997;
Mioduszewski, Hughes
\&
Duncan 1997; Komissarov \& Falle 1997; Aloy \etal\ 1999). Initial
comparisons of the
dynamics and morphology of relativistic and non-relativistic jets
suggested that
relativistic jets develop a much smaller cocoon, exhibit a reduced
Kelvin-Helmholtz
type
instability and propagate far more efficiently than their Newtonian
counterparts
(Duncan
\& Hughes 1994; Marti, M\"uller \& Ibanez 1994; Marti et al. 1997).
However, it was
later
pointed out by Komissarov \& Falle (1998) that it is necessary to
properly define
the
equivalent non-relativistic and relativistic jets in order to obtain a
meaningful
comparison.
 This is achieved by using the appropriate expressions for the
density ratio and Mach number. For example, the two types of jet should
have the
same
thrust. In fact,  non-relativistic and relativistic jets that are
equivalent in
thrust are
quite similar in their gross properties. However, because the equivalent
non-relativistic
jet is heavier, it does not exhibit as well developed a cocoon. A
detailed
comparison has
been carried out by Rosen et al. (1999) who show that the overall
dynamics of both
relativistic and non-relativistic jets are the same, although there are
major
differences
of the jet internal velocity field which makes it infeasible to use
non-relativistic
simulations to calculate the expected radio intensity distribution.

%******************************************************************************

%******************************************************************************

\subsection{Non-relativistic Equivalents of Relativistic Jets}
\label{nonrel}

Rosen et al. (1999) compared both equal power and equal thrust
equivalent jets.
Here, for
simplicity, we shall consider only the latter which seems also to be the
more
interesting one since it should give the same jet head to jet velocity
ratio. Rosen
et al.
give the expression for the jet density $\rho_c$ and Mach number $M_c$
of the
classical
simulation in terms of the relativistic jet parameters. They are
\beq
\rho_c = \gamma_j^2 \rho_j \left(1+{\Gamma\over \Gamma-1}{P_j\over
\rho_j
c^2}\right),
        \label{S1}
\eeq
\beq
M_c^2 = (\gamma_j^2-1) {1\over \Gamma-1} \left(1+{\Gamma -1\over\Gamma
}{\rho_j
c^2\over
P_j}\right).
                                \label{S2}
\eeq
where $\gamma_j$ is the jet Lorentz factor, $\Gamma$ is the adiabatic
index and
$\rho_j$
and $P_j$ are its density and pressure, respectively.

Because the astrophysical jet is emitting synchrotron emission, it must
have a
population of relativistic electrons (and/or positrons). As we shall
discuss
in Paper III, the relativistic particle content can be either a pair
plasma,
i.e.,  electrons (and positrons) only (Case 1),  or baryon dominated -
electrons plus protons (Case 2). In this work we shall assume that the
relativistic
component dominates over the thermal component and the pressure is then
given by
\beq
P_j = {1\over 3} (1+k) u_e                         \label{S3}
\eeq
where $ k \simeq 0$ in Case 1 and $k\neq 0$, in Case 2, $k$ is the ratio
between the
proton and electron energy. The energy density of the electrons is
\beq
u_e = \bar{\gamma_e} n_e m_e c^2                  \label{S4}
\eeq
where $n_e$ is the particle density and $\bar{\gamma_e}$ is the average
Lorentz
factor of
the electrons.


We can now calculate the non-relativistic jet parameters $\eta_c$ and
$M_c$ from
(\ref{S1}) - (\ref{S4}) and obtain
\beq
\eta_c = \gamma_j^2 \eta (1+\delta_2)       \label{S5}
\eeq
\beq
M_c = \sqrt{3} \gamma_j \beta_j (1+\delta_2^{-1})^{1/2}   \label{S6}
\eeq
where $\eta = n_j/n_a$, $n_a$ is the ambient medium particle density and
(see Paper
III)
\beq
\delta_2 = {4\over 3} (1+k) \bar{\gamma_e} {m_e\over m_p}q.   \label{S7}
\eeq
The parameter $q=n_e/n_j$ is defined in Paper III and we have adopted
$\Gamma=4/3$
and
assumed that for the ambient medium $P_a/\rho_a c^2 \ll 1$.

In Table 1 we list the values of $\eta_c$ and $M_c$ (Runs A1-A10) for
values of the
jet
Lorentz factor in the range 2 - 10, taking a characteristic density
ratio\footnote{This is almost two orders of magnitude higher than the
value found by, e.g.,  Alexander \&
Pooley (1996) for the jet in Cygnus A ($<4\times 10^{-5}$). However,
such small
values of
$\eta$ make the numerical simulations extremely lengthy.} $n_j/n_a \sim
10^{-3}$.

The values of $\delta_2$ were calculated for $k\simeq 0$ (Case 1) and
$k=100$ (Case 2) neglecting the thermal component, that is, adopting
$q=1/2$. The average Lorentz factor $\bar{\gamma_e}$ was calculated
assuming a power law spectrum for the energy distribution with
energy index 2.5 and with lower and upper cutoff frequencies equal to
$10^7$ and $10^{10}$ Hz respectively. Note that using an upper cutoff
frequency of $10^{11}$ Hz  would change our values of Mach number and
density ratio by typically a few percent and would not significantly
change
our sampling of parameter space or our results (cf. \S~\ref{RunC}).


We see that the main difference between Case 1 and 2 is that, for a
given
$\eta$ and Lorentz factor $\gamma_j$, the equivalent classical jets in
Case 2 are relatively
heavier than in Case 1 and have lower Mach number. This can be
understood
as being due to the presence of the proton component which increases the
total enthalpy since their energy is much larger than that of the
electrons
($k=100$). If we had chosen the thermal component as the dominant
component,
the equivalent density contrast would be simply given by $\eta_c =
\gamma_j^2 \eta$ and the Mach number $M_c={\gamma_j \beta_j \over
\gamma_s
\beta_s}$ with the subscript $s$ referring to the sound speed
(Komissarov \& Falle
1998).  The equivalent classical jet speed would have to be extremely
large even for moderate $\gamma$-factors. The differences between Case 1 and 2
are discussed further in \S~\ref{case1}.

We continued to explore parameter space in a second set of simulations
(Runs B1-B10) which were carried out by first fixing the density
($\eta_c=0.1$) and varying $M_c$ and then fixing the Mach number at
6 and varying $\eta_c$. This allows us to more easily appreciate
how each parameter affects the properties of the resulting simulations.
It also includes a run widely studied in the literature, that is,
$M=6$ and $\eta=0.1$. Finally, in the third run (C1-C6) we vary only
slightly one of the parameters in Runs B3 and B8 while keeping the
other fixed in order to estimate how sensitive the results of the
simulations
are to slight changes in $\eta$ and $M$ (see \S~\ref{RunC}).

For our non-relativistic simulations, we have chosen the values of
$\eta_c$ and
$M_c$
shown in Table 1. These values of $\eta_c$ and $M_c$ are shown in Fig.
\ref{ASU_EtaMach}
where the contours of equal $k$ and equivalent Lorentz factor $\gamma_j$
are drawn.
Run A
is represented by solid circles, Run B by open circles and Run C by
crosses. We also
see
that Run B allows us to more thoroughly cover the $\eta - M$ plane.


%******************************************************************************

%******************************************************************************

\subsection{Numerical Setup}

We performed numerical two-dimensional axisymmetric hydrodynamical
simulations using
the
MHD code ZEUS-3D. This is an Eulerian code developed at the National
Center for
Supercomputing Applications (NCSA). The ZEUS algorithm solves the ideal
MHD
equations
using a time-explicit, multi-step, finite difference scheme (Stone \&
Norman
1992a,b;
Stone, Mihalas \& Norman, 1992). It can be used in the study of a great
variety of
astrophysical systems and has been widely tested.

Since our final objective is to derive the behavior of {\it global\/}
properties of
the jets and
their scaling with source size as a function
of Mach number and density contrast,  we have opted to conduct only
medium
resolution
simulations. Because these are less time consuming, it allows us to
cover a large
range of
jet parameters.

Smith \etal\ (1985) have performed simulations of light ($\eta=0.1$)
jets with Mach
number
$M=6$ at different grid spacing, namely, $r_j/4$,  $r_j/8$ and $r_j/16$,
where $r_j$
is
the jet radius. At the lowest resolution the flow is under-resolved
although the
main
features are distinguishable. Doubling the resolution reveals many
details like the
oblique shocks in the jet and vortices inside the cocoon. Smith \etal\
found that
increasing the number of zones per jet radius to 16 unveils even more
detail but
they
conclude that 8 zones per jet radius gives a converged solution for the
global
properties
(see also K\"ossl \& M\"uller 1988). In this work we shall use this grid
spacing
($r_j/8$)
in both directions. Therefore, we have carried out simulations on a grid
with
dimension
$1280 \times 160$ corresponding to 160 jet radii along the symmetry axis
(which is
aligned
with the jet radius) and 20 jet radii in the perpendicular direction. In
two
examples (Runs A3 and B3) we have run the simulation using a larger grid
(1600 $\times$ 160).  We have used reflective boundary conditions along
the
jet axis, outflow elsewhere with the
exception of those grid cells that correspond to the jet inlet.

While the grid spacing determines how much detail of the flow is
revealed, the size
of the
grid along the jet axis is related to the question of how long to run
the
simulation. A
longer run time is necessary to establish whether or not, for example,
the system
eventually reaches a steady self-similar regime (see \S~\ref{sec-zh}).

We consider only pressure matched jets with adiabatic index $\Gamma =
5/3$. The two
main
dimension-less parameters are $\eta_c$ -- the ratio of the jet density
to ambient
density
-- and $M_c$ -- the jet internal Mach number
-- and are shown in Table 1.

The bow shock is easy to identify because of the clear jump of the
ambient gas pressure
and density. However, some artifice has to be used to locate the
boundaries of the cocoon
and jet. We do this by injecting tracer particles into the flow
throughout the width of
the jet inlet. Following the outermost particles up the point where they
reach the working
surface defines the jet boundary. Subsequently, the particles are
deflected back into the
cocoon, and its boundary, which should coincide with the contact
discontinuity, can also
be traced.  However, the instabilities which develop at the outer
boundary of the cocoon
make it difficult to determine its exact location. Thus, we have also
placed a second set
of tracer particles in the external medium to help us to demarcate the
contact
discontinuity.
We determine the bow shock and the borders of the jet and cocoon from
the positions
of the
tracer particles. We use these curves to find the radius of the jet
($R_j$), cocoon
($R_c$) and
bow shock ($R_b$) at a given distance from the central galaxy (jet
inlet).  For all
simulations and at 40 equally spaced positions along the z-axis, we
determine the
positions of the bow shock, cocoon and jet and calculate their expansion
velocity.
We also
calculate the volume and average energy density of slices perpendicular
to the jet
at all
grid points along the jet axis. These slices have width equal to the
grid spacing.

In the simulations presented here we have studied jets with zero opening
angle,
although
we have conducted a few test runs with non-zero opening angles. These
tests and
other
studies (see, e.g., Falle 1991; Komissarov \& Falle 1998) show that the
jets are
recollimated and the resulting jet structure is very similar to those
with zero initial opening angle.

%******************************************************************************

%******************************************************************************

\section{NUMERICAL RESULTS}
\label{sec-num}

\subsection{Numerical simulations in two and three dimensions}

Here we discuss differences in the results of 3D and 2D simulations with
special emphasis on intrinsic limitations to the results that can be
derived
from axisymmetric 2-dimensional simulations.

The great majority of jet simulations have been conducted in two
dimensions due
to both simplicity and computational power limitations. However, the
emergence of
supercomputers has made it possible to carry out 3D simulations with
resolution
comparable to the two dimensional, low-end resolution simulations. Over
the last
decade, a limited number of studies have now been conducted in three
dimensions
(Cox, Gull  \& Scheuer 1991; Basset \& Woodward 1995; Hardee, Clarke  \&
Howell
1995; Norman 1996; Hooda  \& Wiita 1996, 1998; Clarke 1996; Clarke,
Harris
\& Carilli 1997; Bodo et al. 1998; Aloy et al. 1999; Wang, Wiita
 \& Hooda 2000; Tregillis, Jones  \& Ryu 2001). They have shown
important
differences
between two and three-dimensional simulations. Here we shall comment on
those most relevant to the present discussion.

Basset  \& Woodward (1995) conducted a detailed analysis of jet
stability at high
resolution in two (non-axisymmetric, slab jets) and three dimensions.
They simulated
a short section of a long jet using periodic boundary condition and used
different grid
resolutions to test the convergence of the numerical process. They
find that the nature of turbulence differs in two and three dimensional
simulations. In two dimensions, there is an inverse cascade of energy to
large
scales, and small scale turbulence is not well developed. In the 3D
simulations
there is no inverse cascade, and the turbulence is developed to a
greater extent
with small scale turbulence present at low resolution. On the other
hand, they
conclude that the three dimensional simulations have instability modes
with similar
behavior to those of the 2D simulations.

Bodo et al. (1998) have also carried out three dimensional simulations
of
periodic jets and compared them with two dimensional simulations (Bodo
et al.
1994,
1995). They showed that the rapid growth of small scale structure in the
3D jet,
not seen in the 2D simulations, is responsible for most of the
difference
between the evolution of jets in two and three dimensions. Nevertheless,
they
conclude that the large scale morphological evolution of the jet can be
satisfactorily described by 2D simulations.

Norman (1996) has shown that there are some important differences in the
behavior of 3D simulations when compared with 2-dimensional ones -- due
mainly to
non-axisymmetric instabilities. For example,  as a consequence of the
different
character of turbulence in two and three dimensions, the turbulence in
2D
cocoons is mainly composed of large scale eddies while the 3D cocoons
are
characterized by smaller scale eddies. In 3D, strong backflow is
important only
in a region extending a few cocoon radii from the jet head, unlike the
2D
simulations where the backflow is strong over much of the cocoon length.

Entrainment of ambient gas is observed but, contrary to what happens in
2D,
the dense fingers do not penetrate all the way to the jet axis - which
contributes to greater stability of the jet in 3D. Thus, the 3D cocoon
can be
divided into a disturbed region near the working surface and a quiescent
region, which comprises most of its length.

The jet stability depends upon the conditions in the cocoon; and the jet
remains straight as it propagates through the quiescent region and is
deflected and deformed as
it penetrates the disturbed region. The wandering of the jet working
surface
(dubbed the ``dentist drill effect" by Scheuer (1982))
causes the jet to deposit its
momentum over an area much larger than its cross sectional
area. As we shall see later in the 2D simulations, an increase in the
working surface area is also observed but for a different reason.

The division of the jet into a stable and unstable section as a
consequence of the
existence of the cocoon quiescent and disturbed region, has also been
observed by Clarke
(1996). He has carried out 3-dimensional simulations of a light,
supersonic jet ($M=6$)
and observed the existence of three phases of jet evolution. If the jet
length is less
then $M_jr_j/2$ the jet behaves as a one-dimensional flow and it is
uniform and
featureless. Backflow is unimportant and the jet advance speed is given
by the ram
pressure arguments. This is called the {\it 1D phase}. When
two-dimensional effects such
as pinching mode instabilities and large cocoon eddies become important,
the jet is said
to be in the {\it 2D phase}. During this phase the jet length is less
then the so called
``disruption length" ($5M_jr_j$). In this phase, the jet remains
axisymmetric and
its general behavior is identical to that of the jets studied in the
present work.
Finally, when the jet length is larger
than the disruption length, it enters the {\it 3D phase}. The jet is no
longer
axisymmetric, it develops all modes of instability and becomes highly
disrupted.

Norman  \& Hardee (1988) performed simulations of slab jets to study
spatial stability and
have shown that if the jet is perturbed, it oscillates with a frequency
near the resonance
frequency and becomes disrupted. The results are in good agreement with
the linear
stability analysis carried out by Hardee \& Norman (1988). The
disruption length is $\sim
7$ times the $e$ folding length which  corresponds to an amplification
of the initial
perturbation by a factor of $10^3$. The disruption length, which depends
upon the Mach
number, can be expressed as
$
L_{dsrpt} \approx 10 M r_j.
$
Since the maximum length the jet in our simulations reaches is around
70-80 kpc,
with $r_j=500$ pc the effects of instabilities would be felt by jets
whose Mach
number is $\lapprox 14$.

The work by Clarke (1996) shows that the length of time the jet remains
in the 2D phase
increases with the Mach number. As we can see in Table 1, most of the
jets representing
the main body of our simulations (Run A) have relatively high Mach
number. Thus, given the
large ratio between the maximum size of the jet and its radius ($\sim
140 - 160$) the
assumed axisymmetry does not represent a serious limitation except where
$M$ is less than
$\sim 10$. In low Mach number jets (e.g., Run A6 and A7), the
two-dimensional
simulation will not significantly depart from a full 3D simulation until
the jet reaches
about one quarter of the grid size. As the jet lengthens however, the
differences would be
accentuated mainly in the cocoon width and turbulent structure, vortex
shedding, the
character of jet instabilities and of the working surface. {\it
Therefore, the
results regarding the behavior of the head radius, cocoon volume and
energy distribution as a function of jet size in low Mach number
($M \lapprox 10$) and large source sizes must be viewed with caution.}


On the other hand, it is interesting to note that, the disruption of the
jet and the
increase of the working surface area observed by Clarke (1996) in three
dimensions during
the 3D phase has its correspondence in the 2D simulation -- at least for
low power jets. We
observe various degrees of disruption in low Mach number jets due to the
development of
internal shocks (\S \ref{sec-jetpower}). This occurs for $M \le 3$ and
in a less
dramatic fashion for $3 < M <
10$. For example, the present simulations show that, for low Mach
number, oblique internal
shocks dominate the jet flow, and the energy transfer to the working
surface is highly
disrupted. One is tempted to say that, in two dimensions, high Mach
number jets behave as
2D phase jets whereas low Mach number jets have a behavior closer to
that of 3D phase
jets. In fact, the possibility of a unified model for the FR I and FR II
radio sources
based upon their association with, respectively, the 3D and 2D phases,
as suggested by
Clarke, has its parallel in 2-dimensional simulations and is explored in
Section
\ref{sec-jetpower}.

Although a full 3-D investigation would give us somewhat more
confidence, we believe that,
at least for high Mach number, it would not significantly change the
general
picture emerging from the present results. Though, for low Mach number,
we can trust the
results only while the jet size is smaller than a few tens of
kiloparsec. On the other
hand,  as the jet size increases, the high degree of disruption observed
in
the axisymmetric simulations exhibits many similarities with the 3D
results.

One of the main goals of the present work is to make a critical analysis
of
existing self-similar models and to construct an analytic model of
powerful FR II type
radio sources based on the numerical simulations. In this respect, two
dimensional,
axisymmetric simulations seem quite convenient since these models also
have cylindrical
symmetry. In addition, the great advantage of the more economical 2D
simulations is that
they allow us to explore a wider range of jet parameters and different
environments.

%******************************************************************************

%******************************************************************************

\subsection{Structure and Dynamics of the Simulated Jets}

Here we present selected results of the simulations which illustrate the
global behavior
we shall incorporate into our analytical model (Paper III).
Some of the results shown here have been discussed in previous work.
However, we present
them here in order to emphasize new aspects of the jet expansion
dynamics.  In addition,
we want to clarify some of the results and/or arguments used in the
construction of the
analytical model in Paper III.

The simulations are carried out using dimension-less quantities.
However, we present the
results with a scaling to physical units in order to facilitate close
comparison with
astrophysical radio sources. Thus, the jets presented here have radius
500 pc and the grid
size is 80 kpc $\times$ 10 kpc while the grid spacing is 62.5 pc. (For
Runs A3 and B3 the
grid size is 100 kpc $\times$ 10 kpc.) The density and temperature of
the ambient medium
are 1 cm$^{-3}$ and $\sim 4.8\times 10^7$ K, respectively. Since the jet
is pressure
matched, giving its Mach number and density ratio completely  defines
the system. Since
the velocities $v_j$ of the equivalent classical jets listed in Table 1
are
non-relativistic ($0.013 < \beta_c < .36$), the jet power can be
estimated from the
expression
$$
L_j \sim \pi \rho_j v_j^3 r_j^2.
$$
We have calculated $L_j$ using the above values for the jet radius
$r_j$, ambient gas
density and temperature and $\eta$ and $M$. The results are shown in
Table 1 and we see
that the jet power lies in the range $2.9 \times 10^{44}
- 5.9 \times 10^{48}$ erg\  s$^{-1}$.

We have performed the simulations for the parameters shown in Table 1
and in Fig.
\ref{ASUA_All1} we show the density contours at the time the shock front
has nearly
reached the end of the grid, for Runs A1-A5.

The size of each panel in the figures is 80 kpc $\times$ 20 kpc. Fig.
\ref{ASUA_View1}
shows the time evolution of simulations for jets with $\gamma_j=2$, 2.5,
5 and 7 in Case 1
(Runs A1-A4). In Fig. \ref{ASUA_View2} we show the corresponding
simulations for Case 2
(Runs A6-A9). The time in the panels is given in Myr. All the basic
features present in
previous simulations are also seen here. The most striking feature is
the narrowing of the
bow shock with increasing jet speed. This effect can be seen more
clearly in Runs B1-B5
where we have fixed the density ratio $\eta=0.1$ and varied the Mach
number from $M=2$ to
$M=20$. While the cocoon width remains more or less unchanged, the bow
shock varies from a
very broad to a more narrow nearly cylindrical shape. If we plot the
density contours for
Runs B6-B10 for which the jet speed is fixed ($M=6$) but the density
varies over two
decades between $\eta=0.003$ and $\eta=0.3$, we find that when the jet
density increases,
the cocoon also narrows and becomes less prominent.

As is well known, the back flow of jet material relative to the shocked
ambient gas causes
the contact discontinuity to become Kelvin-Helmholtz unstable (e.g.,
Norman \etal\ 1982).
By the time the jet has reached the end of the grid (and especially for
low Mach number)
these instabilities appear to be fully developed. However, as the jet
speed increases we
find that the instabilities are under-developed and the flow of the
shocked jet material
inside the cocoon seems less chaotic. This is very clear when we compare
Fig.
\ref{ASUA_All1}d,e with, e.g., Fig. \ref{ASUA_All1}b. This implies that
most of
the kinetic energy is in the form of large scale motion instead of
turbulent eddies.

In those runs where the instability is well developed, that is, for low
and moderately
large Mach number, there is strong mixing/entrainment of the shocked
ambient gas in the
roughly half of the source nearest the jet inlet. As a result, the
contact discontinuity
surface becomes ill-defined which makes it difficult to accurately
determine the cocoon
boundary. That is the reason why in Fig.'s \ref{ASUA_View1} and
\ref{ASUA_View2}, as time
increases, the cocoon seems to have collapsed in the runs with $M=27.4$,
9.4 and
$M=13.3$.\footnote{We shall take this into account when calculating the
cocoon volume in
\S~\ref{sec-Vol}.}

%******************************************************************************

%******************************************************************************

\subsection{Internal Shocks and Jet Knots}

It has been suggested that the knots observed in jets such as Cygnus A
and M 87 may be
caused by internal shocks (e.g., Bicknell \& Begelman 1996 and
references therein). In
Fig. \ref{ASUA_All1} we can clearly distinguish the series of internal
oblique shocks that
are formed in the jet. In particular, we observe that the distance of
the first shock from the
nucleus (jet inlet) increases with increasing Mach number (Runs B1-B5).
On the other
hand, it seems
not to depend  strongly upon the jet density ratio (Runs B6-B10). One
can use the argument
presented by Wilson \& Falle (1985) to show that this distance is
proportional to the jet
speed. Wilson and Falle suggested that the shock can be induced by
changes in the external
pressure. If so, the information about any such a change propagates in
and out of the jet
at the sound speed. During this time, the jet travels a distance given
by
\beq
\ell \simeq 2r_jM.    \label{ele}
\eeq

We can tentatively identify the jet response scale with the position
$z_s$ of the first
shock.\footnote{We have examined the jet density along its axis as a
function of the
distance from the central galaxy for Runs A1-A10 and B1-B10. Although
there is a quite
complex series of secondary shocks, it is possible to clearly identify
the first shock.}
We have estimated $z_s$ by determining the position of the first most
significant maximum
of the jet density. Fig.'s \ref{ASU_Shock}a and c show the values of
$z_s$ plotted against
the jet response scale $\ell$ for Run B and A respectively. The value of
the jet radius
used in Eq. (\ref{ele}) is the radius of the jet inlet in the
simulations, that is,
$r_j=500$ pc. The straight line represents $z_s = \ell$. The points tend
to lie near the
line and we see that the agreement is better for small Mach number. As
the Mach number
increases, the distance of the first knot becomes much smaller than that
predicted by
(\ref{ele}). This can be attributed to the fact that the average radius
of the jet tends
to decrease with increasing Mach number. In section \ref{sec-Rh} we
examine the jet radius
at its terminal point ($r_h$) as a function of the source size (see Fig.
\ref{ASUA_Rh_Xh})
and we find  that it is smaller than the initial value (500 pc) for high
values of $M$.
This is generally true for the entire jet - not just at the terminal
point. Thus, we
determined an average value of the actual simulated jet radius over the
distance between
the jet inlet and the jet response scale. Substituting this corrected
value of the jet
radius in (\ref{ele}), a new value of $\ell$ is obtained which is used
to plot Fig.
\ref{ASU_Shock}b and d. We see that, now, the agreement is much better
for the relatively
high Mach number jets in Fig. \ref{ASU_Shock}d but is relatively
unchanged for the
moderate $M$ jets in Fig. \ref{ASU_Shock}b. Thus, based on our broad
exploration of
parameter space, we are able to confirm the suggestion of Wilson \&
Falle that the
location of the first knot can be used to give a crude estimate of the
jet Mach number
which is independent of the density ratio.

We have applied the expression (\ref{ele}) to the jet in Cygnus A in
order to estimate its
Mach number.  According to Perley, Dreher \& Cowan (1984), there are at
least two regions
in the jet: a region of large opening angle up to knot 2 at a distance
2'' from the
nucleus and a region of constant width between knots 2 and 4 (from 2''
to 20'') with a
FWHM = 0.7'' (Carilli et al. 1996). We suppose that the position of the
jet inlet in our
simulations corresponds to knot 2 where the jet becomes collimated with
a zero opening
angle. One can then apply (\ref{ele}), e.g., to knot 3, assuming that it
is the
first most prominent internal shock. For Cygnus A we have 1'' $\sim$ 1
kpc so that, taking
the FWHM as a measure of the jet diameter gives $2r_j\sim 700$ pc and
measuring the
separation between knots 2 and 3 we obtain $M\sim 13$. This is roughly
consistent with the
results obtained using different approaches by Williams (1991) ($M\sim
8$) and  Carilli
\etal\ (1996) ($M\sim 8$ and $M<13$).


%******************************************************************************

%******************************************************************************

\subsection{Pressure and Expansion Speed}

Here we discuss the spatial dependence of (1) the pressure in the cocoon
and the shocked
ambient gas (SAG) region, that is, the volume between the cocoon and the
bow shock and (2)
the lateral expansion velocity of the cocoon and bow shock. In Fig's.
\ref{ASUA_Pres1} and
\ref{ASUA_Pres2} we present the pressure profile in the direction of the
jet axis for the
cocoon (thick line) and in two regions inside the SAG region: (1) the
outer SAG region
defined as $R_b/2 < R <R_b$ (thin line) and (2) the inner region defined
as $R_c<R<R_b/2$
(dashed line). The pressure is averaged over slices perpendicular to the
z-axis. The
horizontal axis represents the distance to the central object (jet
inlet) and the vertical
axis is the logarithm of $P$ relative to the ambient pressure $P_a$. In
those sections of
the source where the bow shock has left the simulation grid, we plot
only the pressure in
its inner region. We see that the pressure is always high near the head
as one would
expect and falls with distance from the head. This is generally true for
the three curves
although they differ in some details. The pressure falls more rapidly
inside the cocoon
than inside the SAG region and the cocoon only stays over-pressured near
the jet head
except in very powerful jets (Runs A3 and A4). At intermediate  Mach
number
($M \sim 10$) the cocoons are in pressure balance with the ambient
medium over much of the
length of the source; while at higher ($M>20$) or lower ($M< 5$) Mach
number, the cocoons
are over pressured or under pressured, respectively. At high  Mach
number, the SAG region
maintains a pressure higher than the ambient pressure over much of its
length -- most
probably because it is under expanded. We also find a pressure gradient
perpendicular to
the jet axis -- i.e., away from the high pressure head, the cocoon and
the inner SAG
region are generally in pressure balance, while the outer SAG region is
at higher
pressure.


Note that the pressure in the cocoon drops with distance from the head
faster
in Case 2 than in Case 1, although the width of the cocoon is
comparable.
This occurs  because in Case 2, the total internal energy inside the
cocoon
is much smaller than in Case 1 since the power supplied by the jet is
lower.
Since the cocoon pressure near the head is relatively high in both
cases,
the pressure has to drop faster in Case 2, in order to yield a smaller
volume integral (total internal energy).


Because of the outflow boundary conditions, when the bow shock reaches
the sides
of the grid (at late times in low Mach number jets), it will continue to
move
outwards and will be lost.
This raises the issue of whether this will introduce a bias in our
estimates
of the volume limited average pressure in the corresponding shocked
ambient
gas (SAG) regions.
Thus, we have estimated the signal crossing time $T_c$ in those slices
of the
SAG region  where the bow shock has left the grid. We have compared this
to the time $T_g$ which these regions spend in the simulation once the
bow
shock has left the grid.
Thus, for a given slice of the source, $T_g$ is the age which
corresponds to
the total age of the source at the end of the simulation minus the age
of the
source at the time  the bow shock for that slice left the grid.
We find that for high and intermediate Mach number, $T_c$ is always
larger or of the order of $T_g$. For low Mach number, by the time the
source
has reached the horizontal end of the grid, $T_c$ is smaller than $T_g$.
However, up to the moment the source is approximately half its final
size, the signal crossing time is at least $T_g/2$.

In addition, as expected, the pressure in the SAG region can become
smaller
than the ambient pressure $P_a$ in the very old, relaxed regions which
 increases the local signal crossing time. We conclude that, since in
the majority of the studied parameter space, $T_c$ is not sufficiently
 smaller than the local SAG age to induce  a significant change in the
pressure due to the fact that the bow shock has left the grid.

Also, the pressure profiles in Figures 6 and 7 and in Figure 5 of Paper
II,
show that the behavior of the pressure in the inner SAG region (dashed
line)
does not exhibit a significant difference from the expected behavior as
one goes from near the head, where the bow shock is still inside the
grid,
to near the jet inlet, where the bow shock has already left the grid.
In general, the pressure falls smoothly from a value many times the
ambient
pressure $P_a$ to a value of the order of $P_a$. This is an indication
that
there is not a substantial bias in the estimated value of the volume
limited average
pressure in the regions under consideration.


A plot of the lateral expansion speed of the cocoon normalized to its
maximum value (which occurs near the jet head) shows that initially it
rises and then falls sharply and smoothly away from the head. At later
times it starts varying randomly as Kelvin-Helmholtz
instabilities develop at the contact discontinuity. Excluding a small
portion near the head, the average expansion speed of the cocoon is
nearly
zero and does not change appreciably for most of its length.
Consequently,
the cocoon assumes an essentially slowly expanding cylindrical shape.
The deceleration of the cocoon lateral expansion is also
reported by Reynolds \etal\ (2001a).

Similar to the cocoon, the bow shock lateral expansion speed initially
increases at the
head after which it falls sharply. Away from the head, the bow shock
lateral expansion
speed decreases very slowly and remains above zero -- giving the bow
shock its nearly
conical shape -- with an opening angle which decreases at high Mach
number.

In Figure~\ref{ASUA_Vbc_x} we show the ratio of the lateral expansion
speed of the bow
shock and cocoon boundary to the ambient sound speed (i.e., the sound
speed of the medium
it is expanding into). The bow shock is initially expanding
supersonically near the jet
head, but rapidly slows to a quasi-sonic velocity in the jets with
moderate Mach number.
(The lateral expansion speed of the high Mach number jets remains
supersonic during our
simulations). On the other hand, the cocoon is expanding laterally
subsonically (in the
post bow shock medium) over nearly its entire length.

\subsubsection{Implications for Emission Line Nebulae in Radio Galaxies}

\label{sec-clouds}

There is much interest in using the properties of optical emission line
 clouds in radio galaxies to probe the nature of the interaction between
the radio source and its environment (e.g., Heckman, van Bruegel \&
Miley
1984; Baum, Heckman \& van Bruegel 1992; Neeser, Meisenheimer \&
Hippelein
1997; Koekemoer \& Bicknell 1998; De Vries et al. 1999; Villar-Mart\'in
\etal\
1999; Best, R\"ottgering \& Longair 2000; Tadhunter et al. 2000).
The current picture of cloud acceleration by the radio source is given
by,
e.g., Klein, McKee \& Colella (1994), Anderson \etal\ (1994), (and see
the
discussion in O'Dea \etal\ 2001). The dense clouds in the ambient medium
will
be overrun first by the bow shock which sends a shock into the cloud
and due to the pressure difference across the cloud, accelerates the
cloud to the
velocity $v_{cl} \simeq v_{b}\chi^{-1/2}$, where $v_{b}$ is the bow
shock expansion
velocity and $\chi$ is the density ratio of the cloud and the ambient
medium.
Once the bow shock passes, the cloud is entrained in the wind-like
expanding
post-shock medium and accelerated further. Klein \etal\ estimate that
the
cloud is accelerated to a large fraction of the wind velocity in a
timescale
which is a few times the shock crossing time of the cloud (they call
this
the cloud crushing time) $t_{cc} \simeq a \chi^{-1/2}
v_{b}^{-1}$ where $a$ is the cloud radius. However, the observed line
offsets and/or line FWHM in the emission line nebulae are in the range
100-1000 \kms\ and are 1-2 orders of magnitude slower than the bow shock
initial velocity. Thus there must be some effect which
limits the velocity to which the clouds can be accelerated.

As we have found here, the properties of the radio source SAG region and
cocoon -
especially the pressure and lateral expansion velocity -- are very time
and position dependent. This time dependence will have a strong effect
on the resulting emission line cloud ionization and kinematics. In
particular,
over most of the life time of the radio source, any overrun or entrained
clouds will be embedded in a medium which is moving close to the sound
speed and has a pressure which is close to the undisturbed ambient
pressure.
Thus, the ultimate velocity of any ambient clouds may be closer to a few
hundred to 1000 \kms\ rather than to the much higher ($\sim 10^{4-5}$
\kms) propagation speed of the radio source.


%****************************************************************************

%****************************************************************************

\section{SELF-SIMILAR MODELS}
\label{sec-similar}

An important question is to what degree (and at what stages) the source
propagation and evolution can be characterized as self-similar (e.g.,
Falle 1991; Komissarov \& Falle 1998; and see \S~\ref{sec-zh}).
In this section we present a brief review of the main self-similar
analytical models for radio source propagation (see also Daly 2001). In
section \ref{sec-zh} we compare the results of our simulations with
these models. It turns out that the existing models fall into three
main categories. The properties of the models including the time
dependence of the main quantities that characterize the source
 expansion are given in Table 2.

The starting assumption of all the models is that there is equilibrium
between the pressure in the advancing head of the jet and the ram
pressure
exerted by the ambient gas.  In these models, the jet kinetic energy is
thermalized at the jet head -- inflating an expanding over-pressured
cocoon.

The simplest model is Type I in which the jet advances with constant
speed in a constant density atmosphere. This is achieved by supposing
that the jet has a zero opening angle and that the head pressure is
kept constant. The average global pressure inside the cocoon
is proportional to $t^{-1}$ since its volume increases as $t^2$
while there is a constant energy input from the jet, causing its total
energy to increase at a rate proportional to $t$ (Begelman \& Cioffi
1989). Therefore, the ratio between the head and global cocoon
pressure varies as $t^{-1}$.

In the Type II model the ambient density falls as $d^{-2}$ but, because
the jet is assumed to have a non-zero constant opening angle, the head
pressure varies as $t^{-2}$ and the head advance speed is constant
(Daly 1990). Interestingly, the cocoon expansion is
identical to that of the Type I model. This is because the density,
as seen by the expanding cocoon, is now taken to be constant because
the distances are small compared with the jet size, that is,
$r_c\ll z_h$. This is not rigorously true in an atmosphere with density
varying as $d^{-2}$. However, if the density profile is of the type of
King's
law, for distances smaller that or of the order of the core radius,
the density does not change appreciably. Therefore, assuming that the
cocoon radius is not much larger than the density distribution core
radius,
the external density can be approximated as constant.

The most general model is Type III which assumes a power law density
distribution with an arbitrary exponent $\delta$. Here, the basic
assumption is that, as the jet head advances, it does so in such way as
to keep the ratio between the head and global cocoon pressure
constant (Falle 1991; Begelman 1996). All parameters varies with time
in a manner that depends on the exponent $\delta$.

The three types are independent in the sense that none is a special
or limiting case of another. For example,  for $\delta=0$, Type III
does not reduce to Type I. Type III also does not reduce to Type II
for $\delta=2$. However, the only difference between Type II
and Type III for $\delta=2$ is the behavior of the cocoon  -- apparently
due to the difference in the external density variation in each model.
In Type III, the density varies on all scales as $d^{-2}$; while in
Type II it is taken to be constant on the small scales of the lateral
expansion of the cocoon. We note that the three types of models give
different predictions for the ratio of the head and global cocoon
pressure - which we shall compare with the results of our simulations
in Paper II.


%******************************************************************************

%******************************************************************************

\section{GLOBAL PROPERTIES OF THE SIMULATIONS}
\label{sec-global}

In this Section we describe the results from our simulations on the
behavior of the source size, head radius, advance speed, and
source volume. We compare them to the predictions of the self-similar
models. Our results are summarized in Table 3.

\subsection{Is the  Jet Propagation Self-similar?}
\label{sec-zh}

Falle (1991) and Komissarov \& Falle (1998) suggested that even if jets
did not initially propagate in a self-similar manner, they would
eventually
evolve to do so. They derived the following criterion to determine
{\it a priori} when, if ever, the jet expansion will
become self-similar. The transition to a self-similar propagation should
occur when the jet length is much larger than a characteristic length
$\ell_c \sim
\sqrt{L_j/c^3\rho_a}$
where $L_j$ is the jet luminosity and $\rho_a$ the ambient gas density.
Using, $ L_j
\sim \rho_j v_j^3 r_j^2$ we obtain
$$
\ell_c \sim \sqrt{\eta \beta^3} r_j.
$$
In their simulations, Falle (1991) and Komissarov \& Falle (1998) used
jet parameters such that $\ell_c$ is slightly larger than the jet 
radius. They show that the self-similar
regime is only attained when the jet has propagated a distance equal to
a few tens of $\ell_c$.

For the simulated jet properties listed in Table 1 we estimate a
characteristic
length 0.4
pc $< \ell_c < 60$ pc. This is much smaller (by $\sim 1-2$ orders of
magnitude) than
the
radius of the jet used in these simulations (500 pc). Therefore, we
would expect
that as
long as the length of the simulated jet is a few times its radius it
will display
self-similar behavior. Our grid is 160 times the jet radius and, in the
worse case,
that
is, for the most powerful jet, it is greater than $10^3 \ell_c$. We also
note that
the jet
of Cygnus A is believed to be $\sim 1$ kpc in diameter while the overall
projected
size of
the source is about 120 kpc, giving a projected length to diameter ratio
for each
jet of
$\sim 60$.  Therefore, by the above criterion, our simulations should be
large
enough to
address the question of self-similarity in sources with typical observed
properties.


In Fig. \ref{ASUA_Xh_Tt} we show the source size (the distance from the
source head
to the
nucleus (jet inlet)) $z_h$ as a function of time. In order to determine
whether  the
source has reached a self-similar stage we compare the slope of this
plot with the
predictions of the Type I and III models. Figure  \ref{ASUA_Xh_Tt} shows
a solid
line
representing the least squares fit of the function $z_h\propto t^m$ for
the points
with
$z_h > 10$ kpc. The dashed  and dotted lines denote
 $z_h\propto t^{3/5}$ (Type III) and $z_h\propto t$, (Type I, i.e.,
constant advance speed). We see that, the fit lies near the two
theoretical
predictions.
However, a single fit over all points may not be the best way to
determine the
exponent
since it can vary with source size. In fact, in order to determine
whether the
system is
evolving towards a self-similar regime, we will examine how the exponent
$m$ changes with time.

We have devised two different ways to do this, which give consistent
results. (1) We first
fit the exponent using all the points (in Fig. \ref{ASUA_Xh_Tt}). Next
we start at the
second point from the jet inlet and fit all the remaining points. Next
we start at the
third point from the jet inlet and fit all the remaining points, and so
on, each time
fitting a decreasing number of points and starting the fit further from
the jet inlet. The
process continues until there are only 5 points left to fit. The
exponents obtained in
this way are plotted as a function of the position $z_h$ of the starting
point.  Thus, for
smallest value of $z_h$ the slope $m$ represents a global average value
and gives no
information on possible local variations. As one approaches the highest
values of $z_h$
and, if the system tends to a self-similar phase, $m$ should approach
the asymptotic
value, either 1 or 3/5 depending on which model best describes the
simulated jet.
(2) In the second method, the exponent $m$ is estimated by calculating the
logarithmic derivative
of $z_h$ with respect to $t$, that is, $d(\log{z_h})/d(\log{t})$. After
applying a moving
window average around each point, we plot $m$ as a function of $z_h$.
This gives us
information on the local variations of $m$. Again, we expect that, as we
come closer to
the outermost points, the exponent should approach the self-similar
limit if the jets
become self-similar before the end of the simulation. Fig.
\ref{ASUA_Similar_Xh}
shows the
results for Runs A1-A10. The
derived values of $m$ are denoted by open circles and a solid line for
the first and
second methods, respectively.

The simulated jets do not uniformly favor any particular self-similar
model.
Instead, the
derived values of the exponent $m$ tend to be bounded by the predictions
of the Type
I and
III models. There is some indication (Fig. \ref{ASUA_Similar_Xh}) that
the lower
power
jets tend to have exponents closer to 0.6 (i.e., decelerating) and the
higher power
jets
have exponents closer to (or even greater than) unity (i.e., propagating
with nearly
constant velocity). We also  find that in a given simulation, the
exponent tends to vary
as a function of source size. However, in many runs, even after the
source has propagated
$\sim 160$ jet radii, the value of the exponent $m$ does not seem to
have settled
down completely. Exceptions are Runs A5, A9, B2, B5, B7 and B8. At
constant
$\eta=0.1$ the asymptotic value of $m$ tends to decrease from $\sim 1$
to 3/5 as the Mach number increases, although Run B4 is still unsettled.


The behavior of the source size as a function of size is complex and
sources appear to go through different ``phases" which depend on their
Mach number.
In \S~\ref{sec-Vh} we examine the advance speed of the source, and
find that, for high Mach number, the advance speed remains  more or
less constant for the length of the simulation. For moderate Mach
numbers, the advance speed starts decreasing long before the source
reaches the end of the grid. Finally, at very low $M$, the advance
speed remains almost constant -- near the ambient sound speed. Thus
the exponent $m$ approaches 1 for very low and also for very high
values of $M$, while for moderate $M$ it oscillates between 1 and 3/5.

We thus suggest that the source expansion may go through to up three
distinct phases, depending on the jet initial Mach number. The head
of a highly supersonic jet advances first at an almost constant speed
(i.e., the exponent $m=1$) and eventually starts decelerating
(due to the increase in size of the jet head with time -
\S~\ref{sec-Rh}), thereby entering a second phase. During this phase,
the jet advance decelerates and $m$  approaches the value
3/5. Much later, its advance speed will reach the ambient sound speed
after which it remains constant (again $m=1$) until the jet is
eventually disrupted (see \S~\ref{sec-jetpower}).
Moderately supersonic jets start decelerating earlier, as soon as they
have
propagated a few tens of jet radii and, therefore, enter phase two
quickly (with $m$ approaching 3/5).  Later, when its speed reaches
the external sound speed it enters phase three with $m=1$.
For very low Mach number jets the propagation speed is almost sonic
from the start so the source enters the ambient medium with a constant
speed and, therefore, is essentially already in phase three ($m=1$).

We note that Komissarov \& Falle (1998) observed what they call an
``intermediate phase"
during which the bow shock advances at nearly a constant speed before it
reaches a
self-similar regime when $z_h \propto t^{3/5}$. This phase, which in
their simulations can
be very long, may be interpreted as phase one described above, with the
jet later making a transition to our phase two.

We conclude that within the parameter space of the simulations carried
out here
(which we have argued include parameters typical of powerful radio
galaxies),
the behavior of the source size as a function of time is not always
self-similar. Instead, the simulations show more complex behavior,
in some runs with transitions between different limiting values of the
exponent of size vs. time.


%******************************************************************************

%******************************************************************************

\subsection{The Head Radius}
\label{sec-Rh}

Here we discuss the size of the working surface of the jet  and its
dependence on
source
size and Mach number. From Fig. \ref{ASUA_View1} and \ref{ASUA_View2} we
see that the
radius of the jet at the working surface as well as the cocoon radius at
$z_c$, which
coincides approximately with the position of the working surface, vary
as the source
expands. This can be seen more clearly in Fig. \ref{ASUA_Rh_Xh} where we
have plotted both
quantities as a function of the source length $z_h$. Note that there is
a transition (or
transient) regime between the jet inlet and $\sim 20 r_j$ ($\sim 10$
kpc).

We notice a quasi-periodic modulation of the head radius that can be
attributed to
``vortex-shedding" as discussed by Norman et al (1982), K\"ossl
\& M\"uller (1988), Falle (1991) and also observed by Cioffi \& Blondin
(1992). Vortices are generated near the head as the jet material is
deflected
sideways and
are periodically shed from the head into the cocoon. According to Falle
(1991), for
large
jet speed, the vortex is not immediately swept away into the cocoon by
the back flow
and starts growing in size by accumulating vorticity. In this process,
the size of head
consequently increases. As the vortex core becomes diffuse it is finally
swept away
from the head, the head size decreases and a new vortex begins to form.
This process
modulates the bow shock such that in many runs it displays an undulating
shape (Fig.
\ref{ASUA_All1}; and see also Fig. 11 of Komissarov \& Falle 1998).
A detailed description
of the vortex shedding phenomenon and how it
modulates the jet advance speed (see section \ref{sec-Vh}) is given by
Lind et al. (1989).


In addition to (or in spite of) the quasi-periodic oscillation, there
are trends in the
average value of the radius of the jet and cocoon near the head, (which
we shall call
generally $r_h$) as a function of both source size and Mach number.  In
order to examine
this, we plotted separately the average value of $r_h$ for Runs B1-B5 as
a function of the Mach number (constant $\eta$) and for Runs B6-B10 as a
function of $\eta$
(constant M). We have also included the third set of simulations,
referred to as Runs C1-C6 in Table 1 which differ only slightly in their
input parameters.

We observe that that both the jet and cocoon head radii remain constant
for low and intermediate $M$ and decrease at high $M$. In addition, the
rate of increase of the jet and cocoon head radii with source size
(growth rate) is negative for low and intermediate $M$ and
positive for high $M$. On the other hand, at $M=6$ the average
value of the jet and cocoon head radius and their growth rates do not
depend significantly on the density contrast $\eta$.

Since we have established that the behavior of $r_h$ against $M$ does
not depend strongly
upon the value of the jet density ratio, we include in the analysis the
results of Runs
A1-A10 which encompass a larger range of Mach numbers. In Fig.
\ref{ASU_Rmed_ABC} we plot
the results of all 26 runs combined illustrating the dependence of the
head radius on the
Mach number. The cocoon and the jet terminal radius are plotted in Fig.
\ref{ASU_Rmed_ABC}a as open circles and triangles, respectively while
the growth rate of the
radius is shown in Fig. \ref{ASU_Rmed_ABC}b. We find that $r_h$ remains
more or less
constant for $M\lae 10$ and then decreases for large Mach number. The
growth rate of the
head as a function of source size is negative for $M < 10$ and positive
for $M\gae 10$. A
possible explanation for the non-linear behavior of the jet radius with
$M$ is the
following. At high Mach number, downstream from the re-confinement
point, the pressure
surrounding the jet head is very high and the jet is strongly focused
due to what
Komissarov \& Falle (1998) call ``super-collimation". This causes the
average radius to
decrease with increasing $M$. Because the pressure scales as $M^2$, this
process, which is
very effective at high Mach number, becomes less efficient at low Mach
number ($M<10$) and
the jet head radius remains approximately constant at low $M$.


%****************************************************************************

%****************************************************************************

\subsection{The Advance Speed of the Source}
\label{sec-Vh}

Here we discuss the behavior of the advance speed of the source as a
function of
source
size and examine the dependence on head radius. We show the advance
speed of the
source as
a function of source size in Fig. \ref{ASUA_Vh_Xh} for Run A. After the
transient
region that
extends to $\sim 10$ kpc, the advance
speed decreases slowly for low Mach number and stays approximately
constant at high
Mach
number.

Assuming ram pressure balance between the head of the source and the
ambient medium, the velocity of the jet head is given by (Cioffi \&
Blondin
1992; Loken et al. 1992; Hardee et al. 1992; Frank et al. 1998)
\beq
v_h=v_j \left({\mu\sqrt{\eta}\over 1+\mu\sqrt{\eta} }\right),
\label{vh}
\eeq
where $$\mu=r_j/r_h$$ and $$\eta =\rho_j/ \rho_a.$$ We calculated the
required
radius of
the working surface $r_h^*$, such that, for a given jet radius $r_j$ and
density
contrast
$\eta$, the jet advances at the ``observed'' speed shown in Fig.
\ref{ASUA_Vh_Xh}.
The result is
shown as a solid line in Fig. \ref{ASUA_Rh_Xh} for Run A.
We see that $r_h^*$ closely follows both the jet and the cocoon radius
at the working surface.

Taking a complementary approach,  we have also calculated the expected
advance speed
using
Eq. (\ref{vh}) and adopting the values of jet and cocoon radius  $r_h$
derived from
the
simulations (Fig. \ref{ASUA_Rh_Xh}). In Fig. \ref{ASUA_Vh_Xh} we
compare the results with the observed speed. The solid and dotted lines
represent
the
advance speed calculated using the cocoon and jet radius, respectively.
Except for
low
Mach number ($M \lae 10$), we see that they agree well with the measured
speed,
showing
that Eq. (\ref{vh}) gives a good description of the head motion for high
$M$ jets.
For
small values of $M$, the advance speed is small and, consequently, the
ram pressure
is
very low, of the order of the external pressure. Thus, the assumption
that that the
external pressure is negligible compared to the ram pressure breaks down
and Eq.
(\ref{vh}) is no longer valid.

In summary, we find that (1) low $M$ jets decelerate in a constant
density
atmosphere,
while high $M$ jets propagate at nearly constant velocity; (2) the
oscillations of
the
head area modulate the speed of the head; and (3) the variation of the
head speed is
anti-correlated with the head radius consistent with the hypothesis of
balance
between the
pressure in the head and the ambient gas ram pressure (Begelman \&
Cioffi, 1989).

%****************************************************************************

%****************************************************************************

\subsection{Source Volume}
\label{sec-Vol}

Here we examine the time dependence of the source volume ($V\propto
t^n$) and determine to
what extent the volume behaves in self-similar manner. The simple model
proposed by
Begelman \& Cioffi (1989) predicts that the overall volume of the cocoon
increases with
time as $t^2$ (Type I models) while in the improved models of e.g.,
Begelman (1996),
Bicknell, Dopita \& O'Dea (1997) and Kaiser \& Alexander (1997), the
volume is
proportional to $t^{9/5}$ (Type III models). In order to test these
predictions, we have
determined the volume of both the cocoon and of the region limited by
the bow shock as a
function of source size in our simulations. Figs. \ref{ASUA_Vb_T} and
\ref{ASUA_Vc_T} give
the results for the whole source and cocoon respectively. We see that
the expansion of the
source is generally consistent with the  prediction of the models; while
the cocoon volume
tends to increase more slowly than predicted.

However, as discussed in \S~\ref{sec-zh}, it is useful to examine the
dependence of the
exponent $n$ on the source size in order to determine whether the source
evolves to a
self-similar regime. Therefore, we have used the same method employed in
section
\ref{sec-zh} to examine in more detail how the exponent $n$ varies. The
results are
presented in Fig. \ref{ASUA_Similar_Vb} and Fig. \ref{ASUA_Similar_Vc}
for the whole
source and cocoon respectively.

Unfortunately, meaningful results are only
possible for the whole source as long as
the bow shock  remains within the grid boundaries. This is true  only
for the high Mach
number jets where we see that the exponent $n$ is in  good agreement
with the predictions
of Type I models ($n=2$) for $M\gapprox 27$ (Runs A3, A4 and A5 - Case
1) and above the
predicted value for $M=13.3$ and $M=19.2$ (Runs A9, A10 - Case 2) -
i.e., the source
volume increases faster with increasing size in Case 2. These results
are, to some extent,
independent of those obtained for the source size.

We see from Fig. \ref{ASUA_Similar_Xh} that in only two of the above
examples
(Runs A4 and
A5) is the size of the source proportional to $t$ as required by Type I
models. Therefore,
the fact that $n=2$ does not always mean that the bow shock radius is
proportional to
$t^{1/2}$ as predicted by the models -- i.e.,  the source volume
increases in an
approximately self-similar way even though the source length and the bow
shock radius may not.

The behavior of the cocoon is more complex due to the K-H instabilities
which make the definition of the contact discontinuity more 
problematical and the exponent $n$ exhibits
large oscillations as a function of source size. However, we do note (1)
the cocoon volume expands more slowly than the source volume;  (2) 
there is a suggestion of a trend for $n$
to increase with source size; and (3) in general $n<2$ and is closer to
9/5 than to 2.

%******************************************************************************

%******************************************************************************

\section{LOW POWER JETS: SONIC BOOMS, THE FR II-FR I TRANSITION, AND
RELICS}
\label{sec-jetpower}


One of the fundamental phenomenological results in the study of radio
galaxies is the differences in properties of the radio sources
in the powerful FRII and weaker FRI radio sources (e.g., Fanaroff
\& Riley 1974; Bridle 1984; O'Dea 1985; Laing 1993; Owen \& Ledlow
1994; Baum, Zirbel \& O'Dea 1995). A working paradigm is that the lower
luminosity plume-like FRIs have transonic non-relativistic jets
on large scales, while the higher luminosity edge-brightened
FRIIs have supersonic relativistic jets on large scales. However,
the reasons for these differences is still a matter of discussion
(e.g.,  De Young 1993; Bicknell 1994, 1995; Komissarov
1994; Baum \etal\ 1992,1995; Gopal-Krishna, Wiita and Hooda 1996;
Gopal-Krishna \& Wiita 1988, 2000, 2001; Meier 1999).

The FRI/FRII dichotomy has been addressed with numerical simulations.
Hardee, Clarke  \& Howell (1995) argued that in a three dimensional
simulation at Mach 5 the jet is disrupted on a scale length shorter
than observed in FR II type radio sources and, therefore, could not
possibly describe them but, instead, could describe low power FR I
sources.
The deceleration of initially relativistic jets due to entrainment
has been studied by Bowman \etal\ (1996).

Here we discuss the behavior of relatively low Mach number and low
power jets ($L_j \lapprox 3 \times 10^{45}$ erg.s$^{-1}$, Table 1).
Fig.'s \ref{ASUA_Vh_Xh} shows that at low Mach number, $v_h$ decreases
asymptotically to the sound speed of the external
medium.\footnote{The implications of this for the behavior of the source
size are discussed in \S \ref{sec-zh}.} For $M=2$, the speed of the head
becomes transonic by the time it has propagated just over 20 jet radii
($\sim 10$ kpc). In these low power sources, it is clear that, if
the bow shock does not further decelerate and its velocity is kept
constant, the source will expand in an apparently self-similar fashion
with $z_h\propto t$.

The results of the simulation of the low power jets in Run B1
demonstrate
that the jet does not keep up with the bow shock and essentially
stops feeding the source. This can also be seen for Run B2 and, in
Figs. \ref{ASUA_View2}, for Runs A6 and A7 although, here
the effect is less dramatic. The disruption of the jet may be due
to the many internal shocks in the jets. Thus, at these low Mach
numbers, the jet loses its ability to drive the shock at the working
surface and its speed is no longer supersonic with respect to the
ambient gas. The bow shock continues advancing but at a speed that
approaches the ambient sound speed, i.e.,  behaving as a sonic wave.
The pressure profile in Fig.  \ref{ASUA_Pres2} shows that the bow
shock has become a small perturbation. This can be seen more clearly
in Fig. \ref{ASUA_LowPower} where we have plotted the advance speed of
the jet head and the bow shock just in front of the jet as a function
of time for the low power jets (Runs B1, B2, A6 and A7). We see that
for $M=2$, the bow shock speed rapidly reaches the ambient sound
speed while the jet loses power and becomes subsonic.
Eventually, the jet goes through phases of negative velocity and
recedes. We also note that, by the end of the simulations, these
are relatively old sources with ages between 5 and 7 $\times
10^7$ yr and have low cocoon pressure. Our results are consistent with
those
of Reynolds, Heinz, \& Begelman (2001a) where they turn off the
jet in their simulations and find the bow shock becomes a ``sonic boom"
which continues to propagate through the ambient medium. Note that we
find that that this phenomenon will also occur in low power
jets which ``stall" rather than turning off altogether and thus,
this may be a somewhat more general phenomenon in the environments
of radio galaxies.

In general, at late stages in the low power jets, the pressure inside
the cocoon is very low compared with the ambient pressure (Fig.
\ref{ASUA_Pres2}). In addition,  at late times these sources are
not over-pressured near the head. Therefore, we do not expect
these sources to exhibit an edge-brightened radio image characteristic
of FR-II sources, although they may have appeared so at earlier
stages when the head pressure was higher. To illustrate this more
clearly, we have created a synthetic radio map of the low Mach number
Run A6 ($M=3.3$). We assume that only the jet and the cocoon, formed
by the shocked jet material, contain a population of emitting
relativistic electrons and magnetic field $B$
and that they are in equipartition. The synchrotron emissivity at a
given frequency is proportional to $\propto n_e B^{(b+1)/2}$,
where $n_e$ is the electron density and $b$ the electron energy
spectral index. We then take the electron number density and the
magnetic
energy proportional to the hydrodynamical pressure $P$ (Mioduszewski,
Hughes \& Duncan 1997) and write the emissivity as $\propto
P^{(b+5)/4}$.
The radio emissivity calculated in this way was attributed to each
grid cell according to the gas pressure and, after rotating the
simulation grid around the jet axis, we integrated over the entire
volume of
the source to obtain the radio intensity.

The results are shown Fig. \ref{RadioFR1_2G} where we present,
separately,
gray-scale and contour maps of the jet, cocoon and jet plus cocoon
for Run A6 at five stages of the source evolution. We can see how the
bright end of the jet in (a) and (b) associated with
the hot spot becomes less luminous in (c) and (d) until it disappears
in (e). At this stage, the jet itself has lost most of its power and
has retracted towards the nucleus (jet inlet). Initially, while the jet
powers it, the cocoon maintains a classical lobe configuration;
however, by the end of the simulation it has acquired more of a
plume-like
structure with diffuse ends. Thus, we  observe a change from an
edge-brightened morphology in (a) and (b) characteristic of FR-II
double sources, to an edge-darkened morphology in (d) and (e) as
if it had transitioned to a FR-I type source. This is consistent with
the scenario discussed by Gopal-Krishna \& Wiita (2001).
We suggest that these morphological changes will also occur in sources
when  the jet has turned off and the source is becoming a ``relic" (see
also Reynolds, Heinz, \& Begelman 2001b; Br\"uggen \& Kaiser 2000).

These results are suggestive and argue in favor of a scenario in
which (1) the difference in jet power (and/or Mach number) is
responsible
for the differences between the morphology of FR I and FR II sources
and (2) FR IIs could evolve into FR Is as the source ages and the
jet becomes disrupted and can no longer produce a hot spot.


%******************************************************************************

%******************************************************************************

\section{DEPENDENCIES ON THE RELATIVISTIC PARTICLE CONTENT OF THE JETS}
\label{case1}


We have considered two cases for the relativistic particle content of
the jets: Case 1
with relativistic electrons (and positrons) and Case 2 with relativistic
electrons and protons (see \S~\ref{nonrel}). The main difference between
Case 1 and 2 is that, for a given
$\eta$ and Lorentz factor $\gamma_j$, the equivalent classical jets in
Case 2 are
relatively heavier than in Case 1 and have lower Mach number. This can
be
understood as being due to the presence of the proton component. The
protons
increase the total enthalpy since their energy is much larger than that
of
the electrons ($k=100$).

Comparing Case 1 and 2 jets with identical equivalent Lorentz factor  we
see
that Case 1 jets: (1) propagate faster (Fig.  \ref{ASUA_View1} and Fig.
\ref{ASUA_View2}); (2) have cocoons with a smaller aspect ratio
(Fig.  \ref{ASUA_View1} and Fig. \ref{ASUA_View2});
(3) have higher pressure in the cocoons (Figs. \ref{ASUA_Pres1} and
\ref{ASUA_Pres2}), i.e., the pressure drops more slowly with distance
behind the source head (4) maintain higher lateral expansion velocity in
the bow shocks (Fig. ~\ref{ASUA_Vbc_x}). In Paper II
of this series, where we study jets propagating in an isothermal
atmosphere
with a density gradient, we find similar differences between Case 1 and
2 jets.

These differences generally become less important as the jet Lorentz
factor
increases and at the highest value we have studied ($\gamma_j = 10$),
the properties of the source with the heavy Case 2 jet (Run A10), do not
differ significantly from the analogous lighter
Case I jet  (Run A5). In addition, the degree to which the source
expansion is
self-similar is not very dependent on the particle content (e.g., Fig.
\ref{ASUA_Similar_Xh} and \ref{ASUA_Vh_Xh}).

We find that some aspects of the source energetics (discussed in Paper
II) are indeed different for Case 1 and 2. However, there is good
agreement in the evolution of the cocoon pressure with volume, the
ratio between the cocoon internal and total energy, and the ratio
between the head and cocoon pressure (see Fig.'s 15 and 19 of Paper II)
at large $\gamma_j$.

Thus, the material content of the relativistic jet does affect its
propagation and the overall structure and energetics (see Paper II)
of the source.

%******************************************************************************

%******************************************************************************

\section{CONSISTENCY AND ROBUSTNESS OF THE NUMERICAL RESULTS}
\label{RunC}

Here we comment on the results of the set of simulations (Run C) that
were performed to investigate the sensitivity of the numerical results
to small changes in the initial jet parameters. The input parameters
of Run C were chosen  to be only slightly different from
those of Runs B3 and B8 (see Table 1).


Comparing the density contour maps for Runs C3-C6 and Run B3 we
notice some variations from run to run in the cocoon internal
structure and, in a few runs, on its overall shape. This shows that
the turbulent behavior of the shocked jet material inside the cocoon
is quite sensitive to the initial conditions (which is expected).
However, the {\it
global\/} properties of the source do not show any substantial
differences. We have
plotted many quantities such as head radius, advance speed, volume,
pressure, and we find no significant variation in these and other
parameters such as the exponents of the scaling laws. For example, a
plot of the exponent $m$ in $z_h\propto t^m$ for Runs C1-C6,
B3 and B8 shows that, while there are minor variations, the general
behavior of the exponent $m$ is quite similar among the runs in a given
class. There is good agreement in the behavior of all the other
quantities we have examined. The main conclusion is that small
changes in either the Mach number $M$ or the density contrast $\eta$
do not result in major changes in the overall properties of
the source. Thus, we are confident that the variations in behavior
observed among the other runs are physically meaningful.

Next, we anticipate some results to be presented in Paper II for a
jet propagating in an isothermal atmosphere with a density gradient
to show that those  results are also robust against small changes
in input parameters. For example, a similar plot of the exponent $m$
 in $z_h\propto t^m$ for Runs C11-C16, B13 and B18 for jets propagating
in an isothermal atmosphere with a density gradient (to be discussed
in Paper II) shows that, for runs within a given environment, no
major differences are found in the behavior of the exponents
of the scaling laws.

%******************************************************************************

%******************************************************************************

\section{SUMMARY }
\label{sec-summary}

In this paper we have undertaken a broad parameter space survey of
two-dimensional
hydrodynamic simulations of powerful jets propagating into a constant
density atmosphere.
We examine the evolution of the global properties of the sources as a
function of source
size, for a range of Mach number, density contrast, and jet power. We
consider the effect
on the morphology and dynamics of whether relativistic protons are
present in the jet. We
also compare our results with the expectations of current analytical
self-similar models
(see Table 3). The main results can be summarized as follows:

The material content of the relativistic jet does effect its propagation
and the overall
structure and energetics (see Paper II) of the source. The jet in Case 1
(pressure
dominated by relativistic electrons) propagates faster, the cocoon has
an smaller aspect
ratio, the cocoon maintains higher pressure, and the bow shock expands
at higher lateral
velocity than Case 2 jets (pressure dominated by relativistic electron
and protons).
However, these differences generally become less important as the jet
Lorentz factor increases.

Small variations in input Mach number and density contrast have
significant effects on the
resulting morphology of the cocoon (e.g., due to the development of K-H
instabilities).
However, the derived global parameters are not sensitive to small
changes in input
parameters.

The cocoon boundary does not have a smooth regular shape as assumed in
the analytic
models, but is very irregular due to K-H instabilities and mixing with
the shocked gas
inside the bow shock.  In some runs, the region of the cocoon near the
inlet has lost its
definition and is completely mixed with the shocked gas.

We have confirmed the suggestion of Wilson \& Falle (1985) that the
location of the first
shock in the jet depends linearly on the jet radius and Mach number
($\ell \simeq 2r_jM
$). We show that the result is largely independent of density contrast
and that over a
wide range in Mach number the relation can be used to make a crude
estimate of the jet Mach number.

The pressure in the cocoon is highest at the jet head and declines with
distance away from
the head. At intermediate  Mach number ($M \sim 10$) the cocoons are in
pressure balance
with the ambient medium over much of the length of the source; while at
higher ($M>20$) or
lower ($M< 5$) Mach number, the cocoons are over pressured or under
pressured,
respectively. Away from the high pressure head, the cocoon and the inner
SAG region are
generally in pressure balance, though the outer SAG region can be at
higher pressure.

The bow shock lateral expansion decelerates quickly behind the jet head,
but the velocity
remains somewhat supersonic. Similarly, the cocoon lateral expansion
velocity drops
quickly behind the jet head and becomes subsonic in the post bow shocked
gas.  This causes
the bow shock and cocoon to diverge in shape with the cocoon becoming
quasi-cylindrical,
and the bow show becoming quasi-conical.

The strong time dependence of the cocoon and SAG region pressure and
lateral expansion
speed is expected to influence the properties of the emission line
nebula associated with
the radio source. The slow expansion of the cocoon and bow shock over
most of their length
will limit the maximum velocity to which ambient clouds can be
accelerated.

We review the properties of self-similar models and classify them into
three classes: (1)
Type I in which the jet advances with constant speed in a constant
density atmosphere; (2)
Type II in which the ambient density falls as $d^{-2}$ but, because the
jet has a non-zero
constant opening angle, the head advance speed is constant; and (3) Type
III which assumes
a power law density distribution with an arbitrary exponent $\delta$ and
which advances
with a speed which drops as $t^{-2/5}$ in a uniform atmosphere.

The radius of the jet head is a weak function of density ratio and a
stronger function of
Mach number, with higher Mach number jets tending to have smaller jet
heads. For low Mach
number ($M < 10$) the jet head radius tends to decrease slightly with
time, while for
$M\gae 10$, the jet head radius increases with time. Superimposed on
these trends are
large fluctuations in the head radius due to the well known phenomenon
of vortex shedding.

The advance speed of the head decreases slowly for low Mach number and
remains roughly
constant for high Mach number. However, once the head speed reaches the
sound speed of the
external medium, it does not decrease further. The advance speed of the
jet head depends
on the head radius in agreement with the assumption of balance between
the head pressure
and the ram pressure from the ambient medium.

The volume of the SAG region increases with time at a rate which is
close, but not
identical to the predictions of the Type I self-similar models. The bow
shock expansion
rate seems to increase with time for small Mach number but is relatively
more stable at
high Mach number. The cocoon volume is harder to define due to the K-H
instabilities at
the contact discontinuity, but it generally increases more slowly than
the SAG region
volume.

The behavior of the time evolution of the source size ($z_h \propto
t^m$) is in general
more complex than the predictions of the self-similar models. The time
dependence of the
source size exhibits up to  3 phases whose relative length depends on
the jet Mach number.
Initially, the jet will propagate at constant velocity, then it starts
to decelerate due
to expansion of the jet head -- during which $m=3/5$. Once the head
speed slows to the
external sound speed, it does not decelerate further and continues to
propagate at
constant velocity until the jet is finally disrupted.

In low Mach number jets, the dissipation of energy through internal
shocks is able to
disrupt the jets and turn off the supply of energy to the cocoon. At
this stage, the bow
shock propagates as a sound wave and the morphology of the inferred
radio emitting
material becomes faint and diffuse. Thus, we suggest that low Mach
number jets may be able
to transition between FR II and FR I type radio morphologies.


%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

\acknowledgements

JCC acknowledges the financial support of PRONEX/Finep and CNPq and
the hospitality of STScI where this work was carried out. CPO received
partial
support from NASA through grant GO-07855.01-96A from the Space Telescope
Science
Institute, which is operated by the Association of Universities for
Research in Astronomy, Inc., under NASA contract NAS 5-26555.
We thank Paul Wiita for helpful discussions and the referee for detailed
and helpful comments.

\begin{thebibliography}{}

\bibitem[]{} Alexander, P., 2000, MNRAS, 319, 8

\bibitem[]{} Alexander, P., \& Pooley, G. G. 1996, in Cygnus A - Study
of a Radio Galaxy, eds. C. L. Carilli \& D. E. Harris, 
(Cambridge: Cambridge University Press), 149

\bibitem[]{} Aloy, M. A., Ibanez, J. M$^a$., Marti, J. M$^a$., Gomez,
J.-L., M\"uller, E., 1999, \apjl, 523, 125L

\bibitem[]{} Anderson, M. C., Jones, T. W., Rudnick, L., Tregillis, I.
L.,
\& Kang, H., 1994, \apjl, 421, L31

\bibitem[]{} Baldwin, J. E., 1982, in Extragalactic Radio Sources,
Proceedings of IAU Symposium No. 97, ed. D. S. Heeschen \& C. M. Wade,
(Dordrecht: Reidel), 21

\bibitem[]{} Baum S. A., Heckman, T. M., \& van Breugel, W., 1992, ApJ,
389, 208

\bibitem[]{} Baum S. A., Zirbel, E. L., \& O'Dea C. P., 1995, ApJ, 451,
88

\bibitem[]{} Basset, G. M., \& Woodward, P. R. 1995, \apj, 441, 582

\bibitem[]{} Begelman, M. C. 1996, in Cygnus A: Study of a Radio Galaxy,
ed.  C. L. Carilli \& D. A. Harris (Cambridge: Cambridge University Press),
209

\bibitem[]{} Begelman, M. C., \& Cioffi, D. F. 1989, ApJ, 345, L21

\bibitem[]{} Best, P. N., R\"ottgering, H. J. A., \& Longair, M. S.,
2000, \mnras,  311, 23

\bibitem[]{} Bicknell, G. V., 1994, ApJ, 422, 542

\bibitem[]{} Bicknell, G. V., 1995, ApJS, 101, 29

\bibitem[]{} Bicknell, G. V., Dopita, M. A., \& O'Dea, C. P., 1997, ApJ,
485,
112

\bibitem[]{} Bicknell, G. V., \& Begelman, M. C. 1996, ApJ, 467, 597

\bibitem[]{} Blandford, R. D. 1994, in The Physics of Active Galaxies,
ed. G. V. Bicknell, M. A. Dopita, \& P. J. Quinn, (San Francisco: ASP), 23

\bibitem[]{} Blundell, K. M., Rawlings, S., \& Willott, C. J., 1999,
\aj, 117, 677

\bibitem[]{} Bodo, G., Massaglia, S., Ferrari, A., \& Trussoni, E. 1994,
A\&A, 283, 655

\bibitem[]{} Bodo, G., Massaglia, S., Rossi, P., Rosner, R., Malagoli,
A., \& Ferrari, A. 1995, A\&A, 303, 281

\bibitem[]{} Bodo, G., Rossi, P., Massaglia, S., Ferrari, A., Malagoli,
A., \& Rosner, R. 1998, A\&A, 333, 1117

\bibitem[]{} Bowman, M., Leahy, J. P., \& Komissarov, S. S. 1996,
\mnras, 279, 899

\bibitem[]{} Bridle, A. H., 1984, AJ, 89, 979

\bibitem[]{} Bridle, A. H., \& Perley, R. A., 1984, \araa, 22, 319

\bibitem[]{} Br\"uggen, M., \& Kaiser, C. R., 2000, \mnras, in press

\bibitem[]{} Carilli, C. L., Perley, R. A., Bartel, N., \& Sorathia, B.
1996, in ASP Conf.  Ser. 100, Energy Transport in Radio Galaxies and
Quasars, ed. P. E. Hardee, A. H. Bridle, \& J. A. Zensus, 
(San Francisco: ASP), 287

\bibitem[]{} Carvalho, J. C., 1985, \mnras, 215, 463

\bibitem[]{} Carvalho, J. C., \& O'Dea, C. P., 2001a, ApJS, in press
(Paper II)

\bibitem[]{} Carvalho, J. C., \& O'Dea, C. P., 2001b,  in preparation
(Paper III)

\bibitem[]{} Cioffi, D. F., \& Blondin, J. M. 1992, ApJ, 392, 458

\bibitem[]{} Clarke, D. A., Norman, M. L., \& Burns, J. O. 1989, \apj,
342, 700

\bibitem[]{} Clarke, D. A., \& Burns, J. O. 1991, \apj, 369, 308

\bibitem[]{} Clarke, D. A. 
1996, in ASP Conf.  Ser. 100, Energy Transport in Radio Galaxies and
Quasars, ed. P. E. Hardee, A. H. Bridle, \& J. A. Zensus, 
(San Francisco: ASP), 311

\bibitem[]{} Clarke, D. A., Harris, D. E., \& Carilli, C. L. 1997,
\mnras, 284, 981

\bibitem[]{} Cox, C. I., Gull, S. F., \& Scheuer, P. A. G. 1991, \mnras,
252, 558

\bibitem[]{} Daly, R. A. 1990, ApJ, 355, 416

\bibitem[]{} Daly, R. A., 1994, \apj, 426, 38

\bibitem[]{} Daly, R. A., 2001, in Life Cycles of Radio Galaxies,
ed. J. Biretta, A. Koekemoer, C. O'Dea, \& E. Perlman (New Astronomy
Reviews), in press.

\bibitem[]{} De Vries, W. H., O'Dea, C. P., Baum, S. A., \& Barthel, P.
D., 1999, ApJ, 526, 27

\bibitem[]{}  De Vries, W. H., O'Dea, C. P., Barthel, P. D.,
Fanti, C., Fanti, R., \& Lehnert, M. D., 2000, \aj, 120, 2300


\bibitem[]{} De Young, D. S., 1993, \apj, 405, L13

\bibitem[]{} Duncan, G. C., \& Hughes, P. A. 1994, ApJ, 436, L119

\bibitem[]{} Falle, S. A. E. G. 1991, MNRAS, 250, 581

\bibitem[]{} Fanti, C., Fanti, R., Dallacasa, D., Schilizzi, R. T.,
Spencer, R. E., \& Stanghellini, C., 1995, A\&A, 302, 317

\bibitem[]{} Ferrari, A., Massaglia, S., \& Bodo, C. 1996, in
Extragalactic Radio Sources, IAU Symposium No. 175, ed, 
R. D. Ekers, C. Fanti, L.  Padrielli (Dordrecht: Kluwer Academic 
Publishers), 449

\bibitem[]{} Frank, A., Ryu, D., Jones, T. W., \&
Noriega-Crespo, A. 1998, ApJ, 494, L79.

\bibitem[]{} Garrington, S. T., Leahy, J. P., Conway, R. G., \&
Laing, R. A., 1988, Nature, 133, 147

\bibitem[]{} Gomez, J. L., Marti, J. M$^a$., Marscher, A. P., Ibanez, J.

M$^a$., \& Marcaide, J. M. 1995, \apjl, 449, L19

\bibitem[]{} Gopal-Krishna, \& Wiita, P. J., 1988, Nature, 333, 49

\bibitem[]{} Gopal-Krishna, Wiita, P. J., \& Hooda, J. A., 1996, A\&A,
316, L13

\bibitem[]{} Gopal-Krishna, \& Wiita, P. J., 2000, A\&A, 363, 507

\bibitem[]{} Gopal-Krishna, \& Wiita, P. J., 2001, A\&A, 373, 100

\bibitem[]{} Guerra, E. J., Daly, R. A., \& Wan, L., 2000, \apj, 544,
659

\bibitem[]{} Hardee, P. E., \& Norman, M. L. 1988, ApJ, 334, 70

\bibitem[]{} Hardee, P. E., \& Norman, M. L. 1990, ApJ, 365, 134

\bibitem[]{} Hardee, P. E., White III, R. R., Norman, M. L., Cooper, M.
A., \&
Clarke, D. A. 1992, ApJ, 387, 460

\bibitem[]{} Hardee, P. E., Clarke, D. A., \& Howell, D. A. 1995, \apj,
441, 644

\bibitem[]{} Heckman, T. M., van Breugel, W. J. M., \& Miley, G. K.
1984,
ApJ, 286, 509

\bibitem[]{} Hooda, J. S. \& Wiita, P. J. 1996, ApJ, 470, 211

\bibitem[]{} Kaiser, C. R., \& Alexander, P. 1997, \mnras, 286, 215

\bibitem[]{} Kaiser, C. R., Dennett-Thorpe, J., \& Alexander, P., 1997,
MNRAS, 292, 723

\bibitem[]{} Klein, R. I., McKee, C. F., \& Colella, P., 1994, ApJ, 420,
213

\bibitem[]{} Koekemoer, A. M., \& Bicknell, G. V., 1998, ApJ, 497, 662

\bibitem[]{} K\"ossl, D., \& M\"uller, E. 1988, A\&A, 206, 204

\bibitem[]{} K\"ossl, D., M\"uller, E., \& Hillebrandt, W. 1990, A\&A,
229, 397

\bibitem[]{} Komissarov, S. S., 1994, MNRAS, 269, 394

\bibitem[]{} Komissarov, S. S., \& Falle, S. A. E. G. 1997, MNRAS, 288,
833

\bibitem[]{} Komissarov, S. S., \& Falle, S. A. E. G. 1998, MNRAS, 297,
1087

\bibitem[]{} Laing, R. A., 1988, Nature, 331, 149

\bibitem[]{} Laing, R. A., 1993, in Astrophysical Jets, 
ed. D. Burgarella, M. Livio, \& C. O'Dea,
(Cambridge: Cambridge University Press), 95


\bibitem[]{} Lind, K. R., Payne, D. G., Meier, D. L., \& Blandford R. D.

1989, ApJ, 344, 89

\bibitem[]{} Loken, C., Burns, J. O., Clarke, D. A., \& Norman, M. L.
1992, ApJ,
392, 54

\bibitem[]{} Marti, J. M$^a$., M\"uller, E., \& Ibanez, J. M$^a$. 1994,
A\&A, 281, L9

\bibitem[]{} Marti, J. M$^a$., M\"uller, E., Font, J. A., Ibanez, J.
M$^a$., \&
Marquina, A. 1997, ApJ, 479, 151

\bibitem[]{} McCarthy, P. J., 1993, \araa, 31, 639

\bibitem[]{} Meier, D. L., 1999, \apj, 522, 753

\bibitem[]{} Mioduszewski, A. J., Hughes, P. A., \& Duncan, G. C. 1997,
ApJ,
476, 649

\bibitem[]{} Nath, B. B. 1995, MNRAS, 274, 208

\bibitem[]{} Neeser, M. J., Meisenheimer, K., \& Hippelein, H., 1997,
ApJ, 491,
522

\bibitem[]{} Norman, M.L., Smarr, L., Winkler, K.-H.A., \& Smith, M.D.
1982,
A\&A, 113, 285

\bibitem[]{} Norman, M.L., \& Hardee, P.E. 1988, ApJ, 334, 80

\bibitem[]{} Norman, M.L. 
1996, in ASP Conf.  Ser. 100, Energy Transport in Radio Galaxies and
Quasars, ed. P. E. Hardee, A. H. Bridle, \& J. A. Zensus, 
(San Francisco: ASP), 319

\bibitem[]{} O'Dea C.P., 1985, ApJ, 295, 80

\bibitem[]{} O'Dea C.P., \& Baum S.A. 1997, AJ, 113, 148

\bibitem[]{} O'Dea C.P., 1998, PASP, 110, 493

\bibitem[]{} O'Dea C.P., \etal\ 2001, in preparation %CSS paper

\bibitem[]{} Owen, F. N., \& Ledlow, M. J., 1994,
in Proceedings of ``The First
Stromlo Symposium: The Physics of Active Galaxies," ASP Conf. Ser. 54, 
ed. G. V. Bicknell, M. A. Dopita, \& P. J.  Quinn, 
(San Francisco: ASP), 319

\bibitem[]{} Perley, R.A., Dreher, J.W., \& Cowan, J.J. 1984, ApJ, 285,
L35

\bibitem[]{} Rayburn, D. R. 1977, MNRAS, 179, 603

\bibitem[]{} Readhead, A. C. S., Taylor, G. B., Pearson, T. J., \&
Wilkinson, P. N., 1996, ApJ, 460, 634

\bibitem[]{} Reynolds, C. S., Heinz, S., \& Begelman, M. C., 2001a,
\apjl, 549,
179L

\bibitem[]{} Reynolds, C. S., Heinz, S., \& Begelman, M. C., 2001b,
\mnras, in press

\bibitem[]{} Rosen, A., Hughes, P. A., Duncan, G. C., \& Hardee, P. E.
1999,
ApJ, 516, 729

\bibitem[]{} Scheuer P.A.G. 1974, MNRAS, 166, 513

\bibitem[]{} Scheuer P.A.G. 1982, in Extragalactic Radio Sources,
Proceedings of IAU Symposium No. 97, ed. D. S.
Heeschen \& C. M. Wade, (Dordrecht: Reidel), 163

\bibitem[]{} Smith, M. D., Norman, M. L., Winkler, K.-H. A., \& Smarr,
L. 1985, \mnras, 214,67

\bibitem[]{} Spinrad, H., Dey, A., \& Graham, J. R., 1995, \apj, 438,
L51

\bibitem[]{} Stern, D., \& Spinrad, H., 1999, PASP, 111, 1475

\bibitem[]{} Stone, J. M., \& Norman, M. L. 1992a, ApJS, 80, 753

\bibitem[]{} Stone, J. M., \& Norman, M. L. 1992b, ApJS, 80, 791

\bibitem[]{} Stone, J. M., Mihalas, D., \& Norman, M. L. 1992, ApJS, 80,
819

\bibitem[]{} Tadhunter, C. N., Villar-Mart\'in, M., Morganti, R.,
Bland-Hawthorn,
J., \& Axon, D., 2000, MNRAS, 314, 849

\bibitem[]{} Taylor, G. B., Barton, E. J., \& Ge, J., 1994, \aj, 107,
1942

\bibitem[]{} Tregillis, I. L., Jones, T. W., \& Ryu, D. 2001,
astro-ph/0104305

\bibitem[]{} Villar-Mart\'in, M., Tadhunter, C., Morganti, R., Axon, D.,
\&
Koekemoer, A., 1999, \mnras, 307, 34

\bibitem[]{} Williams, A. G. 1991, in Beams and Jets in Astrophysics,
ed. P. A. Hughes, (Cambridge: Cambridge Univ.  Press),  342

\bibitem[]{} Wilson, M. J., \& Falle, S. A. E. G., 1985, MNRAS, 216, 971

\bibitem[]{} Wilson, M. J. 1987, MNRAS, 226, 447

\bibitem[]{} Wang, Z., Wiita, P. J., \& Hooda, J. S. 2000, \apj, 534,
201

\bibitem[]{} Zensus, J. A., 1997, \araa, 35, 607

\end{thebibliography}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%******************************************************************************

%******************************************************************************

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\clearpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{deluxetable}{rrrrrrr}
\tablecolumns{7}
\tablewidth{0pc}
\tablecaption{Simulation Parameters}
\tablehead{
\colhead{Run} & \colhead{Case}   & \colhead{$\gamma_j$} &
\colhead{$\beta_j$} &
\colhead{$\eta_c$}    & \colhead{$M_c$}   & \colhead{Jet
Power\tablenotemark{a}}}
\startdata
A1  &    1 &    2.0    &    0.866  &   0.0044  &   9.7  & 1.56 (47)\\
A2  &    1 &    2.5    &    0.917  &   0.0069  &  12.8  & 2.86 (47)\\
A3  &    1 &    5.0    &    0.980  &   0.0276  &  27.4  & 1.40 (48)\\
A4  &    1 &    7.0    &    0.990  &   0.0542  &  38.8  & 2.84 (48)\\
A5  &    1 &    10.0   &    0.995  &   0.1106  &  55.7  & 5.88 (48)\\
A6  &    2 &    2.0    &    0.866  &   0.0210  &  3.3   & 2.80 (45)\\
A7  &    2 &    2.5    &    0.917  &   0.0328  &  4.4   & 5.32 (45)\\
A8  &    2 &    5.0    &    0.980  &   0.1310  &  9.4   & 2.60 (46)\\
A9  &    2 &    7.0    &    0.990  &   0.2568  & 13.3   & 5.25 (46)\\
A10 &    2 &    10.0   &    0.995  &   0.5241  & 19.2   & 1.11 (47)\\
\cline{1-7}
B1  &\nodata&  \nodata  &   \nodata &   0.100   &  2.0   & 2.86 (44)\\
B2  &\nodata&  \nodata  &   \nodata &   0.100   &  4.0   & 2.29 (45)\\
B3  &\nodata&  \nodata  &   \nodata &   0.100   &  6.0   & 7.73 (45)\\
B4  &\nodata&  \nodata  &   \nodata &   0.100   & 12.0   & 6.18 (46)\\
B5  &\nodata&  \nodata  &   \nodata &   0.100   & 20.0   & 2.86 (47)\\
B6  &\nodata&  \nodata  &   \nodata &   0.003   &  6.0   & 4.46 (46)\\
B7  &\nodata&  \nodata  &   \nodata &   0.010   &  6.0   & 2.44 (46)\\
B8  &\nodata&  \nodata  &   \nodata &   0.030   &  6.0   & 1.41 (46)\\
B9  &\nodata&  \nodata  &   \nodata &   0.100   &  6.0   & 7.73 (45)\\
B10 &\nodata&  \nodata  &   \nodata &   0.300   &  6.0   & 4.46 (45)\\
\cline{1-7}
C1  &\nodata&  \nodata  &   \nodata &   0.029   &  6.0   & 1.43 (46)\\
C2  &\nodata&  \nodata  &   \nodata &   0.035   &  6.0   & 1.31 (46)\\
C3  &\nodata&  \nodata  &   \nodata &   0.110   &  6.0   & 7.37 (45)\\
C4  &\nodata&  \nodata  &   \nodata &   0.100   &  6.1   & 8.12 (45)\\
C5  &\nodata&  \nodata  &   \nodata &   0.100   &  6.3   & 8.94 (45)\\
C6  &\nodata&  \nodata  &   \nodata &   0.100   &  5.9   & 7.35 (45)\\
\tablenotetext{a}{$erg.s^{-1}$}
%%%
\tablecomments{In Case 1 the jet relativistic particle content is
a pair plasma  - electrons (and positrons) - and in Case 2,
is baryon dominated - electrons plus protons.}
%%%
\enddata
\end{deluxetable}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\clearpage


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
\begin{deluxetable}{cccccc}
\tablecolumns{5}
\tablewidth{0pc}
\tablecaption{Types of Self-similar Models}
\tablehead{
%%\colhead{}    &  \multicolumn{3}{c}{Non-shell Stars} & \colhead{} &
%%\multicolumn{3}{c}{Shell Stars} \\
%%\cline{2-4} \cline{6-8} \\
\colhead{ } & \colhead{ }   & \colhead{Type I} & \colhead{Type II} &
\colhead{Type III}}
\startdata
density         &$\rho$ & {\it constant} & $d^{-2}$       &
$d^{-\delta}$\\
head speed      &$v_h$  & {\it constant} & {\it constant} &
$t^{\delta-2\over
5-\delta}$\\
head area       &$A_h$  & {\it constant} & $t^2$          &
$t^{\delta+4\over
5-\delta}$ \\
source size     &$z_h$  &      $t$       & $t$            & $t^{3\over
5-\delta}$
     \\
head pressure   &$P_h$  & {\it constant} & $t^{-2}$       &
$t^{-{\delta+4\over
5-\delta}}$\\
cocoon radius   &$r_c$  &     $t^{1/2}$  & $t^{1/2}$      & $t^{3\over
5-\delta}$
       \\
cocoon pressure &$P_c$  &      $t^{-1}$  & $t^{-1}$       &
$t^{-{\delta+4\over
5-\delta}}$  \\
            &$P_h/P_c$  &      $t$       & $t^{-1}$       & {\it
constant}
         \\
\tablecomments{Type I: Begelman \& Cioffi (1989); Loken et al. (1992);
Cioffi
\& Blondin (1992); Nath (1995) \\
Type II : Daly (1990) \\ Type III: Falle (1991); Begelman (1996);
Bicknell, Dopita
\&
O'Dea (1997); Kaiser \& Alexander (1997)}
\enddata
\end{deluxetable}


\clearpage

\begin{deluxetable}{cccccccccccc}
\tabletypesize{\small}
\rotate
\tablecolumns{11}
\tablewidth{0pc}
\tablecaption{Comparison of the Uniform Atmosphere Simulations with
Self-Similar
Models }
\tablehead{
\colhead{} &
\colhead{} & \multicolumn{2}{c}{Self-Similar Model} &
\colhead{} & \multicolumn{3}{c}{Simulation}
           & \multicolumn{3}{c}{Simulation} &
\colhead{} \\
\colhead{} & \colhead{} & \colhead{} & \colhead{} & \colhead{} &
 \multicolumn{3}{c}{case 1} & \multicolumn{3}{c}{case 2} & \colhead{}\\
%%\multicolumn{3}{c}{} \\
%\cline{2-3} \cline{5-7}\\
\colhead{Relationship} & \colhead{exponent} & \colhead{Type I} &
\colhead{Type III}
 & \colhead{}
 & \multicolumn{2}{c}{exponent}  & \colhead{note}
 & \multicolumn{2}{c}{exponent}  & \colhead{note} & }
\startdata
$z_h\propto t^m$ &m& 1 & 3/5 &&$\sim 0.6$& low $\gamma$ &  a,c  &$\sim
0.8$& low
$\gamma$&b,c &
Sec. 5.1 \\
                 & &   &     &&$\sim 1.0$&high $\gamma$ &  a,c  &$\sim
1.0$&high
$\gamma$&b,c &
        \\
\cline{1-12}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% bow shock
&&&&&&&&&&&\\ $V\propto t^n$   &n& 2 & 9/5 &SAG region
&2.0                &high
$\gamma$ & d,f &$>2$
  &high $\gamma$&e,f& Sec. 5.4\\

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% cocoon
                 & &   &     &cocoon   &$\sim1.0$ to $2.0$ &all
$\gamma$  & g,i
&$\sim1.8$ to
$2.0$&low  $\gamma$&h,i&         \\
                 & &   &     &         &
&              &
&$<1.8$
  &high $\gamma$&h,i&         \\

\tablecomments{
a - $m$ decreases with increasing source size and asymptotes to final
value\\
b - $m$ is relatively constant with source size\\
c - $m$ is in the interval $0.8-1.0$ for low $M$, $0.6-0.8$ for moderate
$M$ and
$\sim
        1.0$ for high $M$\\
d - $n$ is relatively constant with source size\\
e - $n$ increases moderately with increasing source size\\
f - bow shock runs off the grid at low $\gamma$\\
g - $n$ increases with increasing source size for high $\gamma$\\
h - $n$ increases with increasing source size\\
i - $n$ shows large oscillations as a function of source size}

\enddata
\end{deluxetable}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%   \end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[ht]
\figurenum{1}
\epsscale{0.7}
\plotone{f1.eps}
\figcaption{Position in the $\eta-M$ plane of the jets in Run A (solid
circles)
together with the additional sets: Run B (open circles) and Run C
(crosses).
The curves are contours of equal $k$ and equivalent Lorentz factor
$\gamma_j$.
The parameters of the runs are given in Table 1.
\label{ASU_EtaMach}
} \end{figure} \clearpage



\begin{figure}[ht]
\figurenum{2}
\epsscale{0.9}
\plotone{f2.eps}
\figcaption{Contours of density of the simulation of a light jet in Case
1
(Runs A1-A5) with the following values of $\eta_c$ and $M_c$: (a)
0.0044, 9.7;
(b) 0.0069, 12.8; (c) 0.0276, 27.4; d) 0.0542, 38.8 and (e) 0.1106,
55.7.
\label{ASUA_All1}
} \end{figure} \clearpage



\begin{figure}[ht]
\figurenum{3}
\epsscale{0.5}
\plotone{f3.eps}
\figcaption{
Time evolution of the jet simulations for Case 1 (Runs A1-A4). The
curves represent
the
outline of the bow shock, cocoon and jet. The jet Mach number and
density ratio are
indicated and the time in each panel is given in Myr.
\label{ASUA_View1}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{4}
\plotone{f4.eps}
\figcaption{
As in Fig. \ref{ASUA_View1} but for Case 2 (Runs A6-A9).
\label{ASUA_View2}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{5}
\epsscale{0.7}
\plotone{f5.eps}
\figcaption{
Position $z_s$ of the first shock inside the jet (vertical axis) plotted
against the
jet
response scale $\ell$ given by Eq. (\ref{ele}) (horizontal axis). Panels
(a) and (c)
correspond to Runs B and A, respectively, where we have adopted  the jet
inlet
radius (500
pc) as the jet radius in the expression of $\ell$. Panels (b) and (d)
correspond to
Runs B
and A, respectively, where the jet response scale has been calculated
using the
actual
average radius for the simulated jet.
\label{ASU_Shock}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{6}
\epsscale{0.8}
\plotone{f6.eps}
\figcaption{
Pressure along the direction of the jet axis (Case 1, Runs A1-A4)
inside the cocoon (thick line), in the outer SAG region $R_b/2 < R <R_b$
(thin line) and in the inner SAG region $R_c<R<R_b/2$ (dashed line).
The horizontal axis is the distance from the galaxy nucleus
(jet inlet) and the vertical axis is the logarithm of $P$ relative to
the ambient
pressure
$P_a$. The range on the vertical axis is from $P/P_a = 0.1$ to $P/P_a =
10^3$. The
horizontal dotted line represents $P/P_a=1$ and the time is given in
Myr.
\label{ASUA_Pres1}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{7}
\plotone{f7.eps}
\figcaption{
As in Fig. \ref{ASUA_Pres1} but for Case 2 (Runs A6-A9).
\label{ASUA_Pres2}
} \end{figure} \clearpage




\begin{figure}[ht]
\figurenum{8}
\epsscale{0.8}
\plotone{f8.eps}
\figcaption{
The ratio of the lateral expansion speed of the bow shock (thin line)
and cocoon boundary
(thick line) to the ambient sound speed (dotted line)
-- i.e., the sound speed of the medium the bow shock or cocoon is
expanding
into - for the cocoon this is the post bow shock gas). For each
simulation, the
results
are shown at the time step just before the bow shock leaves the lateral
boundary of the
grid, which is not necessarily the end of the simulation.
 \label{ASUA_Vbc_x}
} \end{figure} \clearpage

\begin{figure}[ht]
\figurenum{9}
\plotone{f9.eps}
\figcaption{
Source size ($z_h$) as a function of time (open circles) for Runs
A1-A10. The solid
line
is the least squares  fit to the data, and the dotted and dashed lines
 indicate $z_h\propto t$ (Type I) and $z_h\propto t^{3/5}$ (Type III),
respectively.
\label{ASUA_Xh_Tt}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{10}
\plotone{f10.eps}
\figcaption{
Exponent $m$ of the relationship $z_h\propto t^m$ estimated using
the two methods described in the text (\S~\ref{sec-zh}) for Runs A1-A10.
The open circles are the values derived using the first method, and
the horizontal axis represents the position of the first fitted point.
The solid line is the logarithmic derivative of $z_h$
with respect to $t$ (derived using the second method) and plotted as
a function of $z_h$.
The dotted lines indicate the values predicted by the Type I (1) and
Type III (3/5) models.
\label{ASUA_Similar_Xh}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{11}
\plotone{f11.eps}
\figcaption{
Average radius of the jet (solid circles) and cocoon (open circles) at
the source
head as
a function of the size of the source. The solid line is described in
section
\ref{sec-Vh}.
\label{ASUA_Rh_Xh}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{12}
\epsscale{0.5}
\plotone{f12.eps}
\figcaption{
(a) Average radius of the jet head (triangles) and cocoon (open circles)
and (b) the
rate of increase in the size of the jet (triangles) and cocoon (open
circles) head radius
with source size (growth rate) as a function of the Mach number for all
26 runs (A, B and
C). The solid circles represent the head radius calculated using the
``observed" head speed from Eq. (\ref{vh}).
\label{ASU_Rmed_ABC}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{13}
\epsscale{0.8}
\plotone{f13.eps}
\figcaption{
Advance speed of the source head $v_h$ as a function of source size for
Runs A1-A10 (open circles). The solid and dotted lines are the advance
speed calculated using Eq.  (\ref{vh}) and adopting the derived radius
of the cocoon head and jet head, respectively, given in
Fig. \ref{ASUA_Rh_Xh}. We show the sound speed in the ambient medium
(dashed lines) except where it is ($< 10^{-2} v_j$).
\label{ASUA_Vh_Xh}
} \end{figure} \clearpage



\begin{figure}[ht]
\figurenum{14}
\plotone{f14.eps}
\figcaption{
Volume of the region limited by the bow shock ($V_b$) as a function of
time. The solid line
represents the least squares fit to the simulation results where we have
avoided the
transient regime $x\lapprox 10$ kpc. The dotted and dashed lines denote
$V\propto t^{2}$
(Type I) and $V\propto t^{9/5}$ (Type III) respectively.
\label{ASUA_Vb_T}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{15}
\plotone{f15.eps}
\figcaption{
Volume of the cocoon as a function of time. The solid line represents
the least squares
fit to the simulation results where we have avoided the transient
regime  $x\lapprox 10$
kpc. The dotted and dashed lines denote $V\propto t^{2}$ (Type I) and
$V\propto t^{9/5}$ (Type III) respectively.
\label{ASUA_Vc_T}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{16}
\plotone{f16.eps}
\figcaption{
Exponent $n$ of $V_b\propto t^n$ calculated according the two methods
described in
the text (\S~\ref{sec-zh}) for Runs A1-A10.
The open circles are the values derived using the first method, and
the horizontal axis represents the position of the first fitted point.
The solid line is the logarithmic derivative $d(\log{V_b})/d(\log{t})$
(derived using the second method) and plotted as a function of $z_h$.
The dotted lines are the values predicted by Type I (2) and Type III
(9/5) models.
\label{ASUA_Similar_Vb}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{17}
\plotone{f17.eps}
\figcaption{
As in Fig. \ref{ASUA_Similar_Vb} but for the cocoon volume.
\label{ASUA_Similar_Vc}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{18}
\epsscale{0.6}
\plotone{f18.eps}
\figcaption{
Advance speed of the bow shock (open circles) and of the jet head
(crosses) relative to
the jet initial speed as a function of time for the four lowest power
jets we have
simulated (a) Run A6, (b) Run A7, (c) Run B1 and (d) Run B2. The
horizontal line denotes the sound speed of the ambient medium.
\label{ASUA_LowPower}
} \end{figure} \clearpage


\begin{figure}[ht]
\figurenum{19}
\epsscale{0.3}
\plotone{f19.eps}
\figcaption{
Contour and gray-scale synthetic radio images of Run A6. The upper panel
shows the radio
emission of the jet, the central panel the cocoon and the bottom panel
the sum of the jet and cocoon. Each panel shows the same five 
stages of evolution which are shown in Fig.
\ref{ASUA_View2}. The darkest area corresponds to the maximum intensity
and the same intensity scale is used in each panel.
\label{RadioFR1_2G}
} \end{figure} \clearpage



\end{document}

