A&A 427, 415-429 (2004)
DOI: 10.1051/0004-6361:20040349

Stability of hydrodynamical relativistic planar jets

I. Linear evolution and saturation of Kelvin-Helmholtz modes

M. Perucho¹ - M. Hanasz² - J. M. Martí¹ - H. Sol³

1 - Departamento de Astronomía y Astrofísica, Universidad de Valencia, 46100 Burjassot, Spain
2 - Torun Centre for Astronomy, Nicholas Copernicus University, 97-148 Piwnice k. Torunia, Poland
3 - LUTH, Observatoire de Paris-Meudon, Pl. J. Jansen 5, Meudon 92195, France

Received 27 February 2004 / Accepted 11 July 2004

Abstract
The effects of relativistic dynamics and thermodynamics in the development of Kelvin-Helmholtz instabilities in planar, relativistic jets along the early phases (namely linear and saturation phases) of evolution has been studied by a combination of linear stability analysis and high-resolution numerical simulations for the most unstable first reflection modes in the temporal approach. Three different values of the jet Lorentz factor (5, 10 and 20) and a few different values of specific internal energy of the jet matter (from 0.08 to 60.0 c²) have been considered. Figures illustrating the evolution of the perturbations are also shown.

Our simulations reproduce the linear regime of evolution of the excited eigenmodes of the different models with a high accuracy. In all the cases the longitudinal velocity perturbation is the first quantity that departs from the linear growth when it reaches a value close to the speed of light in the jet reference frame. The saturation phase extends from the end of the linear phase up to the saturation of the transversal velocity perturbation (at approximately 0.5c in the jet reference frame). The saturation times for the different numerical models are explained from elementary considerations, i.e. from properties of linear modes provided by the linear stability analysis and from the limitation of the transversal perturbation velocity. The limitation of the components of the velocity perturbation at the end of the linear and saturation phases allows us to conclude that the relativistic nature of the flow appears to be responsible for the departure of the system from linear evolution.

The high accuracy of our simulations in describing the early stages of evolution of the KH instability (as derived from the agreement between the computed and expected linear growth rates and the consistency of the saturation times) establishes a solid basis to study the fully nonlinear regime, to be done elsewhere. The present paper also sets the theoretical and numerical background for these further studies.

Key words: galaxies: jets - hydrodynamics

1 Introduction

Relativistic jets are found in many astrophysical objects like active galactic nuclei (see e.g. Zensus 1997), microquasars (Mirabel & Rodriguez 1999) and gamma-ray bursts (Sari et al. 1999). The basic physical property of jets is their huge flux of kinetic energy, taken away from compact unresolved or marginally resolved central regions. The presence of jets in such objects provides a unique opportunity to study the physics of these central regions through investigations of the jet flows. However, the extraction of information about the central object is not a simple task, because jets interact in complex ways with the surrounding interstellar and intergalactic medium. On one hand the dynamical interaction of the jet matter with the ambient medium leads to the formation of shocks, turbulence, acceleration of charged particles and subsequent emission of a broad-range electromagnetic radiation, which makes jets observable. On the other hand the complex nature of the interaction which is due to the presence of flow instabilities, especially the Kelvin-Helmholtz (KH) instability, makes it difficult to distinguish those properties which are directly related to the central object from those which are due to the interaction of jets with the ambient medium. As an example one can invoke the wavelike patterns in jets, which may result from the precession of the rotation axis of the accretion disk at the jet base, or from KH instability, which finds favorable conditions at the interface of jet and external medium (Trussoni et al. 1983). Very recently, the KH instability theory has been successfully used to interpret the structure of the pc-scale jet in the radio source 3C 273 (Lobanov & Zensus 2001).

At kiloparsec scales, the surprisingly stable propagation of relativistic jets in some sources (e.g., Cyg A) contrasts with the deceleration and decollimation observed in other sources (e.g., 3C 31). These two sources are representative examples of two kinds, which are classified into Fanaroff-Riley type I (FRI) and Fanaroff-Riley type II (FRII) radio sources (Fanaroff & Riley 1974). The morphological dichotomy of FRI and FRII sources may be related to the stability properties of relativistic jets with different kinetic powers.

This complex situation motivated us to study the interaction of relativistic jets with their ambient medium and more specifically the KH instability in detail, by applying a linear stability analysis along with numerical simulations.

The linear analysis of KH instability in relativistic jets started with the work of Turland & Scheuer (1976) and Blandford & Pringle (1976) who derived and solved a dispersion relation for a single plane boundary between the relativistic flow and the ambient medium. Next, Ferrari et al. (1978) and Hardee (1979) examined properties of KH instability in relativistic cylindrical jets by following the derivation of the dispersion relation in the nonrelativistic case, in the vortex sheet approximation, done by Gill (1965). They numerically solved the dispersion relation, found unstable KH modes and classified them into the fundamental (surface) and reflection (body) family of modes. The classification is related to the number of nodes, across the jet, of sound waves reflecting in between jet boundaries. An internal wave pattern is formed by the composition of oblique waves, for which the jet interior is a resonant cavity. The physical meaning of KH instability in supersonic jets has been discussed by Payne & Cohn (1985), who have shown that the presence of instability is associated with the overreflection of soundwaves (the modulus of reflection coefficient is larger than 1) on the sheared jet boundaries.

Subsequent studies include effects of magnetic field (Ferrari et al. 1981), the effects of the shear layer, replacing the vortex sheet in the nonrelativistic planar case (Ferrari et al. 1981; Ray 1982; Choudhury & Lovelace 1984), and conical and cylindrical jet geometry (Birkinshaw 1984; Hardee 1984, 1986, 1987). The effects of cylindrical a shear layer have been examined by Birkinshaw (1991) in the nonrelativistic case, attempted by Urpin (2002) in the relativistic case and the presence of multiple components (jet + sheath + ambient medium) was investigated by Hanasz & Sol (1996, 1998). Other authors have investigated current-driven modes in magnetized jets (Appl & Camenzind 1992; Appl 1996) in addition to KH modes.

An extension of the linear stability analysis in the relativistic case to the weakly nonlinear regime has been performed by Hanasz (1995) and led to the conclusion that Kelvin-Helmholtz instability saturates at finite amplitudes due to various nonlinear effects. An explanation of the nature of the mentioned nonlinearities has been proposed by Hanasz (1997). The most significant effect results from the relativistic character of the jet flow, namely from the fact that the velocity perturbation cannot exceed the speed of light.

A more recent approach is to perform a linear stability analysis in parallel with numerical simulations and to compare the results of both methods in the linear regime and then to follow the nonlinear evolution of the KH instability resulting from numerical simulations. Hardee & Norman (1988) and Norman & Hardee (1988) have made such a study for nonrelativistic jets in the spatial approach and Bodo et al. (1994) in the temporal approach. In the relativistic case this kind of approach was applied for the first time by Hardee et al. (1998) in the case of axisymmetric, cylindrical jets and then extended to the 3D case by Hardee et al. (2001).

Similarly to the linear stability analysis, the numerical simulations of jet evolution can be performed following both the spatial and the temporal approach, depending on the particular choice of initial and boundary conditions. In the temporal approach one considers a short slice of jet limited by periodic boundaries along the jet axis, and adds some specific perturbation, e.g. an eigenmode resulting from the linear stability analysis. Due to the periodic boundary conditions the growing perturbations can only be composed of modes having a wavelength equal to the length of the computational box and/or its integer fractions (Bodo et al. 1994, 1995, 1998). Whereas the spatial approach appears more appropriate to analyze the global dynamics and morphology of the whole jet, the temporal approach is suitable for the comparison between the numerical results and analytical studies of the jet stability because, due to the fact that only part of the jet is simulated, a high effective numerical resolution is achievable with limited computer resources.

Numerical simulations (Martí et al. 1997; Hardee et al. 1998; Rosen et al. 1999) demonstrate that jets with high Lorentz factors and high internal energy are influenced very weakly by the Kelvin-Helmholtz instability. Moreover, Hardee et al. (1998), Rosen et al. (1999) note that contrary to the cases with lower Lorentz factors and lower internal energies, the relativistically hot and high Lorentz factor jets do not develop modes of KH instability predicted by the linear theory. They interpret this fact as the result of a lack of appropriate perturbations generating the instability in the system. In the limit of high internal energies of the jet matter the Kelvin-Helmholtz instability is expected to develop with the highest growth rate.

In this paper we focus on the role of internal energy and the Lorentz factor for the stability of relativistic planar two-dimensional jets in the temporal approach. We chose this simple configuration in order to enhance our ability to understand the complex nonlinear evolution of the unstable KH modes in relativistic jets. For this purpose, we perform a linear analysis of the KH instability and construct perturbations analytically: the eigenmodes of the linear problem representing the most unstable first reflection mode for each set of jet parameters. These modes are then incorporated as initial states of numerical simulations. The choice of planar two-dimensional jets allows us to study symmetric and antisymmetric modes of KH instability which resemble, respectively, pinch and helical modes of cylindrical jets. Only the flutting modes of cylindrical jets, corresponding to azimuthal wavenumbers larger than 1, have no counterparts among the eigenmodes of planar two-dimensional jets. In the present paper we investigate only the symmetric modes.

	Figure 1: The geometry of the flow considered in the linear stability analysis and the numerical simulations, including a description of the boundary conditions.
Open with DEXTER

Our aim is to investigate physical details of the process of transition of KH instability modes from the linear to the nonlinear regime. In the next step (Paper II) we shall investigate the relations of long-term nonlinear evolution of KH instability (including the mixing phase and the properties of the evolved flow at the end of our simulations) on the initial jet parameters.

The paper is organized as follows. In Sect. 2 we describe our method, i.e. perform the linear analysis of the KH instability for the relativistic planar case, then define perturbations and describe numerical algorithm and other details of numerical simulations. In Sect. 3 we present and discuss results of the both linear stability analysis and the numerical simulations. In Sect. 4 we summarize and conclude our paper.

   
2 Method

Following the standard procedure (see e.g. Gill 1965; Ferrari et al. 1978; Hardee 1979) we derive the dispersion relation for the Kelvin-Helmholtz modes. We focus on the simplest geometrical configuration of two-dimensional planar relativistic jets and apply the temporal stability analysis.

2.1 Linear stability analysis

The full set of equations describing the current problem consists of the set of relativistic equations of hydrodynamics for a perfect fluid, (e.g. Ferrari et al. 1978)

$$\gamma^2 \left( \rho + {p\over c^2} \right) \left[ {\partial ... ...abla p + {\vec v \over c^2} {\partial p \over \partial t} = 0,$$ (1)

$$\gamma \left( {\partial \rho\over \partial t} + \vec v \cdotp... ...r \partial t} + \vec\nabla \cdotp (\gamma \vec v) \right] = 0,$$ (2)

and an (adiabatic) equation of state

$$p \rho_0^{-\Gamma} = \mbox{{\rm const}}.$$ (3)

In the preceding equations, c is the speed of light, $\rho_0$ is the particle rest mass density (i.e., $\rho_0 \equiv m n$, where m is the particle rest mass and n the number density in the fluid rest frame). $\rho$ stands for the relativistic density which is related to the particle rest mass density and the specific internal energy, $\varepsilon$, by $\rho = \rho_0(1+\varepsilon/c^2)$. The enthalpy is defined as $w = \rho + p/c^2$, the sound speed is given by $c_{\rm s}^2 = \Gamma p/w$, where $\Gamma$ is the adiabatic index. The relation between pressure and the specific internal energy is $p =(\Gamma - 1) \varepsilon \rho_0$. The velocity of the fluid is represented by $\vec v$ and $\gamma$ is the corresponding Lorentz factor.

The assumed geometry of the jet considered in the forthcoming linear stability analysis and the numerical simulations is sketched in Fig. 1. First of all the considered jet is 2D-planar and symmetric with respect to the x=0 plane. The flow in the jet moves in the positive z direction and its matter forms a contact discontinuity (a vortex sheet) with the matter of external medium at $x= -R_{\rm j}$ and $x=R_{\rm j}$. From now on all quantities representing the jet will be assigned with the "j'' subscript and the quantities representing the ambient medium will be assigned with "a''.

The following matching conditions are imposed on the interface between the jet and the ambient medium

  

$$p_{\rm a} = p_{\rm j} \mbox{for} \vert x\vert=R_{\rm j}$$ (4)

$$h_{\rm a} = h_{\rm j} \mbox{for} \vert x\vert=R_{\rm j}.$$ (5)

The matching conditions express the assumption of equality of the jet and ambient pressures and the equality of transversal displacements of jet ($h_{\rm j}$) and ambient ($h_{\rm a}$) fluid elements adjacent to the jet boundary (at $\vert x\vert = R_{\rm j}$).

In addition, the Sommerfeld radiation condition (expressing the disappearance of perturbations at the infinity) will be applied for linear perturbations.

   
2.1.1 Equilibrium state for the linear stability analysis

We assume that the initial state is an equilibrium configuration. The initial equilibrium can be described by the following set of independent parameters: the Lorentz factor corresponding to the unperturbed longitudinal jet velocity, $v_{\rm j}$, $\gamma = (1 - v_{\rm j}^2/c^2)^{-1/2}$, the particle rest mass density of the jet $\rho_{\rm0j}$ (the particle rest mass density of the ambient medium is normalized to unity: $\rho_{\rm0a} = 1$) and the specific internal energy of the jet $\varepsilon_{\rm j}$. The ambient medium is assumed to be at rest ( $v_{\rm a} = 0$).

The other dependent parameters describing the equilibrium state are: the internal jet Mach number $M_{\rm j} = {v_{\rm j}}/{c_{\rm sj}}$ corresponding to the initial jet longitudinal velocity, the relativistic rest mass density contrast $${\nu = \frac{\rho_{\rm0j}(1+\varepsilon_{\rm j}/c^2)} {\rho_{\rm0a}(1+\varepsilon_{\rm a}/c^2)}}$$ (or, equivalently, the enthalpy contrast $${\eta= \frac{\rho_{\rm0j}\left(1+(\Gamma_{\rm j}-1)\varepsilon_{\r... ...(\Gamma_{\rm a}-1)\varepsilon_{\rm a}/c^2/ (1+\varepsilon_{\rm a}/c^2)\right)}}$$) and the specific internal energy of the ambient medium $${\varepsilon_{\rm a} = \frac{(\Gamma_{\rm j} - 1) \rho_{\rm0j}}{(\Gamma_{\rm a} - 1) \rho_{\rm0a}} \varepsilon_{\rm j}}$$, which is related to the specific internal energy of the jet through the pressure balance condition.

2.1.2 Dispersion relation

The first step towards the dispersion relation is to reduce the Eqs. (1)-(5) to the dimensionless form through the normalization of spatial coordinates to the jet radius $R_{\rm j}$, velocities to the sound speed of the jet material $c_{\rm sj}$, time to the dynamical time  $R_{\rm j}/c_{\rm sj}$ and pressure to the equilibrium pressure.

The next step is to decompose each dependent quantity into the equilibrium value and the linear perturbation. After the reduction of equations to the dimensionless form and substitution of the perturbed quantities in Eqs. (1)-(5), we obtain the following set of dimensionless linearized equations. In the following the dimensionless quantities will be assigned the same symbols as the previous dimensional ones. The subscript "1'' stands for the linear perturbation of the corresponding variable. For the jet medium we get

$$\Gamma_{\rm j} ~ \gamma^2 \left(\frac{\partial{\mbox{{\boldm... ...$ }}_{\rm j}}{c^2} \frac{\partial p_{\rm j1}}{\partial t} = 0,$$ (6)

 

$$\frac{\partial p_{\rm j1}}{\partial t} + \mbox{{\boldmath$v$ }}_{... ...$v$ }}_{\rm j} \cdot \nabla \mbox{{\boldmath$v$ }}_{\rm j1}\right)\right) = 0 ,$$ (7)

where $\mbox{{\boldmath$v$ }}_{\rm j}$ is the initial (unperturbed) jet velocity in units of jet internal sound speed, c is the speed of light in units of the sound speed and the normalized pressure p₀=1 is omitted in the equations. For the ambient medium we get

$$\frac{\Gamma_{\rm j}}{\eta} \frac{\partial {\mbox{{\boldmath$v$ }}_{\rm a1}}}{\partial t} + \nabla p_{\rm a1} = 0,$$ (8)

$$\frac{\partial p_{\rm a1}}{\partial t} + \Gamma_{\rm a}\nabla \cdot \mbox{{\boldmath$v$ }}_{\rm a1} = 0.$$ (9)

The linearized matching conditions (4) and (5) at the jet interface read

  

$$p_{\rm a1} = p_{\rm j1} \mbox{for}$$ |x|=1       (10)

$$h_{\rm a1} = h_{\rm j1} \mbox{for}$$ |x|=1, (11)

where the displacements of fluid elements adjacent to the contact interface in the linear regime are related to transversal velocities by the following formulae

$$v_{{\rm j}x1} = \left( \frac{{\partial}^{} {}}{{\partial} {t}... ...\frac{{\partial}^{} {}}{{\partial} {z}^{}} \right) h_{\rm j1},$$ (12)

$$v_{{\rm a}x1} =\frac{{\partial}^{} {}}{{\partial} {t}^{}} h_{\rm a1}.$$ (13)

The following wave equations can be derived from the Eqs. (6)-(9), respectively for the jet

 

$$\gamma^2 \left(\frac{{\partial}^{} {}}{{\partial} {t}^{}}+v_{\rm ... ...\rm j}}{c^2} \frac{{\partial}^{} {}}{{\partial} {t}^{}}\right)^2 p_{\rm j1} =0,$$ (14)

and for the ambient medium

$$\frac{{\partial}^{2} {p_{\rm a1}}}{{\partial} {t}^{2}} - \fr... ...c{{\partial}^{2} {p_{\rm a1}}}{{\partial} {z}^{2}}\right) = 0.$$ (15)

It is apparent that Eqs. (14) and (15) describe propagation of oblique sound waves in the jet and ambient medium respectively.

The next step is to substitute perturbations of the form

$$\delta_{\rm j1} = \left[ \delta^+_{\rm j} F^+_{\rm j} (x) + ... ...right] \exp i \left(k_{\parallel } z - \omega t\right) + c.c.$$ (16)

with $F^{\pm}_{\rm j} (x) = \exp (\pm i~ k_{{\rm j\perp}} x)$ to describe waves propagating in positive and negative x-directions in the beam. Here we use $k_{\parallel}$ and $k_{a,{\rm j}\perp}$ for longitudinal and transverse wavenumbers. Perturbations in the external medium are of the form

$$\delta_{\rm a1} = \delta^+_{\rm a} F^+_{\rm a} (x) \exp~ i (k_{\parallel } z-\omega t) +c.c.$$ (17)

where $F^+_{\rm a} (x) = \exp (i~ k_{{\rm a}\perp} x)$ for |x|>1. Assuming that ${\rm Re} (k_{{\rm a}\perp}) > 0$, only outgoing waves are present in the ambient medium. After the substitution of the explicit forms of pressure perturbations to Eqs. (14) and (15) we obtain the following expressions for the perpendicular components of wavevectors

$$k_{\rm a \perp} = \left(\frac{\Gamma_{\rm j}}{\eta \Gamma_{\rm a}} \omega^2 - k_\parallel^2\right)^{1/2}$$ (18)

$$k_{{\rm j} \perp} = \left({\omega'}^2-{k'_\parallel}^2\right)^{1/2}$$ (19)

which are standard relations for linear sound waves in both media (note that $\eta \Gamma_{\rm a}/\Gamma_{\rm j}$ is the squared sound speed of ambient medium in units of the sound speed of jet), $\omega' = \gamma(\omega - v_{\rm j} k_\parallel)$ and $k'_{\parallel} = \gamma(k_\parallel - \frac{v_{\rm j}}{c^2} \omega)$ are frequency and wavenumber of the internal sound wave in the reference frame comoving with the jet. The wavevectors $\mbox{{\boldmath$k$ }}_{\rm j} = (k_{{\rm j}\perp}, k_{{\rm j}\parallel})$ and $\mbox{{\boldmath$k$ }}_{\rm a} = (k_{\rm a\perp}, k_{\rm a\parallel})$ determine the direction of propagation of sound waves in the jet and ambient medium respectively. Vanishing of $k_{{\rm j}\perp}$, for instance, would mean that the jet internal sound waves move in the axial direction. In cases $k_{{\rm j}\perp} \ne 0$ and/or $k_{\rm a\perp} \ne 0$ the propagation of the sound waves is oblique with respect to the jet axis.

The corresponding perturbation of the jet surface can be written as

$$h_{\rm a1}=h_{\rm j1}= h \exp~ i (k_{\parallel } z-\omega t) +c.c. ,$$ (20)

which after the substitution to Eqs. (12) and (13) reads

$$v_{{\rm j}x1} = -i \left(\omega - v_{\rm j} k_{\parallel}\right) h_{\rm j1} ,$$ (21)

$$v_{{\rm a}x1} = -i \omega h_{\rm j1}.$$ (22)

With the aid of the above perturbations the whole system of partial differential equations is reduced to a set of homogeneous linear algebraic equations. The dispersion relation appears as solvability condition of the mentioned set of equations, namely it arises from equating the determinant of the linear problem to zero. Within the present setup, the dispersion relation for the Kelvin-Helmholtz instability in supersonic, relativistic, two-dimensional slab jets can be written as

 

$$\frac{1}{\nu \Gamma_{\rm j}} \frac{\omega}{\omega'} \frac{\left({... ...a'}^2 - {k'}^2_{\parallel}\right)^{1/2} & {\rm for}& s=-1\\ \end{array}\right.$$ (23)

where $s = \pm 1$ is the symmetry of perturbation. In the present paper we focus on symmetric perturbations (s=1).

The above dispersion relation is solved numerically with the aid of the Newton-Raphson method (see, e.g., Press et al. 1992). In this way we obtain the complex frequency as a function of the parallel component of the wavenumber $k_{\parallel}$, since we choose the temporal stability analysis for the present investigations.

   
2.1.3 Eigenstates of the system

Once the solutions of the dispersion relation $\omega(k)$ are found, the derivation of eigenstates of the system is straightforward. By utilizing the set of relativistic Eqs. (6)-(9) one can relate perturbations of gas density and velocity to the perturbation of pressure.

First, we choose the amplitude of the pressure wave in the jet, $p_{\rm j}^+$, as a free constant and then relate the other constants to its value. The corresponding "-'' amplitudes are computed by multiplying these values by s in case of pressure and longitudinal velocity and by -s in the case of transversal velocity.

The amplitude of pressure wave in the ambient medium,  $p_{\rm a}^+$, is found from the pressure matching condition at x=1:

$$p_{\rm a}^+\exp\left(i~k_{\rm a\perp}\right) = p_{\rm j}^+ \e... ...p}\right)+s~p_{\rm j}^- \exp\left(-i~k_{{\rm j}\perp}\right).$$ (24)

The remaining amplitudes of velocity perturbations are

$$v_{{\rm a}x}^{+}= \frac{\eta}{\Gamma_{\rm j}} \frac{k_{\rm a\perp}}{\omega} p_{\rm a}^+,$$ (25)

$$v_{{\rm j}x}^{+}= \frac{k_{{\rm j}\perp}}{\Gamma_{\rm j} \gamma \omega'} p_{\rm j}^{+},$$ (26)

$$v_{{\rm a}z}^{+} = \frac{\eta}{\Gamma_{\rm j}} \frac{k_{\parallel}}{\omega} p_{\rm a}^+,$$ (27)

$$v_{{\rm j}z}^{+} = \frac{k'_{\parallel}}{\Gamma_{\rm j} \gamma^2 \omega'} p_{\rm j}^+,$$ (28)

where $\eta$ is the ratio of enthalpies as defined earlier.

Regarding specific internal energy and rest-mass density, perturbations are calculated as follows:

$$\frac{\varepsilon_{a,\rm j1}}{\varepsilon_{a,\rm j}}= \frac{\Gamma_{a,\rm j}- 1}{\Gamma_{a,\rm j}} p_{a,\rm j1}$$ (29)

$$\frac{\rho_{a,\rm j1}}{\rho_{a,\rm j}+1/c^2}=\frac{1}{\Gamma_{a,\rm j}} p_{a,\rm j1},$$ (30)

where $\rho_{a,\rm j}$ refer to the relativistic density in ambient and jet media in dimensionless units, respectively, and $\rho_{a,\rm j1}$is the corresponding perturbation.

2.1.4 Lorentz transformation of velocity perturbations

In the forthcoming discussion of results of numerical simulations we shall analyze the velocity components in the reference frame comoving with the jet. The standard Lorentz transformation formulae for the velocity components are

 

$$v'_z = \frac{v_z - v_{\rm j}}{1- v_z v_{\rm j}/c^2}, v'_x = \frac{{v_x}\sqrt{1-v_{\rm j}^2/c^2}}{1-v_z v_{\rm j}/c^2}\cdot$$ (31)

Assuming that the Lorentz transformation is applied for linear perturbations of sufficiently small amplitude, these transformation rules can be approximated by the following formulae

 

$$v'_{z1} = \gamma^2 v_{z1}, v'_{x1} = \gamma v_{x1},$$ (32)

where v_x1, v_z1 are the velocity perturbations in the reference frame of the ambient medium and v'_x1, v'_z1 are the velocity perturbations in the jet reference frame. These transformation rules tell us that the longitudinal and transversal velocity perturbations in the jet reference frame are larger by factors of $\gamma^2$ and $\gamma$ respectively than the corresponding components in the rest (ambient medium) frame.

   
2.2 Numerical simulations

Now we proceed by performing numerical simulations of jet models to study the growth of unstable modes through the linear and non-linear phases. To this aim, we adopt as an equilibrium initial state the one described in the former section. In order to avoid the growth of random perturbations with wavelengths of the order of the cell size, we replace the transverse discontinuous profiles of equilibrium quantities by smooth profiles of the form

$$v_z(x) = \frac{v_{\rm j}}{\cosh(x^m)}$$ (33)

$$\rho_0(x) = \rho_{\rm0a} - \frac{\rho_{\rm0a} - \rho_{\rm0j}}{\cosh(x^m)}$$ (34)

(see Bodo et al. 1994), where m is the steepness parameter for the shear layer. Typically we used m=40, for which the shear layer has a width ${\sim} 0.1 R_{\rm j}$ (see Appendix A).

Table 1: Equilibrium parameters of different simulated jet models along with solutions of the dispersion relation (23), corresponding to fastest growing first reflection mode, taken as input parameters for numerical simulations. The meaning of the symbols is as follows: $\gamma$: jet flow Lorentz factor; $\varepsilon$: specific internal energy; $c_{\rm s}$: sound speed; p: pressure; $\nu$: jet-to-ambient relativistic rest mass density contrast; $\eta$: jet-to-ambient enthalpy contrast; $M_{\rm j}$: internal jet Mach number; $k_{\parallel , \perp }$: longitudinal/transverse wavenumbers; $\omega _{r,i}$: real/imaginary part of the wave frequency. Labels aand j refer to ambient medium and jet, respectively. The primes are used to assign wavenumber and complex frequency in the reference frame comoving with jet. The last column shows the linear growth rate of KH instability in the jet reference frame expressed in dynamical time units, i.e. in which time is scaled to $R_{\rm j}/c_{\rm sj}$. All the quantities in the table, except the last column, are expressed in units of the ambient density, $\rho _{\rm0a}$, the speed of light, c, and the jet radius, $R_{\rm j}$.

On the other hand the linear solutions derived in Sect. 2.1.3 correspond to an idealized case with a contact discontinuity at the jet interface. Finding of analogous solutions for the corresponding sheared boundary problem is more difficult from the numerical point of view (see Ray 1982; and Choudhury & Lovelace 1984), therefore we decided to implement corresponding smoothed eigenstates obtained for the vortex sheet limit as approximate solutions for the case of sheared boundary. The validity of this procedure is proven a posteriori by the convergence of theoretical growth rates and the ones determined for simulated KH modes. Since, in case of the vortex sheet approximation, all first order dynamical variables (except pressure and the transversal displacement of fluid element adjacent to the jet boundary) are discontinuous, we apply the following smoothing of the equilibrium quantities and linear perturbations at the thin layer surrounding the jet boundary

$$\delta_1=\delta_{\rm a1}-\frac{\delta_{\rm a1}-\delta_{\rm j1}}{\cosh(x^m)},$$ (35)

where $\delta$ stands for any perturbed variable. Such smoothed eigenstates are subsequently implemented as small-amplitude perturbations in the initial equilibrium model of the simulations.

Our first aim is to reproduce the linear evolution of unstable modes by means of hydrodynamical simulations. For this reason and in order to make the picture of evolution of the KH instability as clear and simple as possible, we perturb the equilibrium with the most unstable first reflection mode, which is smoothed at the jet interface according to formula (35). The length of the computational grid is in each simulation equal to the wavelength of the mentioned mode. Defining the initial state in such a way, we know precisely which growth rate to expect in numerical simulations. The difference between the expected and computed growth rates is a measure of the performance of the code in describing the linear regime and gives us confidence on the numerical results from the non-linear evolution.

The simulations have been performed using a high-resolution shock capturing code to solve the equations of relativistic hydrodynamics in Cartesian planar coordinates. The code is an upgrade of that used to study large scale (Martí et al. 1997) as well as compact relativistic jets (Gómez et al. 1997). The equations of relativistic hydrodynamics are written in conservation form and the time variation of the proper rest mass, momentum and total energy within numerical cells are calculated through the fluxes across the corresponding cell interfaces. Fluxes are calculated with an approximate Riemann solver that uses the complete characteristic information contained in the Riemann problems between adjacent cells (Donat & Marquina 1996). It is based on the spectral decomposition of the Jacobian matrices of the relativistic system of equations derived in Font et al. (1994) and uses analytical expressions for the left eigenvectors (Donat et al. 1998). The spatial accuracy of the algorithm is improved up to third order by means of a conservative monotonic parabolic reconstruction of the pressure, proper rest-mass density and the spatial components of the fluid four-velocity, as in GENESIS (see Aloy et al. 1999, and references therein). Integration in time is done simultaneously in both spatial directions using a total-variation-diminishing (TVD) Runge-Kutta scheme of high order (Shu & Osher 1988). See Martí et al. (1997) for the differential equations as well as their finite-difference form, and a description of exhaustive testing of the hydrodynamical code. Besides relativistic density, momentum and energy, the code also evolves a passive scalar representing the jet mass fraction. This allows us to distinguish between ambient and jet matter helping us to characterize processes like jet/ambient mixing or momentum exchange.

The parameters of the simulations presented in this paper are listed in Table 1. Values of the parameters were chosen to be close to those used in some simulations by Hardee et al. (1998) and Rosen et al. (1999) and are chosen to span a wide range in thermodynamical properties as well as beam flow Lorentz factors. In all simulations, the density in the jet and ambient gases are $\rho_{\rm0j}=0.1$, $\rho_{\rm0a} = 1$ respectively and the adiabatic exponent $\Gamma_{\rm j,a}=4/3$.

Since the internal rest mass density is fixed, there are two free parameters characterizing the jet equilibrium: the Lorentz factor and the jet specific internal energy displayed in Cols. 2 and 3 of Table 1. Models whose names start with the same letter have the same thermodynamical properties. Beam (and ambient) specific internal energies grow from models A to D. Three different values of the beam flow Lorentz factor have been considered for models B, C and D. The other dependent parameters mentioned in Sect. 2.1.1 are displayed in Cols. 4-10 of Table 1. Note that given our choice of $\rho_{\rm0j}$, the ambient media associated to hotter models are also hotter. The next three columns show the longitudinal wavenumber together with oscillation frequency and the growth rate of the most unstable reflection mode. The following three columns display the same quantities in the jet reference frame. Next two columns show the perpendicular wavenumbers of linear sound waves in jet and ambient medium respectively. The last column shows the linear growth rate of KH instability in the jet reference frame expressed in dynamical time units, i.e. in which time is scaled to $R_{\rm j}/c_{\rm sj}$. All other quantities in the table, are expressed in units of the ambient density, $\rho _{\rm0a}$, the speed of light, c, and the jet radius, $R_{\rm j}$.

As it is apparent in Table 1 growth rates tend to increase with the specific internal energy of the beam and to decrease with the beam flow Lorentz factor. Note that, in the jet reference frame, models with the same thermodynamical properties tend to have (within $\approx$10%) the same growth rates. Note also that in the jet reference frame and in dynamical time units, all the models have comparable (within a factor 1.5) linear growth rates.

The initial numerical setup consists of a steady two-dimensional slab jet model (see Fig. 1). As stated above, a thin shear layer between the ambient medium and the jet is used instead of the vortex sheet. Due to symmetry properties, only half of the jet (x>0) has to be computed. Reflecting boundary conditions are imposed on the symmetry plane of the flow, whereas periodical conditions are settled on both upstream and downstream boundaries. The numerical grid covers a physical domain of one wavelength along the jet (3 to 20 $R_{\rm j}$; see Table 1) and 100 $R_{\rm j}$ across (200 $R_{\rm j}$ in the case of models D). The size of the transversal grid is chosen to prevent loses of mass, momentum and energy through the boundaries. 400 numerical zones per beam radii are used in the transverse direction across the first 3 $R_{\rm j}$. From this point up to the end of the grid, 100 (200, in case of models D) extra numerical zones growing geometrically have been added. The width growth factor between contiguous zones is approximately 1.08 for models A, B and C and 1.04, for models D. Along the jet, a resolution of 16 zones per beam radius has been used. The agreement of the linear stability analysis with the results of numerical simulations depends on the grid resolution (and the width of the initial shear layer). The applied resolution of $400 \times 16$ grid zones per beam radius is chosen on the base of several tests, which are presented in detail in Appendix A.

The steady model is then perturbed according to the selected mode, with an absolute value of the pressure perturbation amplitude inside the beam of $p_{\rm j}^{\pm}=10^{-5}$. This means that those models with the smallest pressure, like model A, have relative perturbations in pressure three orders of magnitude larger than those with the highest pressures, D. However this difference seems not to affect the linear and postlinear evolution.

Table 2: Times for the different phases in the evolution of the perturbed jet models. In Cols. 2-4 we show: $t_{\rm lin}$ - end of linear phase (the ratios of amplitudes of the different quantities are not constant any longer), $t_{\rm sat}$ - end of saturation phase (the amplitude of the transverse speed perturbation reaches its maximum), $t_{\rm peak}$ - the peak in the amplitude of the pressure perturbation is reached. (see Fig. 2). Note that, as a general trend, $t_{\rm lin} < t_{\rm sat} \leq t_{\rm peak}$. In Col. 5 we show the saturation time in dynamical time units and in the jet reference frame.

Figure 2: Evolution of the relative amplitudes of perturbations. Dotted line: pressure perturbation ( $(p_{\max}-p_0)/p_0$). Dashed line: longitudinal velocity perturbation in the jet reference frame ( $0.5~(v'_{z,\max}-v'_{z,\min})$). Dash-dotted line: perpendicular velocity perturbation in the jet reference frame ( $0.5~(v'_{x,\max}-v'_{x,\min})$). The search for maximum $(p_{\max}$, $v'_{x,\max}$, $v'_{z,\max}$) and minimum ( $v'_{x,\min}$, $v'_{z,\min}$) values have been restricted to those numerical zones with jet mass fraction larger than 0.5. Solid line: linear analysis prediction for the growth of perturbation. Note that abscissae in the last column plots extend up to $t= 1000~R_{\rm j}/c$ and that ordinate values adapt to fit the scale of each plot. Arrows in the plot of Model A05 point to specific stages of evolution used to define $t_{\rm lin}$, $t_{\rm sat}$ and $t_{\rm peak}$ (see text).

Open with DEXTER

   
3 Results

Following the behaviour of simulated models we find that the evolution of the perturbations can be divided into a linear phase, a saturation phase and a mixing phase. This section is devoted to describing the properties of the linear and saturation phases in our models, focusing on the influence of the basic parameters. The remaining fully nonlinear evolution will be discussed in Paper II. Our description shares many points with that of Bodo et al. (1994) for the case of classical jets.

In order to illustrate the growth of perturbations and determine the duration of the linear and saturation phases in our simulations, we plot in Fig. 2 the amplitudes of the perturbations of the longitudinal and transversal velocities, inside the jet and in the jet reference frame, together with the pressure oscillation amplitude. We plot also the growth of the imposed eigenmodes resulting from the linear stability analysis. Both the velocity perturbations are transformed from the ambient rest frame to the unperturbed jet rest frame using the Lorentz transformation rules for velocity components, given by Eq. (31).

We define the characteristic times in the following way. During the linear phase the ratios of the exponentially growing amplitudes of pressure, longitudinal velocity and transversal velocity remain constant by definition. We define the end of the linear phase ( $t_{\rm lin}$) as the moment at which one of quantities starts to depart from the initial exponential growth. Within the set of our models the first quantity to break the linear behaviour is the longitudinal velocity perturbation.

Later on the transversal velocity saturates, i.e. stops growing at the saturation time ( $t_{\rm sat}$). We call the period between $t_{\rm lin}$ and $t_{\rm sat}$ the saturation phase. We find also that the pressure perturbation amplitude reaches a maximum. This moment is denoted by $t_{\rm peak}$. In Paper II we will see that this peak announces the entering of the fully nonlinear regime. The choice of $t_{\rm lin}$, $t_{\rm sat}$, $t_{\rm peak}$ has been illustrated in Fig. 2 (top panel).

Table 2 collects times of the linear and saturation phases in the different models (Cols. 2-4). The last column shows the saturation time in dynamical units and in the jet reference frame. The change of reference frame eliminates the effects coming from the jet flow Lorentz factor that stretches out the rhythm of evolution in the ambient rest frame. Dynamical time units are adapted to the characteristic time of evolution of each model. Hence this change of units and reference frame allows us to compare the relevant scales of evolution of all the models directly.

   
3.1 Linear phase

In the linear phase the ratios of oscillation amplitudes of different dynamical variables (density, pressure and velocity components) are constant. We have used this property to identify the end of the linear phase with the time $t_{\rm lin}$ at which one of the variables deviates from linear growth in a systematic way. Note that our definition is different from that of Bodo et al. (1994) who associate the end of the linear phase with the formation of the first shocks inside the jet.

	Figure 3: Pressure distribution at the end of linear phase for models B05 ( left panel) and D05 ( right panel). The continuous line corresponds to jet mass fraction equal to 0.5 and serves to distinguish jet and ambient media.
Open with DEXTER

We note that in our simulations, during the linear phase, the growth of perturbations of each dynamical variable follows the predictions of the linear stability analysis with relatively good accuracy. On average the growth rates measured for all our numerical experiments are about 20% smaller than corresponding growth rate resulting from the linear stability analysis. This small discrepancy may be partially a result of the application of the shearing layer in simulations and the vortex sheet approximation in the linear stability analysis, and also to a lack of transversal resolution in the numerical simulations (see Appendix A).

In all cases the amplitude of the longitudinal velocity oscillation is the first quantity to stop growing. The reason for the limitation of the velocity oscillations is obvious: the oscillations of velocity components (corresponding to sound waves propagating in the jet interior) cannot exceed the speed of light. This limitation (specific to relativistic dynamics) is easily noticeable in the jet reference frame, but it is obscured by the Lorentz factor (in the first and second power) in the reference frame of the ambient medium, as it follows from the formulae (32).

As seen in Fig. 2, model C20 evolves distinctly from the other models: at t = 270, the linear phase ends before the saturation of longitudinal speed. See the middle panels of Fig. 7 (especially the pressure panel) which could be responsible for the peculiar evolution of model C20. Note also that model C20 is the only one that has a theoretical perturbation growth rate smaller than the numerical one (probably associated with the excited short wavelength mode).

At the end of the linear phase, the values of quantities like density, pressure and flow Lorentz factor are still very close to the corresponding background values (the perturbation in pressure is between 10% and 50% of the background pressure in all the models). Figure 3 shows the distribution of pressure at the end of the linear phase for two representative models B05 and D05.

3.2 Saturation phase

As it has been mentioned, the end of the linear phase coincides with the limitation of the longitudinal oscillation velocity. At time $t_{\rm lin}$ the transversal velocity component is smaller than the speed of light, by an order of magnitude approximately, so the transversal velocity perturbation has still room to grow. We defined $t_{\rm sat}$ as the time corresponding to the saturation of the transversal velocity and saturation phase as the period between  $t_{\rm lin}$ and $t_{\rm sat}$.

We find that the longitudinal velocity perturbation amplitude reaches almost the speed of light (more than 0.9 c) while the transversal perturbation amplitude stops its growth at the level of approximately 0.5c for all presented simulations. This can be explained in the following way. The eigenmodes are built of oblique sound waves overreflected at the jet/ambient medium interface. This means that the amplitude of the reflected wave is larger than the amplitude of the incident one. Sound waves are longitudinal waves, so the gas oscillation velocity can be split into the longitudinal and transversal components separately for the incident and reflected waves. Let us consider locally the incident wave and the reflected wave in the jet medium, in a fixed point close to the jet boundary. Then the longitudinal velocity components of the incident and reflected waves sum in phase, while the transversal ones sum in counterphase. This means that while the amplitude of the total longitudinal velocity oscillation component grows and approaches the speed of light, the oscillation amplitude of the total transversal component is a difference of the two velocities, each smaller than the speed of light. Therefore the oscillation amplitude of the total transversal component has to take values significantly smaller than the longitudinal one.

	Figure 4: Snapshot around saturation of logarithmic maps of pressure, rest-mass density and specific internal energy and non-logarithmic Lorentz factor for model A05.
Open with DEXTER

	Figure 5: Same as Fig. 4 for models B05 ( upper), C05 ( middle) and D05 ( lower).
Open with DEXTER

	Figure 6: Same as Fig. 4 for models B10 ( upper), C10 ( middle) and D10 ( lower).
Open with DEXTER

	Figure 7: Same as Fig. 4 for models B20 ( upper), C20 ( middle) and D20 ( lower).
Open with DEXTER

The duration of the saturation phase depends on the Lorentz factor and the specific internal energy with a tendency to increase with the former and to decrease with the latter (exception made of models C10 or D10 and C20). It ranges between a few tens of (absolute) time units to a few hundreds.

We note that the saturation times expressed in dynamical time units and in the jet reference frame (see Table 2) are almost equal for all models. This similarity can be explained by the following argument. First, the amplitude of pressure perturbations $p_1^\pm=10^{-5}$ is the same for all simulations. Second, the amplitudes of pressure and transversal velocity perturbations are related by formula (26). The perturbation of perpendicular velocity is then transformed to the jet reference frame following Eq. (32). Following the results of numerical experiments, the transversal velocity grows with the linear growth rate until the transversal velocity perturbation reaches the upper limit, i.e. $\vert v_{{\rm j}\perp}^{+}\vert \exp (\omega_i t_{\rm sat}) \simeq 0.5c$. If we express $t_{\rm sat}$ in dynamical time units then in the jet reference frame

$${t'}_{\rm sat}^{\rm dyn} \simeq \frac{\gamma}{c_{\rm sj} \ome... ...t\vert\frac{\omega'}{k_{{\rm j}\perp}}\right\vert\right) \cdot$$ (36)

Since the term $\omega'/k_{\rm j\perp}$ under the logarithm is close to the jet sound speed and is varying only slightly in between models, while the other factors under the logarithm are constants, the saturation time is proportional to the inverse growth rate. From the linear stability analysis, the latter one measured in the jet reference frame and expressed in dynamical time units is almost equal for all the models (see last column in Table 1). Substitution of numbers into the formula (36) provides $t'_{\rm sat}$ in the range $10 \div 15$, which is consistent with the value 13.5 in Table 2. However, the remarkable convergence of saturation times for all the models is difficult to explain in view of the additional randomness in the evolution of perturbations, which is apparent in Fig. 2.

During the saturation phase, the jet inflates (Bodo et al. 1994 called this phase the expansion phase) and deforms due to transversal oscillations. On the other hand, the saturation time  $t_{\rm sat}$coincides (within a few time units) with the appearance of an absolute maximum in the pressure distribution (at $t_{\rm peak}$) at the jet boundary, and with the start of the mixing phase (to be studied in Paper II). Figures 4-7 display snapshots of several quantities close to the end of the saturation phase before the distortion of the jet boundary.

The structure of perturbations at the end of the saturation phase is quite similar in all models. The longitudinal wavelength of perturbations is given by the linear stability analysis and it is constant because of the fixed length of the computational domain. Along with the longitudinal wavelength, the opening angle of oblique waves, far from the jet symmetry plane, is given by the linear analysis. The opening angle of oblique waves is enlarged in the vicinity of the jet interface. Closer to the jet symmetry plane the perturbations are already in the nonlinear phase, which can be recognized by the apparent presence of oblique shock fronts in the jet itself and the ambient medium. These shocks form as a natural consequence of the nonlinear steepening of oblique sound waves as their amplitude becomes large. Practically all the four quantities displayed for a given model (i.e. pressure, rest-mass density, Lorentz factor and internal energy) show the same structure, so there is no need to discuss them separately. Therefore we treat the pressure distribution as a representative quantity. The following properties of the flow patterns can be noticed in Figs. 4-7:

1.

the interface between jet and ambient medium forms a regular sinusoidal pattern. The amplitudes are all comparable and have values of the order of $0.5 {-} 1.0 R_{\rm j}$. In all cases the departure from a regular sinusoidal shape is rather small;

2.

the structure of oblique shocks crosses the jet interface in such a way that the deformed interface is almost parallel to the shock front. One should notice however that the oblique shock moves in the flow direction in the external medium and in the counter-flow direction in the jet medium, as can be deduced from the distribution of hot post-shock and cold pre-shock material. These two oppositely moving shock waves fit together because of the jet background flow;

3.

the highest gas pressure is always located on the jet symmetry plane (the brightest spots on pressure plots), where the crossing shocks form the familiar x-pattern, but the pressure enhancement on the jet interface is almost as strong as on the jet symmetry plane;

4.

the shocks are much stronger (i.e., the jumps between pre- and post-shock pressure larger) in the cold cases. In the case of hot models the jumps are smoother. This is a straightforward consequence of the larger ratios of the perturbated velocity to the sound speed in the cold cases;

5.

there is no significant variation of the structure of nonlinear perturbations for different values of the Lorentz factor;

6.

model C20 displays a structure of internal waves, which is absent in other cases. This fine structure is probably connected to the distinct evolution of this model in the postlinear phase (see Sect. 3.1 and Fig. 2).

Therefore the differences between models just before the end of the saturation phase can be considered as minor.

   
4 Summary

This is the first of a series of papers devoted to the effects of relativistic dynamics and thermodynamics in the development of KH instabilities in planar relativistic jets, accross both linear and fully nonlinear regimes. To this aim, we have performed a linear stability analysis and high-resolution numerical simulations for the most unstable first reflection modes in the temporal approach, for three different values of the jet Lorentz factor $\gamma$ (5, 10 and 20) and a few different values of specific internal energy of the jet matter (from 0.08 to 60.0 c²). In this paper we concentrate in the early stages of evolution of the KH instability, namely the linear and saturation phases. The present paper is also intended to set the theoretical and numerical background for the whole series of papers.

Our simulations describe the linear regime of evolution of the excited eigenmodes of the different models with high accuracy. The growth rates of the perturbed modes in the vortex sheet approximation were determined with an average relative error of 20%.

In all the examined cases the longitudinal velocity perturbation is the first quantity that departs from linear growth when it reaches a value close to the speed of light in the jet reference frame. The reason for this saturation, specific to relativistic dynamics, is not so obvious in the reference frame of the external medium where it saturates at a smaller value ( ${\approx} c/\gamma$, where $\gamma$ is the bulk flow Lorentz factor in the jet).

The saturation phase extends from the end of the linear phase up to the saturation of the transversal velocity perturbation (at approximately 0.5c in the jet reference frame). The saturation times for the different numerical models have been explained from elementary considerations, i.e. from properties of linear modes provided by the linear stability analysis and from the limitation of the transversal perturbation velocity.

The limitation of the components of the velocity perturbation at the end of the linear and saturation phases allows us to conclude that the relativistic nature of the flow appears to be responsible for the departure of the system from linear evolution. This behaviour is consistent with the one deduced by Hanasz (1995, 1997) with the aid of analytical methods.

At the saturation time the perturbation structure is close to the structure of the initial perturbation (the one corresponding to the most unstable first reflection mode), except that the oblique sound waves forming the perturbation became steep due to their large amplitude. It is interesting to note that the oblique shocks are stronger (i.e., the pressure jumps are larger) in the colder cases.

Our simulations, performed for the most unstable first reflection modes, confirm the general trends resulting from the linear stability analysis: the faster (larger Lorentz factor) and colder jets have smaller growth rates. As we mentioned in the Introduction, Hardee et al. (1998) and Rosen et al. (1999) note an exception which occurs for the hottest jets. These jets appear to be the most stable in their simulations (see also the simulations in Martí et al. 1997). They suggest that this behaviour is caused by the lack of appropriate perturbations to couple to the unstable modes. This could be partially true as fast, hot jets do not generate overpressured coocons that let the jet run directly into the nonlinear regime. However, from the point of view of our results, the high stability of hot jets may have been caused by the lack of radial resolution that leads to a damping in the perturbation growth rates. We note that the agreement between the linear stability analysis and numerical simulations of KH instability in the linear range has been achieved for a very high radial resolution of 400 zones/$R_{\rm j}$, which appears to be especially relevant for hot jets. Finally, one should keep in mind that the simulations performed in the aforementioned papers only covered about one hundred time units, well inside the linear regime of the corresponding models for small perturbations.

The high accuracy of our simulations in describing the early stages of evolution of the KH instability (as derived from the agreement between the computed and expected linear growth rates and the consistency of the saturation times) establishes a solid basis to study the fully nonlinear regime, to be done in the second paper of this series. In this paper (Paper II) we will show that the similarities found in the evolution of all the models accross the linear and saturation phases is lost and very different nonlinear evolution is found depending on the initial jet parameters.

Acknowledgements

This work was supported in part by the Spanish Dirección General de Enseñanza Superior under grant AYA-2001-3490-C02 and by the Polish Committee for Scientific Research (KBN) under grant PB 404/P03/2001/20. M.H. acknowledges financial support from the visitor program of the Universidad de Valencia. M.P. has benefited from a predoctoral fellowship of the Universidad de Valencia (V Segles program). The authors want also to acknowledge financial support from the Spanish-French Picasso program and the Polish-French Joumelage program in the period 1997-98 when the research on relativistic jet stability was initiated.

References

• Aloy, M. A., Ibáñez, J. M $^{\underline{\rm a}}$, Martí, J. M $^{\underline{\rm a}}$, & Müller, E. 1999, ApJS, 122, 151 [NASA ADS] [CrossRef] (In the text)
• Appl, S. 1996, A&A, 314, 995 [NASA ADS] (In the text)
• Appl, S., & Camenzind, M. 1992, A&A, 256, 354 [NASA ADS] (In the text)
• Blandford, R. D., & Pringle, J. E. 1976, MNRAS, 176, 443 [NASA ADS] (In the text)
• Birkinshaw, M. 1984, MNRAS, 208, 887 [NASA ADS] (In the text)
• Birkinshaw, M. 1991, MNRAS, 252, 505 [NASA ADS] (In the text)
• Bodo, G., Massaglia, S., Ferrari, A., & Trussoni, E. 1994, A&A, 283, 655 [NASA ADS] (In the text)
• Bodo, G., Massaglia, S., Rossi, P., et al. 1995, A&A, 303, 281 [NASA ADS] (In the text)
• Bodo, G., Rossi, P., Massaglia, S., et al. 1998, A&A, 333, 1117 [NASA ADS] (In the text)
• Choudhury, S. R., & Lovelace, R. V. E. 1984, ApJ, 283, 331 [NASA ADS] [CrossRef] (In the text)
• Donat, R., Font, J. A., Ibáñez, J. M $^{\underline{\rm a}}$, & Marquina, A. 1998, J. Comput. Phys., 146, 58 [NASA ADS] [CrossRef] (In the text)
• Donat, R., & Marquina, A. 1996, J. Comput. Phys., 125, 42 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
• Fanaroff, B. L., & Riley, J. M. 1974, MNRAS, 167, 31 [NASA ADS] (In the text)
• Ferrari, A., Massaglia, S., & Trussoni, E. 1981, MNRAS, 198, 1065 [NASA ADS] (In the text)
• Ferrari, A., Trussoni, E., & Zaninetti, L. 1978, A&A, 64, 43 [NASA ADS] [MathSciNet] (In the text)
• Ferrari, A., Trussoni, E., & Zaninetti, L. 1981, MNRAS, 196, 1051 [NASA ADS]
• Font, J. A., Ibáñez, J. M $^{\underline{\rm a}}$, Marquina, A., & Martí, J. M $^{\underline{\rm a}}$ 1994, A&A, 282, 304 [NASA ADS] (In the text)
• Gill, A. E. 1965, Phys. Fluids, 8, 1428 [NASA ADS] (In the text)
• Gomez, J. L., Martí, J. M $^{\underline{\rm a}}$, Marscher, A. P., Ibáñez, J. M $^{\underline{\rm a}}$, & Marcaide, J. M. 1995, ApJ, 449, L19 [NASA ADS]
• Gómez, J. L., Martí, J. M $^{\underline{\rm a}}$, Marscher, A. P., Ibáñez, J. M $^{\underline{\rm a}}$, & Alberdi, A. 1997, ApJ, 482, 33 [NASA ADS] [CrossRef] (In the text)
• Hanasz, M. 1995, Ph.D. Thesis, Nicholas Copernicus University, Torun (In the text)
• Hanasz, M., & Sol, H. 1996, A&A, 315, 355 [NASA ADS] (In the text)
• Hanasz, M., & Sol, H. 1998, A&A, 339, 629 [NASA ADS] (In the text)
• Hanasz, M. 1997, in Relativistic jets in AGNs, ed. M. Ostrowski, M. Sikora, G. Madejski, & M. Begelman, Kraków, 85 [arXiv:astro-ph/9711275] (In the text)
• Hardee, P. E. 1979, ApJ, 234, 47 [NASA ADS] [CrossRef] (In the text)
• Hardee, P. E. 1984, ApJ, 287, 523 [NASA ADS] [CrossRef] (In the text)
• Hardee, P. E. 1986, ApJ, 303, 111 [NASA ADS] [CrossRef] (In the text)
• Hardee, P. E. 1987, ApJ, 318, 78 [NASA ADS] [CrossRef] (In the text)
• Hardee, P. E. 2000, ApJ, 533, 176 [NASA ADS] [CrossRef]
• Hardee, P. E., & Norman, M. L. 1988, ApJ, 334, 70 [NASA ADS] [CrossRef] (In the text)
• Hardee, P. E., Rosen, A., Hughes, P. A., & Duncan, G. C. 1998, ApJ, 500, 559 [NASA ADS] (In the text)
• Hardee, P. E., Hughes, P. A., Rosen, A., & Gomez, E. A. 2001, ApJ, 555, 744 [NASA ADS] [CrossRef] (In the text)
• Lobanov, A. P., & Zensus, J. A. 2001, Science, 294, 128 [NASA ADS] [CrossRef] (In the text)
• Martí, J. M $^{\underline{\rm a}}$, Müller, E., Font, J. A., Ibáñez, J. M $^{\underline{\rm a}}$, & Marquina, A. 1997, ApJ, 479, 151 [NASA ADS] [CrossRef] (In the text)
• Mirabel, I. F., & Rodriguez, L. F. 1999, Ann. Rev. A&A, 37, 409 [NASA ADS] (In the text)
• Norman, M. L., & Hardee, P. E. 1988, 334, 80 (In the text)
• Payne, D. G., & Cohn, H. 1985, ApJ, 191, 655 [NASA ADS] (In the text)
• Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. 1992, Numerical recipes (Cambridge University Press) (In the text)
• Ray, T. P. 1982, MNRAS, 198, 617 [NASA ADS] (In the text)
• Rosen, A., Hughes, P. A., Duncan, G. C., & Hardee, P. E. 1999, ApJ, 516, 729 [NASA ADS] [CrossRef] (In the text)
• Roy Choudhury, S., & Lovelace, R. V. E. 1984, ApJ, 283, 331 [NASA ADS] [CrossRef]
• Sari, R., Piran, T., & Halpern, J. P. 1999, ApJ, 519, L17 [NASA ADS] [CrossRef] (In the text)
• Shu, C. W., & Osher, S. J. 1988, J. Comput. Phys., 77, 439 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
• Trussoni, E., Ferrari, A., & Zaninetti, L. 1983, in Astrophysical jets. Proceedings of the International Workshop, Turin, Italy (Dordrecht: Reidel), 281 (In the text)
• Turland, B. D., & Scheuer, P. A. G. 1976, MNRAS, 176, 421 [NASA ADS] (In the text)
• Urpin, V. 2002, A&A, 385, 14 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Zensus, J. A. 1997, ARA&A, 35, 607 [NASA ADS] (In the text)

 

  
Online Material

Appendix A: Influence of the numerical resolution in the description of the linear regime of KH modes of relativistic flows

Figure A.1: Linear growth of the amplitude of the pressure perturbation (in logarithmic scale) versus time (in units of $R_{\rm j}/c$) for different resolutions. rx stands for transversal resolution, rz for longitudinal resolution, and m is the value which gives the shear layer steepness.

Growth of the instability depends critically on the numerical viscosity of the algorithm. Hence our first aim was to look for suitable numerical resolutions by comparing numerical and analytical results for the linear regime. We performed a number of simulations based on model A05 changing longitudinal and radial resolution, and also the exponent m for the shear layer steepness (see Eqs. (33) and (34)). Figure A.1 shows results for several of those simulations.

A shear layer was included in order to avoid non-steadiness in initial conditions given by discontinuous separation between the beam and the external medium. Therefore, if for some transversal resolution we introduced a shear layer which was too thin, we found a similar non-steady behavior. However, in order to reproduce the linear regime, we needed a steep enough shear layer, due to the fact that theory was developed for a discontinuous separation between both media. We can see in Fig. A.1 that m=10 models are considerably damped with respect to the theoretical growth, independent of the transversal resolution. Hence resolution perpendicular to the flow appeared to be essential, requiring very high resolutions (400 zones/$R_{\rm j}$) and thin shear layers (m=40) with 40 to 45 zones, i.e., an equilibrium between steepness and number of cells. Lower transversal resolutions and/or thicker shear layers led to non-satisfactory results, with a slow or damped growth. A very low longitudinal resolution results in damping of growth rate, too, as can be seen from comparison of m=20 models, but r_z=16 seems to give a reasonable growth rate compared to theory. This (small) resolution of 16 zones/$R_{\rm j}$ along the jet was taken as a compromise between accuracy and computational efficiency.