{'fraction_non_alphanumeric': 0.02719264649559556, 'fraction_numerical': 0.016006000255330013, 'mean_word_length': 4.503688740558581, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interactions of Jets with Inhomogeneous Cloudy Media

We present two-dimensional slab-jet simulations of jets in inhomogeneous media consisting of a tenuous hot medium populated with a small filling factor by warm, dense clouds. The simulations are relevant to the structure and dynamics of sources such as Gigahertz Peak Spectrum and Compact Steep Spectrum radio galaxies, High Redshift Radio Galaxies and radio galaxies in cooling flows. The jets are disrupted to a degree depending upon the filling factor of the clouds. With a small filling factor, the jet retains some forward momentum but also forms a halo or bubble around the source. At larger filling factors channels are formed in the cloud distribution through which the jet plasma flows and a hierarchical structure consisting of nested lobes and an outer enclosing bubble results. We suggest that the CSS quasar 3C48 is an example of a low filling factor jet - interstellar medium interaction whilst M87 may be an example of the higher filling factor type of interaction. Jet disruption occurs primarily as a result of Kelvin-Helmholtz instabilities driven by turbulence in the radio cocoon not through direct jet-cloud interactions, although there are some examples of these. In all radio galaxies whose morphology may be the result of jet interactions with an inhomogeneous interstellar medium we expect that the dense clouds will be optically observable as a result of radiative shocks driven by the pressure of the radio cocoon. We also expect that the radio galaxies will possess faint haloes of radio emitting material well beyond the observable jet structure.


INTRODUCTION
Previous models of jets propagating through the interstellar medium of radio galaxies have almost always assumed a uniform, usually hot (T ∼ 10 7 K) interstellar medium. These models have contributed substantially to our understanding of the propagation of jets in evolved radio galaxies wherein much of the dense matter that could have previously obstructed the jets has been cleared away. However, it is now clear that in many young radio galaxies, i.e. Gigahertz Peak Spectrum (GPS) and Compact Steep Spectrum (CSS) sources, substantial interaction between the expanding radio plasma and dense interstellar clouds occurs. This interaction has a number of effects: Two of the most important are that it significantly distorts the radio morphology and it also gives rise to shock excited line emission. Strong observational evidence for the latter effect has been obtained from combined Hubble Space Telescope (HST) and radio observations of a number of CSS sources (de Vries et al. 1999;Labiano et al. 2003;O'Dea et al. 2003) showing a strong alignment effect between optical and radio emission and clear evidence for disturbance of the optically emitting gas by the expanding radio lobe. On the other hand, the dense gas cannot have a large filling factor in the case of GPS and CSS sources since the expansion speed of the radio plasma is usually quite high: ∼ 0.1 − 0.4 c (Conway 2002b;Murgia et al. 1999;Parma et al. 1999;Murgia et al. 2002;Murgia 2003). Thus the picture that is conveyed by the combined radio and optical observations is of radio plasma propagating through a predominantly hot and tenuous interstellar medium in which dense clouds are embedded with a low filling factor. Nevertheless, as we show in this paper, even a low filling factor of dense gas can have a pronounced effect on the evolution of the radio lobes.
High redshift radio galaxies provide another example in which radio plasma -cloud interactions are clearly important. For example, in the high redshift radio galaxy 4C41.17, interactions between the jets and large clouds in the Lyman-α emitting halo of the parent galaxy lead to shock-excited line emission and to accelerated rates of star formation (Bicknell et al. 1998). Returning to the low-redshift Universe, radio plasma -clumpy ISM interactions may be responsible for the heating of the ISM that is necessary to counteract the otherwise large mass inflow rates that are inferred from X-ray observations of cooling flows.
Thus, there is ample physical motivation for detailed investigations of the manner in which jets and the related radio lobes interact with realistic interstellar media. One of the most significant drivers of this research is that all of the scenarios discussed above represent examples of AGN-galaxy feedback which, it has been argued (e.g. Silk & Rees 1998;Kawakatu et al. 2003), is a significant determinant of the ultimate structure of galaxies and in particular, the relation between the mass of the central black hole and the mass of the galaxy bulge (Magorrian et al. 1998;Tremaine et al. 2002).
In this first paper in this series we explore the nature of the interaction between jets and their lobes through a series of two-dimensional slab jet simulations. The advantages of this approach is that two dimensional simulations offer the prospect of high resolution and an initial rapid investigation of parameter space, which capture some of the features of a more realistic three dimensional simulations. Fully three dimensional simulations will be presented in future papers.
The choice of coordinate system is non-trivial when representing supersonic jets in two spatial dimensions. Many aspects of a jet burrowing through a uniform external medium have been adequately described in the past by axisymmetric simulations. However an axisymmetric simulation is less appropriate for physical systems where the jet is expected to be bent or deflected. When colliding directly with a dense cloud, an axisymmetric jet must either burrow through the cloud or cease its advance. However, a two-dimensional slabjet can be deflected, although both jet and cloud are essentially infinite in the direction perpendicular to the computational grid. Our present work is concerned with systems where jet-cloud collisions and indirect interactions dominate. We therefore opt for a Cartesian grid and simulate transient slab jets.

Jet energy flux, cloud and ISM density
The choice of jet parameters is based on a self-similar model for lobe expansion (Bicknell et al. (1997), based upon the model by Begelman (1996)). In this model, the lobe is driven by a relativistic jet which expands in such a way that ζ the ratio of hot spot to lobe pressure is constant. A value ζ ∼ 2 was inferred by Begelman from simulations by Lind et al. (1989). Let ρ0 be the atmospheric mass density at a fiducial distance R0 from the core, let δ be the index of the power variation of atmospheric density with radius (ρ = ρ0(R/R0) −δ ) and let the jet energy flux be FE. The temperature of the external medium is unimportant in this model since it is the inertia of the background medium that governs the expansion of the lobe. The advance speed, v b of the head of the bowshock as a function of distance R from the core is given by: (1) where the jet energy flux is 10 46 FE,46 erg s −1 , n(100 pc) is the number density at 100 pc and the numerical values are for ζ = 2 and δ = 1.5. Typical expansion speeds of CSS and GPS sources are of order (0.05 − 0.4) c with the largest velocities measured for the most compact sources (Conway 2002b;Murgia et al. 1999). Conway (2002a) has reviewed estimates of number densities ∼ 1cm −3 for three GPS sources, and so a fiducial value of n(100 pc) = 1 cm −3 is reasonable. Compatibility of the expansion speeds estimated from equation (1) and the observations indicates jet energy fluxes in the range of 10 46−47 erg s −1 . Jets of this power are also consistent with a conventional ratio ∼ 10 −12 − 10 −11 Hz −1 of 1.4 GHz monochromatic radio power to jet energy flux and a median power at 5 GHz ∼ 10 27.5 W Hz −1 . Nevertheless, in other sources which morphologically resemble CSS radio galaxies, such as the less radio powerful Luminous Infrared Galaxies studied by Drake et al. (2003), it is apparent that the jet powers are somewhat less. Thus, in our simulations, we have adopted jet powers of 10 45 and 10 46 erg s −1 .
Hubble Space Telescope observations (O'Dea et al. 2002;de Vries et al. 1997de Vries et al. , 1999 show abundant evidence for dense clouds of emission-line gas associated with CSS radio galaxies and quasars. For a temperature 10 4 Tc,4 K of the dense clouds, a number density nx ∼ 1 cm −3 and temperature 10 7 Tx,7 K for the hot interstellar medium, the number density in the clouds is nc ≈ 10 3 nxTx,7/Tc,4 cm −3 , assuming that the initial, undisturbed clouds are in pressure equilibrium with their surroundings. If we adopt such a value for the density of the interstellar medium then equation (1) implies that the velocity of advance would be ∼ 0.02 c -much lower than what is observed. Hence, as discussed in § 1, the dense clouds must have a small filling factor. Accordingly, in the simulations described below, we create clouds of approximately the inferred density and investigate the effect of various filling factors.
With the inferred radio source and cloud geometry, the line emission orginates from dense, relatively cool clouds that are shocked as they are overrun by the expanding radio source. A model of this type was proposed by de Vries et al. (1999) and elaborated by Bicknell et al. (2003). In both treatments, the excess pressure of the advancing radio lobe drives shocks into the dense clouds. In the simulations described below, this is seen to occur but some direct jet-cloud collisions are also apparent. For radio lobe driven shocks (the main source of shocks), the shock velocity is: where nc and nx are the cloud and ISM number densities respectively. Shocks in the velocity range of 300−950kms −1 are produced for bow-shock advance speeds ∼ 0.1 c and density ratios nc/nx ∼ 10 4 − 10 3 . This is in accord with the observed velocity dispersions and velocity offsets inferred from HST spectra of a sample of CSS sources (O'Dea et al. 2002).
Radiative shocks do not develop instantaneously; the time required for a shock to become fully radiative is effectively the time, t4, required for post-shock gas to cool to ∼ 10 4 K. In order that there be a display of optical emission lines, t4 must be less than a fraction f ∼ 1 of the dynamical time associated with the expansion of the lobe. In estimating t4 we utilize a cooling function, which estimates the non-equilibrium cooling following a high velocity shock, and which is derived from the results of the one-dimensional MAPPINGS III code (Sutherland & Dopita 1993). This cooling finction is also used in our simulations described below. The condition on t4 becomes a condition on the bowshock velocity . This upper limit on v b is insensitive to the various parameters, in particular, f . If condition (3) is not satisfied near the head of the radio lobe it may be satisfied at the sides where the bow-shock expansion speed is about a factor of 2 lower. In this case one would expect to see emission from the sides but not the apex of the radio source. This condition gives a limiting bow-shock velocity v b ∼ 0.1c. Therefore, it would be interesting to see whether the more rapidly expanding sources have a lower ratio of emission line flux to radio power, and whether the emission line flux is distributed in a different manner to that in lower velocity sources. These conditions may also be related to the causes of the trailing optical emission in CSS sources (de Vries et al. 1999). There are two possibilities: (a) The apex of the bow-shock is non-radiative (if v b exceeds the above limit) or (b) The number of dense clouds per unit volume decreases at a few kpc from the core, i.e. beyond the normal extent of the narrow line region. We are inclined to favour option (b) but (a) should also be kept in mind. The simulations capture these characteristics that depend on the time of the effective onset of radiative shocks in the clouds. We see a range of shock conditions varying from almost pure adiabatic shocks to strongly radiative. When the cooling time of gas becomes shorter than the Courant time for the simulation, the integration of the energy equation ensures that the radiative shock occupies approximately one pixel in which the gas temperature is equal to the predefined equilibrium temperature (10 4 K in these simulations).

Overpressured jets -Mach number and density ratio
In order to carry out realistic 2D simulations, we wish to establish jets with powers per unit transverse length that are relevant to the types of radio galaxies mentioned in the introduction, in combination with other reasonable values of parameters such as jet width and density. If the number density of the interstellar medium is ∼ 1 cm −3 on 100 pc scales, then the jet must be overpressured. We first show this in the context of a sub-relativistic jet model and then support the proposition with a relativistic jet model. Let a nonrelativistic jet have a density η times that of the external medium, an internal Mach number M , a specific heat ratio γ, a pressure ξ times the external ISM pressure and a cross-sectional area Aj. Further, let ρx and vx be the density and isothermal sound speed of the ISM. Then the jet mass, momentum and energy fluxes are, collecting the dimensional quantities on the left: Equation (6) clearly shows that the energy flux can be made arbitrarily large by decreasing the density ratio η. However, there is a natural lower limit on η -the value at which an intially relativistic jet becomes transonic. The critical value is defined by the condition ρjc 2 /4pj = 1 (Bicknell 1994), implying that ηcrit = 4ξ(vx/c) 2 ≈ 5.9 × 10 −6 ξTx,7. With the overpressure factor, ξ = 1, η = 10 −5 , a jet diameter of 10 pc, Tx = 10 7 K and M = 20, the energy flux is only ∼ 10 44 erg s −1 . One requires an overpressure factor of about 100 to achieve 10 46 erg s −1 .
This result is easily verified for a relativistic jet whose energy flux is given by where χ = ρjc 2 /4pj, and ρj and pj are the jet rest-mass density and pressure respectively. If ξ = 1, Γ = 10, a jet diameter of 10 pc, Tx = 10 7 K and χ = 0 (the latter being appropriate for a jet dominated by relativistic plasma), then the energy flux is only 1.2 × 10 43 erg s −1 . Hence, for a relativistic jet, we require ξ ∼ 10 2−3 in order to achieve the required jet power. Is an overpressured jet physical? The overpressure ratio refers to the pressure with respect to the undisturbed interstellar medium. However, the lobe itself is over-pressured as a result of the advancing bow shock so that it is possible for the jet to be in approximate equilibrium with the ambient lobe. In a real galaxy the over pressured lobe may not extend back to the core, in which case an overpressured jet would behave initially like a classical overpressured, underexpanded jet.
Therefore, in our simulations, we have adopted low values of the density ratio, η ∼ 10 −3 , Mach numbers, M of the order of 10 and a high overpressure ratio in order to obtain the necessarily high jet energy fluxes (or their equivalent in 2D). Such a low value of η has the additional effect of producing superluminal velocities (β ∼ 5 − 6) in some simulations, highlighting the shortcomings of simulating powerful jets with a non-relativistic code. However, the essential point is that it is mainly the jet power that is physically important in jet-cloud interactions.
One could adopt the strategy of adopting nonrelativistic parameters that correspond more closely to those of relativistic jets, for example by equating the velocity and pressure and matching either the power or thrust (e.g. see Komissarov & Falle (1996) and Rosen et al. (1999)). We have not done this explicitly here since there is no transformation that perfectly maps a non-relativistic flow into a relativistic equivalent. Clearly, this initial investigation of parameter space needs to be followed by more comprehensive simulations utilising a fully relativistic code.

Prescription for cloud generation
In our simulations we prescribe a uniform density for the hot interstellar medium. Within this medium we establish an initial field of dense clouds that mimics the distribution of clouds in the interstellar medium, with a random distribution of positions, shapes and internal density structures. We consider that this is more desirable than specifying, say circular clouds, albeit with random centres; in that case the subsequent evolution may be governed by the initial symmetry. An advantage of the approach that we describe here is that the distribution of the sizes, shapes and directions of propagation of radiative shocks is more realistic than that obtained from initially highly symmetric clouds.
We prescribe the density and distribution of clouds using a Fourier approach. In the first step we establish a grid of complex numbers in Fourier space and assign a random phase to each point. The largest wave number is compatible with the spatial resolution of the simulations. An amplitude filter is applied, ∝ |k| −5/3 , in order to mimic the power spectrum of structures formed by turbulence. Taking the Fourier transform of the wavenumber data produces an array of density fluctuations in real space. These fluctuations, σ(r), are normalised so that their mean is 0 and their variance is 1. The final number density of gas is then defined by the relation: where nx is the number density of hot inter-cloud gas, nmin is the minimum number density of gas in a dense cloud, nσ is a scaling factor and σmin is a threshold value determining whether a particular cell is part of a cloud or else hot intercloud medium. In practice we choose the mean density of cloud gas,n and set the value of nσ implicitly. We apply upper and lower limits to the density distribution in k-space. Amplitudes are zeroed for k < kmin, where kmin corresponds to the Jeans wave number. This prevents the insertion of clouds that would have collapsed to form stars and also avoids the generation of a cloud field dominated by a single large cloud-mass. Amplitudes are also zeroed for k > kmax where kmax corresponds to 3 or 4 pixels of spatial resolution. Fine features, at the 3-4 pixel level, can introduce excessive and artificial sensitivity to effects on the level of individual pixels.
Finally, radial limits are applied to the cloud field as a whole. An outer limit, rmax, (typically ∼ 0.5 − 5 kpc) limits the extent of the cloud field. The inner radius, rmin, is defined by the Roche limit in the vicinity of a black hole: If the black hole mass M bh = 10 9 M⊙ and the cloud density ρc ≈ 15M⊙ pc −3 then rmin ≈ 250 pc. Table 1 presents a summary of the properties of the clouds and intercloud medium, and the extent of cloud coverage in our simulations (which are each described in specific detail in §3). Figures 1 and 2 represent the initial distributions of cloud and intercloud gas in three of our simulations. Figures 3 and 4 show how the masses of individual clouds are distributed in the initial conditions of four of our simulations. As the cumulative distribution functions show, most clouds have masses of less than a few thousand solar masses. The scatter plots show the characteristic depth of clouds: in mean terms the distance from an interior point to the surface, r = DI/S where D is the number of spatial dimensions, I is the interior (area of a 2D cloud) and S is the surface (perimeter of a 2D cloud). A small r means that most of the cloud mass is well exposed to surface effects, such as ablation and the penetration of shocks from outside. If r is large than the interior regions are typically deeper from the locally nearest surface, and correspondingly better insulated from effects at the surface.
Large clouds, near the Jeans limit (e.g. thickness 2r ∼ 105 pc in the B series) are rare in every simulation. The simulations with a greater cut-off σmin, and hence greater filling factor of clouds, also feature individual clouds of greater mass than the simulations with lower σmin.

Code
In the simulations reported here, we have used a locally modified version of a Piecewise Parabolic Method (Colella & Woodward 1984) code, known as VH-1, published on the web by the computational astrophysics groups at the Universities of Virginia and South Carolina 1 . Our modifications of the VH-1 code include restructuring of the code for efficient computation on the computers at the Australian National University (formerly a Fujitsu VP300; now a Compaq Alphaserver SC), addition of code for dynamically important cooling and other code for coping with a numerical shock instability . The code has already been used in several adiabatic applications (e.g. Saxton et al. 2002a, Saxton et al. 2002b) and most recently in two-dimensional simulations of fully radiative wall shocks ). In the latter paper, cooling was dynamically important and the code for preventing numerically unstable shocks was implemented and proved to be important in elminating physically unrealistic features. We refer to this augmented code as ppmlr.
We have recently ported the ppmlr code to the Message Passing Interface (MPI) environment (see Appendix A). This port enables us to take advantage of the Compaq Alphasever parallel processor architecture. The MPI port is described in the appendix. The simulations described here typically involve 16 processors. Future three-dimensional simulations will regularly use 256 processors and in exceptional circumstances, up to 512 processors.

SIMULATIONS
In the following subsections we describe a number of simulations that exhibit a wide range of phenomena when we vary fundamental parameters such as the jet energy flux (FE), the jet density ratio (η) and the filling factor of clouds (fc) obstructing the path of the jet. The jet parameters in the simulations are presented in Table 2, and the characteristics of the clouds and intercloud media are given in Table 1. Figure 1. Logarithmically scaled map of density relative to the external medium, ρ/ρx, in the initial cloud field of the simulation A1. The dark gray region represents the intercloud medium and the lighter gray regions represent cloud cells. Each cloud has an internal density structure defined by the Fourier construction described in the text ( §2.3). The region shown extends over 5 kpc × 5 kpc (1600 × 1600 cells), with the jet origin at the middle of the left edge. Table 1. Parameters and properties of the inter-cloud medium and initial cloud fields in the respective simulations. The X-ray emitting inter-cloud medium has temperature Tx and particle number density nx. Clouds have a minimum density n c,min , mean densitync. The cloud generation threshold is σ min and the clouds occupy a total area Ac, corresponding to a filling factor of fc. The total exposed cloud perimeter is Pc.
12 In each simulation the X-ray emitting inter-cloud medium has temperature Tx = 10 7 K, and the corresponding isothermal sound speed is our unit of velocity in the simulations, v0 = vx = 365 km s −1 . Our unit for time measurement is t0 = v0/1 pc = 2.68 × 10 3 yr. The jet is defined by a group of cells on the left boundary where the density, pressure and velocity are kept constant. Elsewhere on the left boundary a reflection condition is applied, representing the effect of bilateral symmetry of a galaxy with equally powerful jet and counter-jet. The boundary cells on other sides of the grid are subject to an open boundary condition: in each cycle of the hydrodynamic computation, the boundary cell contents are copied into the adjacent ghost cells. This one-way boundary allows hydrodynamic waves to propagate off the grid.
In a journal paper, one is necessarily limited to showing illustrative snapshots from simulations. Movies of the simulations described in the following section, which show the detailed dynamical evolution, may be obtained from http://macnab.anu.edu.au/radiojets/.

Model A1: high density, overpressured jet
This simulation exhibits both expected and unexpected features of jets interacting with a distribution of clouds in an active nucleus. The jet is initially unimpeded by clouds and generally behaves like a jet in a uniform medium showing the usual bow-shock, terminal jet-shock and backflow. Then two small clouds are overtaken by the bow-shock and the overpressure in the cocoon drives radiative shocks into them. Subsequently the jet backflow is distorted by the reverse shock emanating from the upper of the two clouds. The interaction with the cloud field becomes more dramatic when the apex of the bow shock strikes one of the larger clouds just below the jet axis. Again this drives a strong radiative shock into the cloud and the larger reverse shock deflects both the jet and backflow sideways. The bow-shock then intersects the larger of the two clouds and this produces an additional reverse shock that deflects the jet further sideways. The jet continues onwards past these two clouds whlist the dense clouds in the wake of the bow-shock (i.e. immersed within the jet cocoon) are gradually disrupted by the gas surrounding them. A large amount of turbulence is generated within the jet cocoon by the reverse shocks from the bow-shock impacted clouds and the jet flow is clearly affected but not disrupted. At late times many of the clouds immersed within the cocoon show trails resulting from the ablation of gas within the turbulent cocoon. However, the details of these features will be uncertain until we are able to resolve them fully and subdue the effects of numerical mass diffusion.
Notable features of this interaction are: (I) The jet does not directly touch any of the clouds. The radiative shocks within the clouds are all produced as a result of the overpressure behind the bow-shock. (II) The importance of the reverse shocks within the cocoon in deflecting the jet and making the cocoon more turbulent. (III) The survival of the jet despite these interactions (aided by heaviness of the jet in this simulation, η = 0.3835).
The left panel of Figure 5 shows the jet interacting with clouds during an early stage of the simulation, t = 5t0. The bow shock and therefore the influence of the jet have only propagated about 1.7 kpc forward and 0.5 kpc laterally, and most of the clouds remain in their initial condition. The upper panels of Figure 6 show details of the subregion immediately surrounding the jet at this time of the simulation. Within the high-pressure region bounded by the bow shock, all clouds are shocked. The cloud shocks are the intensely cooling regions on the cloud perimenters visible in the upperright panel of Figure 6. The cloud density is sufficiently high that the shocks have penetrated only a few pc at this stage -a small depth compared to a cloud thickness. The other noteworthy cooling component visible in the cooling map ( Figure 6) is the X-ray emission from the bow shock. This is more intense near the apex of the bow shock than at its sides.
The lower panels of Figure 6 show the density and cooling at a later time, t = 10t0, when the jet is more advanced and the clouds are more eroded. Some clouds that lay in the path of the jet at time 5t0 (located at (0.85, 0.15) kpc and (0.95, 0.10) kpc) are now completely shocked and destroyed. Other clouds have survived closer to the nucleus but with less damage because they are out of the jet's path, e.g. those at (0.3, 0.4) kpc and (0.3, 0.25) kpc. The jet survives an oblique interaction with one group of eroded clouds near (0.9, 0.15) kpc, but its width increases to the right of this point. A second cloud collision near (1.7, 0.0) kpc deflects the jet into the negative y direction, and causes its decollimation into a turbulent bubble of hot plasma accumulating at the apex of the bow shock. The head of the jet is unsteady and changes direction completely within a few t0.
The density distribution at time t = 19t0, near the finish of the simulation, appears in the right panel of Figure 5. Clouds are visible in all stages of shock compression. The two major cloud collisions, around (1.00, 0.05) kpc and (1.75, 0.00) kpc, now leave the jet well collimated and relatively straight. The jet loses collimation at a third colli- Figure 2. Density maps of the initial cloud fields of the simulations B1, B2 and B3 (from top to bottom respectively). These cloud fields were generated from the same fractal distribution but use a different threshold σ min to define the cloud/intercloud boundaries. The entire grid covers 4 kpc × 4 kpc (2000 × 2000 cells), but the square subregion shown above is 1.2 kpc on a side. The jet origin is at the middle of the left edge of each panel.
sion, around (2.5, − 0.2) kpc. The turbulent bubbles around the head of the jet have grown, and spurt between the clouds in several different directions. As was the case at t = 5t0 and 10t0, the clouds situated off to the sides have survived much better than clouds that lay close to the direct path of the jet.
In summary, whether a cloud behaves radiatively and contracts, or adiabatically and ablates is affected by its size, location and density. Small clouds are more effectively adiabatic and less durable than large clouds. This may be partly because their exposed surfaces are large compared to their internal volume, and partly because they are geometrically narrower and on average less dense than large clouds. Some of the smallest swell and develop ablative tails almost as soon as they enter the cocoon, e.g. those between (1.3, − 0.7) kpc and (2.1, − 0.7) kpc when t = 10t0, in Figure 6. A larger cloud near (0.2, − 0.5) kpc was within the cocoon for most of the simulation and yet it was barely affected, even until t = 19t0 (right panel, Figure 5). Qualitatively, each cloud appears to suffer one of three possible fates. Small clouds are destroyed quickly regardless of their position. Large clouds near the path of the jet affect the jet through reflection shocks (even if the jet doesn't strike directly) but they inflate and ablate away in the long term. Large bystander clouds situated well off the jet axis endure almost indefinitely Interaction with clouds alters the propagation and shape of the bow shock. The bow shock is severely retarded when it strikes a cloud and refracts through the gaps between clouds. Thus the bow shock surface is conspicuously indented ("dimpled") for as long as it remains amidst the cloud field. We would expect this effect to manifest in 3D simulations also, although the extra degree of freedom in 3D is also likely to bring about some differences.
The jet generally remains identifiable and collimated after deflection in two or at most three cloud interactions. These interactions need not be direct: it is sufficient for the jet to interact with the reverse shocks that reflect off a cloud lying within several times rj off the jet's path. These backwashing shocks spread continually in the turbulent cocoon; they appear to act like a "cushion" against collision between jet and cloud. More than three consecutive cloud interactions disrupts the jet.

Low density, overpressured jets
3.2.1 Jet stability properties, model B0: Our main series of simulations depicts a more powerful jet (FE = 10 46 erg s −1 ) with a much lower density contrast (η = 1.16 × 10 −3 ) and a large overpressure (ξ = 100). The initial internal Mach number of the jet is M = 12.72 at the nucleus. The jet radius is rj = 10 pc, the square grid is 4 kpc wide and the resolution is 2000×2000 cells. The background gas again has a temperature Tx = 10 7 K but the density is raised to nx = 1 cm −3 . These parameters are suggested by the observations of GPS sources (e.g. in the summary by Conway 2002a).
One of the aims of this work is to distinguish the effects of an inhomogeneous medium from the phenomena that are intrinsic to the jet itself. Therefore we have performed a control simulation containing only the jet and an initially uniform intercloud medium. The main large-scale features at the end of this simulation (see Figure 7) are: an outer shell of highly overpressured gas processed through the bow shock, enclosing a turbulent cocoon where jet plasma mixes with external gas. The length of the jet is small compared to the diameter of the bubble surrounding it. The jet is highly unstable: it stays recognisably straight throughout its first 0.1 kpc but at greater distances from the origin it flaps transversely with increasing amplitude, losing collimation at positions beyond ∼ 0.5 kpc from the origin. This highlights a feature of two-dimensional simulations that may not be present in three dimensions. The instability in this jet is driven by 2D turbulence in the cocoon that results from the mixing of gas at the edge, although the instability is not as rapid as in the ensuing cases with cloudy media. The turbulent driving is probably stronger here because of the persistence of vortices associated with the well-known inverse cascade in two dimensions. The instability evident here is seeded by a very small amount of numerical noise in the initial data.
Hardee & Norman (1988) performed a linear stability analysis of the instabilities of two-dimensional slab jets. Two qualitative kinds of instabilities occur. An unstable pinching mode is responsible for diamond shocks at quasi-periodic locations along the length of the jet. Sinusoidal instabilities are responsible for transverse oscillations, which disrupt and ultimately decollimate the jet. According to the analytical and numerical results of Hardee & Norman (1988), the maximally growing fundamental mode has wavelength an e-folding length scale of and angular frequency ω * ≈ 1.5 αsjae/rj , where ae is the adiabatic sound speed in the cocoon immediately exterior to the jet surface, the parameter αsj → 1 for M ≫ 1, and η, in this case, is the density contrast between the jet and its immediate surroundings within the cocoon. In simulation B0, the internal Mach number within the jet is reduced to typically M ≈ 6, downstream of the first diamond shock (see Figure 7). The cocoon gas immediately surrounding the jet has a density typically ∼ 4ρj and adiabatic sound speed ae ∼ 800v0, where v0 ≈ 365 km s −1 is the isothermal sound speed of the undisturbed medium outside the bow shock. Upon substitution into equations (10) and (11), these values imply exponential growth of sinusoidal instabilities on scales l ∼ 100 pc and with wavelengths λ ∼ 200 pc. This is consistent with the growing oscillations that are apparent in individual frames of our simulations, e.g. in Figure 7. From equation (12), the implied period of the instability is T * < 0.052 t0. The simulated jet thrashes back and forth transversely within this typical time.
The sinusoidal transverse instability revealed in this homogeneous ISM simulation are similar in wavelength and growth time to instabilities manifest in the inhomogeneous simulations discussed below. However, the appearance of this instability in this "clean" simulation shows that the instability can be properly regarded as a Kelvin-Helmholtz instability of the jet itself -albeit, in the inhomogeneous cases, stimulated by the turbulence in the cloudy medium. Note that growth and wavelength of the transverse KH instability at first sight seems to be suggestive of a transonic jet. However, the comparison with the Hardee & Norman Figure 6. Subregions of model A1 showing the area disturbed by the jet at the times t = 5t 0 (upper row) and 10t 0 (lower row). The left panels show density ρ/ρx and the right panels show velocity vectors superimposed on maps of the volumetric rate of cooling. The cooling is measured in units of erg cm −3 s −1 , but the upper end of the scale is deliberately saturated to keep the bow shock visible. The actual maximum rate of radiative cooling in the cloud shocks is 1.8 × 10 3 erg cm −3 s −1 when t = 5t 0 . Reflected and refracted shocks are visible at distances of order 10 pc around some of the clouds.
(1988) expressions for maximally growing wavelength and growth time shows that our simultions are consistent with the expected behavior of supersonic jets with the quoted Mach numbers.

Small cloud filling factor: model B1
This is our simulation with the fewest clouds: their filling factor is fc = 4.1 × 10 −3 (σmin = 2.6) within a region spanning radially between 0.1kpc and 0.5kpc from the nucleus. All of the clouds are quickly overtaken by the expanding bow shock. As in model A1, each cloud is shocked by the overpressure of the cocoon (up to 5 × 10 4 times greater the ex- ternal pressure), and contracts as it cools radiatively. Some clouds are accelerated away from their initial positions by the thrust of the jet. For some other clouds, radiative cooling eventually proves inadequate; the tendency to contract is overcome, and the cloud's fate is expansion and ablation. Halfway through this simulation at time t = 24t0, (top panel, Figure 8), shows the destruction of one such cloud in the path of the jet, situated at (0.6, 0.0) kpc but ablating away across about 0.2 kpc to the right. Meanwhile, a few extremely dense and compact clouds persist in apparent safety off to the sides of the jet, within a few hundred pc of the nucleus. They appear to survive indefinitely and appear scarcely changed by the time the simulation ends, as seen in the middle panel of Figure 8. In comparison with Figure 7 we see much more dense gas (derived from the ablated clouds) intermixed with jet plasma in the cocoon.
The last panel in Figure 8 represents the appearance of the jet in synchrotron emission at late stages of the simulation. Here we have used the approximation (e.g. Saxton et al. 2001 that the emissivity jν ∝ ϕp (3+α)/2 , where p is the pressure, ϕ is the relative concentration of mass originating from the jet (a scalar tracer that is passively advected in our code) and α is the spectral index that we take to be 0.6. Several diamond shocks are visible in the well collimated section of the jet closest to the nucleus. The brighter, barred "fishbone" shaped pattern seen at distances x > 0.4 kpc is a transient and typical manifestation of the transverse instability of the jet (as in § 3.2.1).

Intermediate cloud filling factor: model B2
In this case the filling factor of clouds is higher with fc = 3.0 × 10 −2 , corresponding to σmin = 1.9. In the first few t0 of simulation time the jet is completely disrupted within < 0.1 kpc of the nucleus as a result of the KH instability induced by waves reflected off the innermost clouds. By t ∼ 12t0 (see top panel, Figure 9), the jet has ablated the nearest clouds and cleared an unobstructed region out to ∼ 0.1 kpc on either side. Within this region the jet remains well collimated during the rest of the simulation. The clouds that define the sides of this sleeve are shocked but they survive until the last frame (middle panel, Figure 9). The flow of plasma from the head of the jet splits into a fork on either side of a cloud group that sits nearly on the jet axis at about x ≈ 0.5 kpc. Throughout the simulation these clouds retain dense, efficiently cooling cores but they continually lose dense gas entraining into the stream of jet plasma passing to either side. Near the end of the simulation (t = 48t0, Figure 9) these obstructive clouds are trailed by dense ablated gas extending throughout a triangular area up to 0.5 kpc further out in the positive x direction.
The radio emissivity (last panel, Figure 9) has some similarities to the appearance of model B1, with a bright jet remaining well collimated out to around 0.5 kpc. However the turbulent plumes of jet-derived plasma beyond the jet's head are more disrupted by the dense obstructions of the cloud field, Radio bright plasma divides in several directions, as Figure 9 shows. The brightest areas of the jet and post-jet flows are typically concentrated closer to the nucleus than in the less cloudy model B1.

Large filling factor: model B3
This is our simulation with largest filling factor of clouds, fc = 0.12 and σmin = 1.2. As a result of the high coverage of clouds, the jet is disrupted when it is less than 0.1 kpc from its origin. Jet plasma streams through the available gaps between the persistent clouds (see the density distribution in right panel of Figure 11). These low-density, radio-bright channels evolve as the cloud distribution changes. A channel may close when nearby clouds expand ablatively or are blown into positions where they create a new obstruction. Turbulent eddies of jet plasma accumulate in pockets within the cloud field, most conspicuously in places where a channel has closed. The density image, Figure 10, illustrates the initial development of intercloud channels.
The left panel of Figure 11 shows simulated radio emissivity in the inner regions of the galaxy at the end of the simulation. The radio-bright material is concentrated near the nucleus within an irregular, branched region that is not directly related to the alignment of the initial jet. This particular configuration of channels is ephemeral, occurring after the opening and closure of previous channels as the cloudfield evolves.
Eventually, one or more channels pushes into the open region outside the cloudy zone. A spurt of turbulent jet plasma accumulates at that point. In later stages of our simulations, these eruptions eventually become large compared to the cloudy zone, and merge with adjacent eruptions to form a diffuse, turbulent, approximately isotropic outer cocoon composed of a mixture of jet plasma, ablated cloud gas and intercloud gas processed through the bow shock.

Effect of cooling: model B3a
We ran a further simulation with the same cloud field as in B3 but with radiative cooling switched off (see Figure 12). In the absence of radiative cooling, the compression of clouds is insignificant. The outflow of jet plasma rapidly turns the clouds into ablating, flocculent banana shapes. Large clouds and small clouds ablate efficiently. The survival time of individual clouds is shorter than in the simulation including the effects of radiative cooling. Dense material is more readily ablated from the clouds, and this may block channels of jet plasma more readily or at an earlier time, but otherwise the radio morphology is qualitatively similar to that of B3.

Evolution of the outer part of the radio source
In each case, as a consequence of the disruption of the jet, its momentum is distributed throughout the cocoon in an isotropic way. Hence, the long term evolution is similar to that of an energy-driven bubble. Figure 13 presents plots of the average radius of the furthermost extent of the radio plasma versus time, for represenative simulations. The plots also show the corresponding analytic solution for a two dimensional bubble with the same parameters of energy flux per unit length and density (grey line, see Appendix B for a derivation). It is evident that the analytic slope of the radius/time curves is replicated very well by the simulations. The offset may be attributed to the time taken for the isotropic phase to be established. This may be related to the finite minimum width of the jet at its nozzle, the size of the innermost nuclear region that is initially cloudless, and the intrinsic length scale of the jet's instability ( § 3.2.1).
We expect that a similar power-law expansion will occur in three dimensions, although this of course is yet to be confirmed. However, given that jet momentum and power can be redistributed isotropically in three dimensions, it is reasonable to expect that the long term evolution of the outer bubble will be described by a three-dimensional energy-driven bubble.

COMPARISON OF MODEL B1 WITH 3C48
Whilst these simulations have been two-dimensional, some of their features should persist in three-dimensional simulations with similar parameters. In particular the trends with filling factor will probably persist, to some extent, so that comparisons with existing observationsla are of interest, even at this stage. We therefore make a tentative comparison of one of these simulations (model B1) with the CSS quasar 3C48. We can also draw some general conclusions from these simulations, concerning the likely overall morphology of young radio sources.
Simulation B1 has the lowest filling factor of dense clouds and impedes the jet the least. The jet is still disrupted but the flow associated with it continues with reduced collimation and with its forward momentum spread over a much larger region. In the disrupted section, transverse oscillations are produced, with shocks created at points where the jet bends significantly, leading to local enhancements in emissivity. A large scale bubble is produced outside the cloudy region but the jet also retains some coherency in the intermediate region. This particular simulation bears a strong qualitative resemblance to 3C48. In Figure 14 a snapshot of the simulated radio emissivity at a particular time is compared with the 1.6 GHz VLBI image of 3C48 (Wilkinson et al. 1991). In two dimensions, the alternating direction of the jet results from the transverse instability that we have already discussed ( §3.2.1). In three dimensions, the corresponding instability would be a helical m = 1 mode. In both cases the alternating direction of the supersonic flow produces shocks. We see evidence for these, in both the snapshot of the emissivity from the 2D simulation and the image of the jet.
This particular simulation predicts that there will be a diffuse halo of emission external to the jet region of the source. This was, in fact, detected in larger scale images of 3C 48. One such image presented in Wilkinson et al. (1991), shows diffuse emission on a 1 arcsecond scale in 3C48 extending well beyond the compact jet structure evident in the VLBI image. More recent higher quality images have been obtained by Ojha (2004) and one of their contour maps is reproduced in Figure 14. These images clearly reveal the extended structure of 3C48. Wilkinson et al. (1991) had already concluded that the morphology of 3C48 is likely the result of jet-cloud interactions and cited evidence for asymmetric emission line gas northwest of the core as well as the presence of dust revealed by the high IRAS flux. They interpreted the radio knot close to the core as the result of a jet-cloud interaction. Whilst this is possible, our simulation shows that the sort of jet disruption observed in 3C48 can occur without direct jet-cloud encounters. The general turbulence in the jet cocoon can excite jet instabilities, pro- ducing shocks in the supersonic flow either as a result of the jet being bent through an angle greater than the Mach angle, or though pinching-mode instabilities. (Both of these occur in the simulations.) An explanation for the first bright knot involving cocoon turbulence may therefore be more favorable since direct jet-cloud encounters are usually highly disruptive and Wilkinson et al. (1991) appealed to a glancing deflection in order to explain the radio structure. Given that the sort of jet-ISM interaction revealed in our simu-lations, is occurring, we expect that the emission line gas would show evidence for shock excitation as a result of radiative shocks induced by the overpressured radio cocoon. Chatzichristou (2001) has suggested that the emission line spectrum of 3C48 may show evidence for a mixture of AGN and shock excitation. Wilkinson et al. (1991) also noted out that 3C48, classified as a CSS quasar, is much less extended than other similar power radio-loud quasars. They speculated that the   confinement results from the emission line gas and dust inferred to be present in the ISM. Qualitatively, we agree with this, except to point out that the dispersion of the momentum flux (by a small number of clouds) prevents a radio source from expanding rapidly. As we have shown, the source becomes an expanding bubble with a much slower rate of expansion than a jet with the same kinetic power. Hence the gas and dust is not required to physically confine the radio plasma. The slower expansion can be effected by a small filling factor of dense clouds that redistribute the jet momentum isotropically.
We also note a qualitative resemblance between the simulation B3 and M87. The channeling of the flow in B3 produces a hierarchy of structures resembling the hierarchy of radio structures in M87 apparent in the low frequency radio image (Owen et al. 2000). However, we have not emphasized this with a detailed comparison since the scale of the source and the simulation are quite different. Nevertheless, the production of the structure observed in simulation B3 depends mainly on the filling factor of the clouds and this mechanism should also work on larger scales. We propose to confirm this in future with a fully relativistic 3D simulation.
Other sources with similar morphology to 3C48 and which may be being disrupted in a similar way include the TeV γ-ray blazars, MKN 501 and MKN 421. The most recent VLBI images of these sources (Giovannini et al. 2000;Piner et al. 1999) show similar structures in the form of disrupted jets and large scale plumes suggesting that these sources are disrupted on the sub-parsec scale by a cloudy medium. Again the scale of the simulation B1 is such that it does not strictly apply. However, the physical principle of these jets being disrupted by a small filling factor of dense clouds is probably still valid with interesting repercussions for the production of the γ-ray emission (Bicknell et al. 2005).

DISCUSSION
We have presented a series of two-dimensional simulations of slab jets propagating through inhomogeneous media. De-spite the limitations of the restrictions to two dimensions, and to non-relativistic flow, these simulations provide useful insight into the characteristics of jets in radio galaxies which have an inhomogeneous interstellar medium. These include GPS and CSS radio sources, high redshift radio galaxies and galaxies at the center of cooling flows. In all of these galaxies there is ample evidence for dense gas clouds that may interact with the radio source producing at least some of the features evident in our simulations. In the case of GPS and CSS sources the inhomogeneous medium may be the result of a fairly recent merger; in the case of high redshift radio galaxies an inhomogeneous medium is likely to be the result of the galaxy formation process at those epochs, as large galaxies are assembled from sub-galaxian fragments.
As well as the limitations to two dimensions and nonrelativistic flow, there is another feature of our simulations that requires some comment and that is the ablation of cloud material that is evident in most cases. There is likely to be a large contribution of numerical mass diffusion due to the very high density contrasts between jet and cloud gas. Although this has some effect on the evolution, the dominant processes in the case of a low filling factor are waves emanating from the clouds impacted by the high pressure cocoon and, in the high filling factor simulations, direct obstruction of the flow plus wave effects.
Despite the above-mentioned caveats, we have been able to deduce these general features from our simulations: (i) Jet instability. When propagating through the type of inhomogeneous interstellar medium that we have imposed, the jets become unstable. However, this instability does not always result from direct jet-cloud collisions but occurs indirectly from the waves generated by the impact of the high-pressure jet cocoon on the dense clouds. This leads to a general level of turbulence in the medium that induces unstable Kelvin-Helmholtz oscillations in the jet. An unstable jet was also produced in the homogeneous simulation (B0) and this was related to turbulence caused by the entrainment of ISM into the cocoon as well as the higher level of large scale turbulence in two dimensions. However, the simulations with an inhomogeneous medium became unsta- Figure 14. The radio features of 3C48 compared with simulation B1. The left panel is a 22 GHz (K-band) VLA A-array image illustrating the large-scale (arcsecond scale) structure of 3C48 (Ojha 2004). This image shows diffuse emission on a larger scale to that of the jet, consistent with the prediction of the simulation. The middle panel is a 1.6 GHz VLBI image featuring the jet and bright knots (Wilkinson et al. 1991). At the redshift of 3C48 (z = 0.3670) an angular scale of 0.1 arcsec corresponds to 0.514 kpc (assuming H 0 = 70 km s −1 Mpc −1 , Ωm = 0.27 and Ω Λ = 0.73). The right panel shows a snapshot of radio emissivity from the B1 simulation, at time t = 48t 0 . Left and middle images are reproduced with the kind permission of R. Ojha and P. Wilkinson respectively. ble in a shorter distance and we attribute this to the effect of the clouds.
(ii) Production of halo/bubble. The instability of the jet is associated with the spreading of the jet's momentum and a turbulent diffusion of the jet material into a halo or bubble. These are characteristics that we expect to be present in three dimensions. The radius-time curve for the outer part of the halo is well-described by the analytic expression for a two-dimensional bubble. In three dimensions a energy driven bubble model should also be appropriate. This is useful since it gives us a way of estimating the energy flux fed into a large scale radio source from the dimensions of the outer halo and the density of the surrounding interstellar medium. The production of a bubble means that the overall increase of size of the radio source with time is diminished with respect to a source that propagates unimpeded. However, this restriction on the size of the source is not due to direct physical confinement of the radio plasma by denser gas. Rather, it is the result of the isotropic spread of the jet momentum by the process of jet instability described above.
(iii) The effect of filling factor. The main parameter that we have varied in these simulations is the filling factor of the dense clouds and we can see that this makes a considerable difference to the morphology of the radio source.
With only a small filling factor of clouds, the jet is readily disrupted, but retains some directed momentum and produces the 3C48-type morphology described in § 4. As the filling factor increases, the jet propagates via a number of different channels which open and close leading to an hierarchical structure of radio plasma that we noted is similar to the large-scale structure of M87.
(iv) The effect of jet density ratio. We also varied the jet density ratio parameter, η, but only in one simulation. The simulation A1 had η = 0.38 whereas in the other simulations η ∼ 10 −3 . The main noticeable effect is that the heavier jet is affected less by instabilities and direct jet-cloud collisions than the lighter jets. This points to an important difference between models with η > ∼ 0.01 and models with η < ∼ 0.001. Of course, such low density ratios are important to provide the relevant power for radio galaxy jets. However, reassuringly for physics but difficult for modelers, such low values of η push models into the relativistic regime, as we have discussed in § 2.2. There are two important points to be made here. The first is that many jet simulations usually involve a modest value of η 0.01 and these may be misleading. Second, the most realistic simulations will be those that directly include relativistic physics and which also incorporate the mixing of relativistic and thermal plasma. This paper complements other work on the physics of direct jet-cloud collisions (e.g. Fragile et al. (2004), Wiita et al. (2002)) in which the focus has been on the effect of a jet on a single cloud. Direct jet-cloud collisions are clearly important in the context of sources such as Minkowksi's object (Van Breugel et al. 1985). In this paper we have clearly been more concerned with the overall radio source morphology, the evolution of an ensemble of clouds, the effects of indirect jet-cloud collisions (i.e. the effect of the high pressured cocoon on dense clouds) and trends with cloud filling factor. However, we note an interesting connection between Wiita et al. (2002) and this work: We have both noted that less jet disruption occurs in the case of slow dense jets. We also note that the morphology of radiative cloud shocks deduced by Fragile et al. (2004) is also reproduced in some of the individual clouds in our simulations.
Our simulations have clearly all been two-dimensional. What are the expected differences between two-dimensional and three-dimensional work? There are at least two (1) In two dimensions, vortices survive for a lot longer as a result of the inverse cascade of turbulent energy. These vortices play a significant role in disrupting the jet and clouds. We expect much less disruption from vortical motions in three dimensions.
(2) The number of independent directions that the jet can take in three dimensions is greater so that we expect that the clouds will impede the jet less severely and that escape from the cloudy region surrounding the nucleus will be easier. These issues as well as the inclusion of special relativistic effects will be addressed in future papers.

ACKNOWLEDGMENTS
Our research has been supported by ARC Large Grant A699050341 and allocations of time by the ANU Supercomputer Facility. Professor Michael Dopita provided financial support from his ARC Federation Fellowship. We are very grateful to P.N. Wilkinson and R. Ojha for their permission to use their respective images of 3C48. We thank the anonymous referee for the several suggestions improving the format of the paper. cate between themselves using functions of the MPI (Message Passing Interface) library 2 .
If the global grid has dimensions of I × J cells in the x and y directions, then we define a "long subgrid" as an area with I × (J/N ) cells (see the left panel of Figure B1), and a "wide subgrid" as an area of (I/N ) × J cells (see the right panel of Figure B1). Within the memory at each processor we maintain representations of one long subgrid or one short subgrid: see for example the shaded long and wide subgrids held by processor 3 in Figure B1. The total number of cells in a long subgrid is the same as that of a wide subgrid. In fact for each hydrodynamic variable on each processor, we identify the long and wide subgrid arrays with the same section of memory, which is feasible and reliable since the program uses only one of the subgrid representations at any time (as we shall detail below).
We then conceptually subdivide each long or wide subgrid into N blocks separated by transverse cuts (e.g. the dotted line divisions in processor 3's subgrids depicted in Figure B1). For each processor p the block numbered p is the only area of overlap between the long and wide subgrids that it is assigned.
At the start of the simulation each processor is assigned long and wide subgrids according to its identity number. These regions are filled according to the initial conditions which are defined in terms of global grid coordinates. The first hydrodynamic sweep (Colella & Woodward 1984) is in the x direction, and is performed by each processor within its long subgrid arrays (see left panel of Figure B2). Thus the x-sweeps do not encounter any of the subgrid interfaces, obviating the need for any communication of flux information between processors at this stage. The one-dimensional array vectors associated with each sweeep are the maximum possible length, which makes optimal use of processor cache.
After the x sweep is complete, the evolved density, velocity, pressure and other zone data are communicated between the processors. By use of the MPI ALLTOALL function, the data of block i of the long subgrid of processor j is copied to block j of the wide subgrid of processor i, for all i and j. These are corresponding areas of the global computational grid. As the communication sequence proceeds, each processor receives data blocks and assembles them into complete, ordered arrays describing the wide subgrid that it is conceptually assigned. The wide subgrid arrays overwrite the memory areas that was previously used for long subgrid data. However the sequence of the MPI ALLTOALL communication provides for a complete redistribution of data across the set of processors without any collective loss of information.
When the communication is complete, each processor performs a hydrodynamic sweep in the y direction within its respective wide subgrid arrays. (See right panel of Figure B2.) Again, as in the x sweep, no processor boundaries are encountered during the computation, and the onedimensional sweep arrays are long. At the end of the y sweep, the processors communicate blocks of their evolved wide grid arrays in another MPI ALLTOALL operation, ensuring that the long subgrid arrays are up to date. At this point program control operations are performed, e.g. generating output files B2 Temporal evolution of bubble

B2.1 Formulation of equations
Consider a cylindrical bubble that is being fed by a jet with an energy density per unit length (transverse to the jet) equal to A. We consider the jet energy to be effectively thermalised so that the pressure derived from the shocked jet plasma drives cylindrical expansion. Let the energy density of the bubble be ǫ and the pressure of the bubble be p.
We represent the bubble in two dimensions in the x − y plane and take a length, L of bubble in the z−direction. The total energy, E = ǫπR 2 L of this length of bubble with volume πR 2 L is governed by the equation: d dt ǫπR 2 L + p d dt πR 2 L = AL (B5) Figure B1. Distribution of the long and the wide subgrid partitions in a simulation run with 8 processors. For illustrative clarification, the long subgrid assigned to processor 3 is highlighted in the left panel, and the wide subgrid assigned to processor 3 is highlighted in the right panel. Each processor's memory needs only to retain one long subgrid array and one wide subgrid array, plus scalar global variables. Collective maps of the global computational grid need only be assembled from output data files after the simulation is run. Figure B2. Hydrodynamic sweeps in the x direction are only performed in the long subgrid arrays (left panel), and sweeps in the y direction are only performed in the wide subgrid arrays (right panel). Thus every sweep operation occurs along the largest possible array dimension and involves only the arrays kept in the memory of the responsible processor. Communication between processors is only required to exchange updated block arrays of the hydrodynamic variables after each sweep is finished.