A&A 459, 1-19 (2006)
DOI: 10.1051/0004-6361:20053898

Supersonic turbulence in shock-bound interaction zones

I. Symmetric settings

D. Folini¹ - R. Walder^2,3

1 - Institut für Astronomie, ETH Zürich, 8092 Zürich, Switzerland
2 - Observatoire de Strasbourg, 67000 Strasbourg, France
3 - Max-Planck-Institut für Astrophysik, 85741 Garching, Germany

Received 24 July 2005 / Accepted 26 June 2006

Abstract
Colliding hypersonic flows play a decisive role in many astrophysical objects. They contribute, for example, to the molecular cloud structure, the X-ray emission of O-stars, differentiation of galactic sheets, appearance of wind-driven structures, or, possibly, to the prompt emission of $\gamma$-ray bursts. Our intention is thorough investigation of the turbulent interaction zone of such flows, the cold dense layer (CDL). In this paper, we focus on the idealized model of a 2D plane parallel isothermal slab and on symmetric settings, where both flows have equal parameters. We performed a set of high-resolution simulations with upwind Mach-numbers, $5 < M_{{\rm u}} < 90$. We find that the CDL is irregularly shaped and has a patchy and filamentary interior. The size of these structures increases with  $\ell _{{\rm cdl}}$, the extension of the CDL. On average, but not at each moment, the solution is nearly self-similar and only depends on $M_{{\rm u}}$. We give the corresponding analytical expressions, with numerical constants derived from the simulation results. In particular, we find the root-mean-square Mach-number to scale as $M_{{\rm rms}} \approx 0.2~M_{{\rm u}}$. The mean density, $\rho_{{\rm m}} \approx 30~\rho_{{\rm u}}$ is independent of  $M_{{\rm u}}$. The fraction $f_{{\rm eff}}$ of the upwind kinetic energy that survives shock passage scales as $f_{{\rm eff}}= 1 - M_{{\rm rms}}^{-0.6}$. This dependence persists if the upwind flow parameters differ from one side to the other of the CDL, indicating that the turbulence within the CDL and its driving are mutually coupled. Another finding points in the same direction, namely that the auto-correlation length of the confining shocks and the characteristic length scale of the turbulence within the CDL are proportional. Larger upstream Mach-numbers lead to a faster expanding CDL, confining interfaces that are less inclined with respect to the upstream flow direction, more efficient driving, and finer interior structure with respect to the extension of the CDL.

Key words: shock waves - instabilities - turbulence - hydrodynamics - ISM: kinematics and dynamics - stars: winds, outflows

   
1 Introduction

Supersonically turbulent, shock-bound interaction zones are important for a variety of astrophysical objects. They contribute, for example, to structure formation in molecular clouds (Ballesteros-Paredes et al. 1999a; Heyer & Brunt 2004; Vázquez-Semadeni 2004; Hartmann et al. 2001; Hueckstaedt 2003; Hunter et al. 1986) and to galaxy formation (Kang et al. 2005; Anninos & Norman 1996). They affect the X-ray emission of line-driven hot-star winds (Feldmeier & Owocki 1998; Feldmeier et al. 1997; Oskinova et al. 2004; Owocki et al. 1988) and contribute substantially to the physics and emitted spectrum of colliding wind binaries (Stevens et al. 1992; Marchenko et al. 2003; Corcoran et al. 2005; Folini & Walder 2000; Nussbaumer & Walder 1993). The currently most promising model for the prompt emission of $\gamma$-ray bursts is based on internal shocks (Fan & Wei 2004; Piran 2004; Panaitescu et al. 1999; Rees & Meszaros 1994). A similar mechanism has been proposed for micro-quasars (Kaiser et al. 2000), BL Lacs and Blazars (Mimica et al. 2004; Ghisellini et al. 2002), and Herbig-Haro objects (Matzner & McKee 1999).

So far, the shape and turbulent interior of shock-bound interaction zones have been mostly studied separately. In this paper we focus on the system as a whole, stressing that upwind flows, confining interfaces of the interaction zone, and the interior structure of this zone form a tightly coupled system. The turbulence within the interaction zone affects the shape of the confining shocks, which in turn determines how much energy is thermalized at these shocks and how much energy remains available for driving the turbulence.

A variety of papers have been written on the shape and stability of 2D interaction zones, of which we mention only a few. Vishniac (1994) shows by analytical means that geometrically thin, isothermal, 2D, planar, shock-bounded slabs are non-linearly unstable, coining the term non-linear thin shell instability, or NTSI, for this instability. Blondin & Marks (1996) essentially reproduce these analytical predictions numerically, also mentioning the occurrence of supersonic turbulence within the slab. Performing 2D radiative and isothermal simulations of colliding molecular clouds, Klein et al. (1998) observe the complex shaping and instability of the collision zone. The role of a radiative cooling layer has been addressed by several authors. Strickland & Blondin (1995) numerically investigated flows against a wall in 2D, finding that an unstable cooling layer introduces disturbances in the interface separating the cooling layer from the cooled matter. Looking at colliding flows instead of a flow against a wall, Walder & Folini (1998) show that one unstable cooling layer is sufficient to destabilize both confining interfaces of the cooled matter. In addition, the cooled matter becomes supersonically turbulent. If self-gravity is included fragmentation of the interaction zone is observed (Anninos & Norman 1996; Hunter et al. 1986).

An overwhelming amount of literature meanwhile exists on supersonic turbulence. At least part of this attention arises because it is thought that supersonic turbulence can explain the structuring and support of molecular clouds and thus that it plays a decisive role in star formation. A comprehensive view of this issue can be found in the recent reviews by Mac Low & Klessen (2004), Elmegreen & Scalo (2004) and Scalo & Elmegreen (2004). Of particular interest for the work we present here is the paper by Mac Low (1999), where Fig. 4 shows that the wave length of the driving is apparent in the spatial scale of the turbulent structure for monochromatically driven turbulence in a 3D periodic box. The possible importance of the finite size of the slab was recently pointed out by Burkert & Hartmann (2004).

We are trying to make four points with this paper. First, we argue that, within the frame of isothermal Euler equations and in infinite space, the solution may be self-similar and dependent only on the upstream Mach-number, at least to first approximation. Based on this assumption, we give expressions for average quantities of the slab. Second, we show that the numerical solution, which is defined only on a finite computational domain and includes (implicit) numerical dissipation, remains close to self-similar, as long as the width of the slab is small and the root-mean-square Mach-number larger than one. Third, we stress the tight mutual coupling between the turbulence and its driving. Fourth, we point out that spatial scales generally grow with extension $\ell _{{\rm cdl}}$ of the interaction zone, but decrease with increasing upstream Mach-number $M_{{\rm u}}$.

Results are based on a set of simulations that differ only in their upwind Mach-numbers. In this paper we restrict the analysis of these simulations to the above-mentioned three objectives. We postpone a more detailed analysis of the interior structure of the interaction zone to a subsequent paper.

In the following, we first give the details of our physical model and numerical method in Sect. 2. In Sect. 3 we derive the self-similar scaling relations. The numerical results are present in Sect. 4. Discussion follows in Sect. 5, and conclusions in Sect. 6.

   
2 Physical model and numerical method

The numerical treatment of supersonic turbulence is an issue in its own right, so we start this section with a brief summary of some results that are relevant to the present work. We then specify the physical model we consider, explain the numerical method we use and the simulations we perform.

   
2.1 Simulating supersonic turbulence

The shock-compressed layer studied in this paper is supersonically turbulent with root-mean-square Mach-numbers between about 1 and 10. An important fraction of the kinetic energy is dissipated in shocks. Euler equations are sufficient for describing this part of the problem. A cascade transfers the remaining energy to higher and higher wave numbers until it is finally destroyed on the viscous dissipation scale. To also capture this part of the problem, the compressible Navier-Stokes equations should be used; however, the range of spatial scales associated with the energy cascade exceeds the capacity of any computer by far.

In subsonic turbulence, one way out is to use a suitable sub-grid scale model. The model is used to compute an effective viscosity coefficient, which should mimic the cascading between the smallest scale still resolved by the numerical grid and the viscous dissipation scale as precisely as possible. This coefficient is then used in the Navier-Stokes equations instead of the physical viscosity (Lesieur 1999). For the approach to work it is essential that the effective viscosity obtained from the sub-grid scale model exceeds the (implicit) numerical viscosity of the overall numerical scheme. This can be achieved in subsonic turbulence by the use of low-dissipation schemes (Lele 1992).

In supersonic turbulence, explicit sub-grid scale modeling so far does not exist in the above sense. The basic reason is that the numerical treatment of supersonic turbulence requires schemes that can treat shocks appropriately, such as the widely used shock capturing schemes. The (implicit) numerical viscosity of such schemes is, however, much too large to match the above requirement, even if the schemes are of a high order (Garnier et al. 1999; Porter et al. 1992). One strategy for this case, the so called MILES approach (monotone integrated large-eddy simulation), was proposed by Boris et al. (1992) and further explored by Porter et al. (1994,1992). The basic claim is that the numerical viscosity inherent to shock capturing schemes (LeVeque 2002; Hirsch 1995) acts already as a physically correct sub-grid scale model. Solving the Euler equations by means of a shock capturing scheme thus should yield the correct physical answer.

The validity of the claim that implicit numerical viscosity alone leads to a correct physical solution was investigated by Garnier et al. (1999) for a selection of shock capturing schemes, among them a MUSCL-scheme (monotone upwind scheme for conservation laws) similar to the one we use (see Sect. 2.3). For the cases considered (essentially decaying subsonic), they find that the scheme indeed acts as a (very dissipative) sub-grid scale model in that it preserves the flow from energy accumulation on small spatial scales. However, they also find that structures defined on less than 5 grid points are affected by substantial numerical damping. Porter et al. (1994) find, in addition, that the dissipation properties of their scheme (MUSCL with PPM) are highly non-linear, and also they depend not only on the grid spacing but also on the wave length of the flow structure. Structures on less than 32 grid points are affected by numerical damping.

We rely on the MILES approach in this paper for the lack of a better model, although, to our knowledge, the validity and quality of the approach has never been tested for supersonic turbulence. The numerical solutions we obtain are thus rather solutions of the Navier-Stokes equations. Nevertheless, as dissipation in shocks by far dominates numerical dissipation, we expect the "Euler character'' of the solution to prevail.

   
2.2 The model problem

The model problem we consider consists of a 2D, plane-parallel, infinitely extended, isothermal, shock compressed slab. A sketch is given in Fig. 1. Two high Mach-number flows, oriented parallel (left flow, subscript l) and anti-parallel (right flow, subscript r) to the x-direction, collide head on. The resulting high-density interaction zone, the shock compressed slab, is oriented in the y-direction. We denote this interaction zone by CDL for "cold dense layer'' to remain consistent with notation used already in Walder & Folini (1998,1996). We investigated this system within the frame of Euler equations (but see also Sect. 2.1), together with a polytropic equation of state,

$$\frac{\partial\rho}{\partial t} + \vec{\nabla} \left( \rho \vec{v} \right)$$ = 0,                  (1)

$$\frac{\partial\rho \vec{v}}{\partial t} + \vec{\nabla} \left( \rho \vec{v} \otimes \vec{v} + \frac{p}{\mu} I \right)$$ = 0, (2)

$$\frac{\partial E}{\partial t} + \vec{\nabla} ( \vec{v} \left(E + p \right) )$$ = 0, (3)

$$p/(\gamma - 1) .$$ e = (4)

Here, $\rho$ is the particle density, $\mu$ the average mass per particle, $\vec{v} = (v_{x},v_{y})$ is the velocity vector, p thermal pressure, I the identity tensor, e the thermal energy density, and $E=\rho \vec{v}^{2}/2 + e$ the total energy density. For the polytropic exponent, we choose $\gamma = 1.000001$. This value guarantees that jump conditions and wave speeds of a Mach-90 shock are within 0.01 per cent of the isothermal values.

Within the frame of this paper we consider only symmetric settings, where the left (subscript l) and right (subscript r) colliding flow have identical parameters (subscript u for upstream): $\rho_{{\rm l}} = \rho_{{\rm r}} \equiv \rho_{{\rm u}}$ and $\vert v_{{\rm l}}\vert = \vert v_{{\rm r}}\vert \equiv v_{{\rm u}}$.

We look at the problem in a dimensionless form and express velocities in units of the isothermal sound speed $a=\sqrt{T k_{{\rm B}}/\mu}$, with T the temperature and $k_{{\rm B}}$ the Boltzmann constant. Densities we express in terms of the upstream density $\rho _{{\rm u}}$. Finally, we express lengths in units of ${Y}_{{\rm0}}$, the smallest y-extent of the computational domain we used. This artificial choice is necessary as there is no natural time-independent length scale to the problem (see Sect. 3).

Figure 1: Sketch of physical model problem. $\rho _{{i}}$, Mi, and si denote the density, Mach-number, and confining shock of the left ($i=\rm l$) and right ($i=\rm r$) flow. $\rho$ and M denote the density and Mach-number of the CDL. $\alpha$ is the absolute value of the angle between the x-axis and the tangent to the shock. CDL is the shock-compressed interaction zone. The dashed rectangle indicates the computational domain with y-extension Y. Periodic boundary conditions in y-direction imply periodic continuation of the solution (dotted continuation of left and right shock).

Open with DEXTER

   
2.3 Numerical method

Our results were obtained with the AMRCART-code. We used the multidimensional high-resolution finite-volume-integration scheme developed by Colella (1990) on the basis of a Cartesian mesh. Tests showed that this algorithm, compared to dimensional splitting schemes, is significantly more accurate in capturing flow features not aligned with the axis of the mesh. In all our simulations we used a version of the scheme that is (formally) second order accurate in space and in time for smooth flows.

We combine this integration scheme with the adaptive mesh algorithm by Berger (1985). While a rather coarse mesh was sufficient for the upwind flows, the turbulent CDL was resolved on a much finer scale.

We found it useful to have our CDL moving in positive x-direction at a speed of about Mach 20-40. If the CDL was essentially stationary with respect to the computational grid, we observed alignment effects of strong shocks that were nearly parallel to a cell interface (in y-direction). Through the global motion of the CDL, which implied supersonic motion of the confining shocks with respect to the computational grid, we got rid of this problem. We checked that this procedure introduced no systematic effects into the solution. The problem of alignment effects when dealing with high Mach-number flows, nearly stationary shocks, and high order upwind schemes is well known and not particular to our scheme (Jasak & Weller 1995; Colella & Woodward 1984; Quirk 1994). Other work arounds exist, such as smoothing of interfaces by additional viscosity, which is often applied in PPM implementations.

   
2.3.1 Numerical settings and integration time

In the x-direction, our computational domain extended over $200~ {Y}_{{\rm0}}$. The y-extent Y of our domain varied between simulations, ${Y}_{{\rm0}} \le {Y} \le 6~ {Y}_{{\rm0}}$ (see Table B.1). Boundary conditions at the left and right boundaries (x-direction) were "supersonic inflow''. In the y-direction we had periodic boundary conditions. The cell size at the coarsest level was $0.2 ~ {Y}_{{\rm0}}$. The cells at the finest level, covering the CDL, were smaller by a factor 2⁶ to 2⁹, yielding between 320 and 2560 cells over a distance  ${Y}_{{\rm0}}$ (depending on the simulation, see Table B.1).

As will be shown, the relevant time-dependent quantity for the evolution of CDL mean quantities is the average x-extension of the CDL, $\ell _{{\rm cdl}}$. We defined it as $\ell_{{\rm cdl}} \equiv V / {Y}$, where V is the 2D volume of the CDL. For later use we also introduce the volume integrated density $m_{{\rm cdl}} \equiv \int_{{\rm V}} \rho$, the mean density $\rho_{{\rm m}} \equiv m_{{\rm cdl}} / V$, and the average column density $N \equiv m_{{\rm cdl}}/{Y} = \rho_{{\rm m}} \ell_{{\rm cdl}}$. The last quantity was made dimensionless by division through $N_{{\rm0}} \equiv \rho_{{\rm u}}{Y}_{{\rm0}}$. We stopped most simulations at $\ell_{{\rm cdl}} = {Y}/2$.

   
2.3.2 Initial conditions

We investigated three different initial conditions, I=0, 1, 2.

I=0: No CDL exists at t=0. The left and right flows are initially separated by a single interface. The interface is wiggled with a single, sinusoidal mode of wave length 0.1  Y and amplitude  $0.0195~ {Y}_{{\rm0}}$ (about 3 to 25 grid cells, depending on the discretization).

I=1: A CDL is present at time t=0. It has a column density of $N = 14 ~ N_{{\rm0}}$ and a thickness of $0.03125~ {Y}_{{\rm0}}$. The confining shocks are both wiggled, with the same sinusoidal mode and amplitude as the interface in the case I=0. The mass within the CDL is at rest and of constant density, $\rho = \rho_{{\rm u}} M_{{\rm u}}^{2}$, the density the CDL would have in 1D. Note that this initialization implies some violation of the Rankine-Hugoniot jump conditions at the interfaces.

I=2: A CDL is present at time t=0, with column density $N = 56~ N_{{\rm0}}$ and a thickness of $0.125~ {Y}_{{\rm0}}$. The right shock is wiggled as for I=1, the left shock is straight. The density and velocity in the CDL are set as for I=1.

We stress that the initial wiggling of the shocks is not compelling. The only effect of this wiggling is to speed up the initial phase of the evolution. Test cases using another wiggling or starting from straight shocks end up like the simulations we are going to present in the following.

We would like to add a side note on this last point, from our observation that the slab is also destabilized when bound by straight shocks. This has already been reported by Blondin & Marks (1996), who ascribed the destabilization to "numerical noise''. Meanwhile, Robinet et al. (2000) have investigated what is called the carbuncle phenomenon in some more detail. They showed that - contrary to what has been believed so far - a single straight shock is linearly unstable for exactly one mode associated to the upstream Mach-number of $M_{{\rm crit}} = [(5+\gamma) / (3-\gamma)]^{1/2}$. For isothermal conditions, this yields $M_{{\rm crit}}= \sqrt{3}$. They also showed that this single unstable mode is sufficient for making straight shocks aligned with the mesh numerically unstable at all Mach-numbers if the computation is done with a low-viscosity, high-order, shock-capturing scheme. To what degree this instability for a straight shock of any Mach-number is really physical seems an open question to us.

2.4 The different runs

The runs we performed differ in their upwind Mach-numbers, which lie in a range $5 \la M_{{\rm u}} \la 90$, as well as in their initialization, numerical discretization, and the y-extent of the domain. The labels of the different runs are built up as $M\_I.R.Y$. Here, M is the upwind Mach-number, I the initialization (0, 1, or 2), R gives the refinement of the spatial discretization, relative to the coarsest grid simulation we performed (1, 2, 4, or 8). R=1 corresponds to a finest cell size of about $3 \times 10^{-3} ~ {Y}_{{\rm0}}$, R=2 indicates a twice smaller cell size. Y is the domain size (1, 2, 4, or 6) in units of  ${Y}_{{\rm0}}$. For example, R22_0.2.4 denotes a run with $M_{{\rm u}}=22$, initialization I=0, finest cell size about $1.5 \times 10^{-3} {Y}_{{\rm0}}$, and y-extent $4 ~ {Y}_{{\rm0}}$.

The runs we performed are listed in Table B.1. Individual columns in Table B.1 contain (column number in square brackets): label of run [1], following the scheme label=  $M_{{\rm u}}\_~I.R.Y$, where I is the initial condition, R the refinement factor such that cell size = $3.125 \times 10^{-3}~{Y}_{{\rm0}} / {R}$, and Y is the y-extension of the computational domain in units of  ${Y}_{{\rm0}}$; Mach-number of upstream flow, $M_{{\rm u}}$ [2]; stopping time of simulation in terms of $\ell ( N )$ [3]; y-averaged x-extension of CDL at stopping time, relative to y-extent of computational domain, $\ell_{{\rm cdl}}/{Y}$ [4]; average quantities [5-8] of: rms Mach-number, $M_{{\rm rms}}$ [5]; mean density in units of upstream density, $\rho _{{\rm m}}/\rho _{{\rm u}}$ [6]; shock length in units of y-domain, $\ell_{{\rm sh}}/Y$ [7]; driving efficiency, $f_{{\rm eff}}$ [8]; averages taken over $10 \le \ell ( N ) \le 70$ for I=0 and over $60 \le \ell( N ) \le 120$ for I=1, for I=2 we give the values at the end of the simulation in parentheses instead.

   
3 Scaling properties of the model problem

Figure 2: The self-similar 1D solution of isothermal colliding supersonic flows in density ( top) and velocity ( bottom). The interaction zone (labeled CDL) is bounded by two shocks, $s_{{\rm l}}$ and  $s_{{\rm r}}$, having speeds $v^s_{\rm l}$ and  $v^s_{\rm r}$ in the rest frame of the CDL. The density and velocity of the 1D interaction zone, we denote by $\rho _{{\rm 1d}}$ and $v_{{\rm 1d}}$, respectively.

Open with DEXTER

Within the frame of Euler equations and in infinite space, the problem of isothermal supersonically colliding flows can be solved analytically in 1D. The solution, sketched in Fig. 2 and Sect. 3.1, is self-similar and depends only on two free parameters, the Mach-numbers of the left and right upwind flow. In 2D the situation is more complicated: the solution is unstable (Blondin & Marks 1996; Vishniac 1994), the shocks confining the CDL are non-stationary and oblique, the interior of the CDL is supersonically turbulent.

Nevertheless, in infinite space it seems reasonable to assume that the solution, on average, may still evolve in a self-similar manner. We base this assumption on the following two observations. First, the isothermal Euler equations are scale-free in infinite space. Second, the free parameters of the problem ( $\rho _{{\rm u}}$, $M_{{\rm u}}$, and a) do not introduce any fixed length or time scale. Under these conditions, it is possible that the solution also does not depend on length or time separately, but only on their ratio. If so, all length scales should evolve equally with time, which implies, in particular, that the solution then should not depend on the extension of the CDL. We stress, however, that we have no proof of the above assumption of self-similarity.

In the remainder of this section, we elaborate a bit further on the implications of the assumed self-similarity. In Sect. 4 we will see that the relations derived here give a good approximation of the numerical results, but we stress already here three important points. The numerical simulations are carried out in finite space (not infinite); numerical dissipation might play a role; and the simulations are stopped for the most part while the CDL is still small, about half the size of the y-extent of the computational domain. Important aspects that can only be obtained from the numerical solution include quantities related to the driving of the turbulence, the values of proportionality constants, and the interior structure of the CDL. We neglect this last aspect, however, in the current paper to focus on mean quantities instead.

   
3.1 Self-similar 1D solution

Denoting the density and velocity of the CDL by $\rho _{{\rm 1d}}$and $v_{{\rm 1d}}$, and those of the left and right upwind flows by $\rho _{{i}}$ and v_i ($i=\rm l,r$), the solution in the rest frame of the CDL is given by

   

$$\rho_{{\rm 1d}} / \rho_{{i}} M_{{i}}^{\rm 2} + 1 \approx M_{{i}}^2,$$ = (5)

$$v_{{\rm 1d}}$$ = 0, (6)

$$aM_{{i}} / (M_{{i}}^{2} - 1) \approx a/ M_{{i}} \ll a.$$ |v^s_i| = (7)

Here, v^s_i is the velocity of the confining shocks and a is again the isothermal sound speed. The approximations hold for large Mach-numbers. The self-similar character is apparent: the solution is not a function of x and t but only a function of x/t through  v^s_i.

A relation between characteristic length and time scales of the solution, the self-similarity variable $\kappa_{{\rm 1d}}$, can be obtained as follows. As a length scale, we take the spatial extension  $\ell_{{\rm 1d}}$ of the CDL, and as a time scale the time $\tau$ needed to accumulate the corresponding column density  ${ N }_{{\rm 1d}}$. From the relations

$$N_{{\rm 1d}} = \rho_{{\rm 1d}} \ell_{{\rm 1d}}.$$ (8)

and

$$N_{{\rm 1d}} = \tau \left( \rho_{\rm l} v_{\rm l} + \rho_{\rm r} v_{\rm r}\right)$$ (9)

and using $\rho_{\rm l}/\rho_{\rm r} = M^{2}_{\rm r}/M^{2}_{\rm l}$ (see Eq. (5)), we obtain

$$\kappa_{{\rm 1d}} \equiv \frac{\ell_{{\rm 1d}}}{\tau} = a \frac{M_{\rm l} + M_{\rm r}}{M_{\rm l} \cdot M_{\rm r}}\cdot$$ (10)

Thus for strong shocks $\kappa_{{\rm 1d}}$ is nothing else than $\vert v^{s}_{{\rm l}}\vert + \vert v^{s}_{{\rm r}}\vert$. Specializing to symmetric settings ( ${\rm l} = {\rm r}$) yields $\rho_{{\rm 1d}} / \rho_{{\rm u}} = M_{{\rm u}}^{2}$ and $\kappa_{{\rm 1d}} =2a/M_{{\rm u}}$.

   
3.2 Scaling properties of the 2D symmetric solution

In the following, we derive scaling relations for the 2D solution, assuming self-similarity. We confront these relations with corresponding numerical results in Sect. 4.

3.2.1 Density, Mach-number, self-similarity variable

In the following, all velocities are again given in the rest frame of the CDL and we assume that a self-similar solution exists. A natural choice for the (constant) self-similarity variable then is again $\kappa_{{\rm 2d}} \equiv \ell_{{\rm cdl}}/\tau$. Using the definitions of Sect. 2.3.1 we must have, as in the 1D case,

  

$$N = \rho_{{\rm m}} \ell_{{\rm cdl}},$$ (11)

$$N = 2 \tau \rho_{{\rm u}} v_{{\rm u}}.$$ (12)

Dividing the two equations through each other yields $\kappa_{{\rm 2d}} = 2 \rho_{{\rm u}} v_{{\rm u}} / \rho_{{\rm m}}$. As $\kappa_{{\rm 2d}}$ is a constant, the CDL mean density  $\rho_{{\rm m}}$must be constant in time. The root-mean-square velocity  $v_{{\rm rms}}^{2}$ then has to be constant in time as well, at least if the CDL density and velocity, $\rho$ and v, are uncorrelated (in which case we can replace the average over the product  $\rho v^{2}$ by the product of the averages of $\rho$ and v²) and if kinetic pressure dominates over thermal pressure. This can be seen from equating the total upwind pressure with the total pressure within the CDL,

$$\rho_{{\rm u}} (a^{2} + v_{{\rm u}}^{2}) = \rho_{{\rm m}} (a^2 + v_{{\rm rms}}^2).$$ (13)

The simplest ansatz for $\rho_{{\rm m}}$ and $v_{{\rm rms}}$is that they only depend on the upstream Mach-number,

  

$$\rho_{{\rm m}}/\rho_{{\rm u}} = \eta_{{\rm 1}} M_{{\rm u}}^{\beta_{{\rm 1}}},$$ (14)

$$v_{{\rm rms}}/a = \eta_{{\rm 2}} M_{{\rm u}}^{\beta_{{\rm 2}}}.$$ (15)

Using the ansatz for $\rho_{{\rm m}}$ we obtain a first expression for $\kappa_{{\rm 2d}}$ from Eqs. (11) and (12),

$$\kappa_{{\rm 2d}} = 2 a \eta_{{\rm 1}}^{-1} M_{{\rm u}}^{1-\beta_{{\rm 1}}} \propto a M_{{\rm u}}^{1-\beta_{{\rm 1}}}.$$ (16)

A second expression for $\kappa_{{\rm 2d}}$, we obtain from Eq. (13)

$$\rho_{{\rm u}}a^{2}(1+M_{{\rm u}}^{2}) = \rho_{{\rm m}} (a^2 ... ... + \eta_{{\rm 2}}^{2} M_{{\rm u}}^{2 \beta_{{\rm 2}}} \right).$$ (17)

Again using Eq. (12) to replace N, one obtains

$$\kappa_{{\rm 2d}} = 2 a M_{{\rm u}} \frac{ 1 + \eta_{{\rm 2... ...ta_{{\rm 2}} -1} \propto a M_{{\rm u}}^{2\beta_{{\rm 2}} -1}.$$ (18)

The approximation is good for high Mach-number flows, with $\eta_{{\rm 2}}^{2} M_{{\rm u}}^{2 \beta_{{\rm 2}}} \gg 1$, and for $\beta_{{\rm 2}} > 0$, which is, however, to be expected for supersonic turbulence. Comparing Eqs. (16) and (18) gives

  

$$\beta_{{\rm 2}} = 1 - \beta_{{\rm 1}}/2,$$ (19)

$$\eta_{{\rm 1}}^{-1} = \eta_{{\rm 2}}^{2}.$$ (20)

   
3.2.2 Driving energy

From energy conservation, we have $\dot{\cal E}_{{\rm diss}} = \dot{\cal E}_{{\rm drv}} - \dot{\cal E}_{{\rm kin}}$. Here  $\dot{\cal E}_{{\rm drv}}$ is the energy flux density entering the CDL per time and per unit length in the y-direction, and $\dot{\cal E}_{{\rm diss}}$denotes the energy density dissipated per time within an average column of length  $\ell _{{\rm cdl}}$ of the CDL. Finally, $\dot{\cal E}_{{\rm kin}}$ is the change per time of the kinetic energy contained within such an average column. We first turn to the driving energy $\dot{\cal E}_{{\rm drv}}$ and come back to $\dot{\cal E}_{{\rm diss}}$ and $\dot{\cal E}_{{\rm kin}}$ in Sect. 3.2.3.

Part of the total (left plus right) upwind kinetic energy flux density, ${\cal F}_{{\rm e_{{\rm kin}},u}} = \rho _{{\rm u}} v_{{\rm u}}^{3}$, is thermalized at the shocks confining the CDL. The remaining part, $\dot{\cal E}_{{\rm drv}}$, drives the turbulence in the CDL. We assume that $\dot{\cal E}_{{\rm drv}}$ and ${\cal F}_{{\rm e_{{\rm kin}},u}}$ are related by a function of the upwind Mach-number only,

$$\dot{\cal E}_{{\rm drv}} = f_{{\rm eff}}(M_{{\rm u}}){\cal F}_{{\rm e_{{\rm kin}},u}}.$$ (21)

We call the function $f_{{\rm eff}}$ the driving efficiency. An expression for $f_{{\rm eff}}$ can be derived by using the jump conditions for strong, oblique shocks,

 

$$\rho_{{\rm d}} = \rho_{{\rm u}} M_{{\rm\perp,u}}^{2} = \rho_{{\rm u}} M_{{\rm u}}^{2} \sin^{2}\alpha,$$

$$v_{{\rm\perp,d}} = v_{{\rm\perp,u}} M_{{\rm\perp,u}}^{-2} = \frac{a}{M_{{\rm u}} \sin \alpha},$$

$$v_{{\rm\parallel,d}} = v_{{\rm\parallel,u}} = a M_{{\rm u}} \cos\alpha.$$ (22)

The subscript d denotes downstream quantities, right after shock passage; the subscripts ${{\rm\perp}}$ and ${{\rm\parallel}}$denote flow components perpendicular and parallel to the shock, respectively; and $\alpha$ is given in Fig. 1. Using Eq. (22) we obtain

 

$$\dot{\cal E}_{{\rm drv}} \frac{1}{{Y}}\int_{s_{{\rm l,r}}} {\rm d}s \frac{\rho_{{\rm d}} v_{{\rm d}}^{2}}{2} v_{{\rm\perp,d}}$$ =

$$\frac{\rho_{{\rm u}}v_{{\rm u}}^{3}}{2Y} \int_{{Y}_{{\rm l,r}}} {\rm d}y \left(1 - \sin^{2}\alpha + \frac{1}{M_{{\rm u}}^{4}\sin^{2}\alpha}\right),$$ = (23)

where the integral over $s_{{\rm l,r}}$ and ${Y}_{{\rm l,r}}$runs over both shocks and where it was used that $\sin \alpha \; {\rm d}s = {\rm d}y$. The last term on the right hand side of Eq. (23) is omitted in the following. This is justified, as the shocks we observe in our simulations fulfill $\sin \alpha \gg M_{{\rm u}}^{-2}$ for the most part (see Sect. 4.2.2). For $f_{{\rm eff}}(M_{{\rm u}})$ we thus obtain

$$f_{{\rm eff}} = \frac{1}{2Y} \int_{{Y_{{\rm l,r}}}} {\rm d}y (1 - \sin^{2}\alpha) \equiv 1 - \sin^{2}\alpha_{{\rm eff}}$$ (24)

where we used the midpoint rule. The angle $\alpha_{{\rm eff}}$ can be interpreted as an average bending angle. As the ansatz for the Mach-number dependence of $f_{{\rm eff}}$ we thus take

$$f_{{\rm eff}} = 1 - \sin^{2}\alpha_{{\rm eff}} = 1 - \eta_{{\rm 3}}M_{{\rm u}}^{\beta_{{\rm 3}}}.$$ (25)

   
3.2.3 Energy dissipation

A first expression for the column-integrated dissipated energy per time can be obtained from energy conservation, $\dot{\cal E}_{{\rm diss}} = \dot{\cal E}_{{\rm drv}} - \dot{\cal E}_{{\rm kin}}$. For $\dot{\cal E}_{{\rm drv}}$ we just derived an expression, Eqs. (21) and (25). For $\dot{\cal E}_{{\rm kin}}$ we obtain, within the frame of self-similarity,

$$\dot{\cal E}_{{\rm kin}} = \frac{\rho_{{\rm m}} v_{{\rm rms}}... ...u}} a^{3} \eta_{{\rm 2}}^{2} M_{{\rm u}}^{3-\beta_{{\rm 1}}},$$ (26)

where we used Eqs. (14), (15), and (18) to (20). Together we get

$$\dot{\cal E}_{{\rm diss}} = \rho_{{\rm u}} a^{3} M_{{\rm u}}^... ...eta_{{\rm 2}}^{2} M_{{\rm u}}^{{\rm -\beta_{{\rm 1}}}}\right].$$ (27)

The energy dissipated per time within an average column of length  $\ell _{{\rm cdl}}$ is thus independent of this length. If energy dissipation occurs only (as within the frame of Euler equations) or at least dominantly in shocks, this implies that the average distance between shocks increases and / or the average strength of the shocks decreases as the CDL grows.

A second expression for $\dot{\cal E}_{{\rm diss}}$ can be obtained from dimensional considerations. The energy dissipated per unit volume per unit time must be proportional to $\rho_{{\rm diss}} v_{{\rm diss}}^{3} \ell_{{\rm diss}}^{-1}$. Here, $\rho_{{\rm diss}}$, $v_{{\rm diss}}$, and $\ell_{{\rm diss}}$ are the characteristic density, velocity, and length scale of the dissipation. The energy dissipation within an average column of length $\ell _{{\rm cdl}}$can thus be written as $\dot{\cal E}_{{\rm diss}} \propto \rho_{{\rm diss}} v_{{\rm diss}}^{3} \ell_{{\rm diss}}^{-1} \ell_{{\rm cdl}}$. As all length scales must evolve equally with time within the frame of self-similarity, $\ell_{{\rm cdl}}/\ell_{{\rm diss}}$ must be constant, thus

$$\dot{\cal E}_{{\rm diss}} \propto \rho_{{\rm diss}} v_{{\rm diss}}^{3}.$$ (28)

Comparison of Eqs. (27) and (28) suggests $v_{{\rm diss}} \propto a M_{{\rm u}}$ and a more complicated Mach-number dependence for $\rho_{{\rm diss}}$. As $v_{{\rm rms}}$ is the only velocity scale we have, it seems natural to assume that $v_{{\rm diss}} \propto v_{{\rm rms}}$. It then follows that $v_{{\rm rms}} \propto a M_{{\rm u}}$ or $\beta_{{\rm 2}}=1$ (and $\beta_{{\rm 1}}=0$). We note that Gammie & Ostriker (1996) even found $v_{{\rm diss}} = v_{{\rm rms}}$ for a 1D case.

   
3.3 Summary of expected scaling relations

If a self-similar solution exists, we expect the following dependencies:

      

$$\rho_{{\rm m}} \eta_{{\rm 1}} \rho_{{\rm u}} M_{{\rm u}}^{\beta_{{\rm 1}}} = \eta_{{\rm 1}} \rho_{{\rm u}},$$ = (29)

$$M_{{\rm rms}} \eta_{{\rm 2}} M_{{\rm u}}^{\beta_{{\rm 2}}} = \eta_{{\rm 1}}^{-1/2} M_{{\rm u}},$$ = (30)

$$\kappa_{{\rm 2d}} \ell_{{\rm cdl}}/\tau = 2 \eta_{{\rm 1}}^{-1} a M_{{\rm u}},$$ = (31)

$$\dot{\cal E}_{{\rm drv}} \rho_{{\rm u}} a^{3} M_{{\rm u}}^{3} \left(1 - \eta_{{\rm 3}} M_{{\rm u}}^{\beta_{{\rm 3}}}\right),$$ = (32)

$$\dot{\cal E}_{{\rm kin}} \rho_{{\rm u}} a^{3} M_{{\rm u}}^{3} \; \eta_{{\rm 2}}^{2},$$ = (33)

$$\dot{\cal E}_{{\rm diss}} \rho_{{\rm u}} a^{3} M_{{\rm u}}^{3} \left(1 -\eta_{{\rm 3}} M_{{\rm u}}^{\beta_{{\rm 3}}} - \eta_{{\rm 2}}^{2}\right).$$ = (34)

Note the differences to the 1D solution: Eq. (29) predicts the CDL mean density to be independent of $M_{{\rm u}}$ and $\kappa_{{\rm 2d}} \propto a M_{{\rm u}}$, in contrast to $\rho_{{\rm 1d}} \propto M_{{\rm u}}^{2}$ and $\kappa_{{\rm 1d}} \propto a / M_{{\rm u}}$.

In deriving the above relations, we made four basic assumptions: a) we have simple Mach-number dependencies of $\rho_{{\rm m}}$, $v_{{\rm rms}}$, and $f_{{\rm eff}}$, Eqs. (14), (15), and (25); b) the CDL density and velocity are uncorrelated; c) we have high Mach-numbers in the sense that $\eta_{{\rm 2}}^{2} M_{{\rm u}}^{2 \beta_{{\rm 2}}} \gg 1$ or $M_{{\rm rms}}^{2} \gg 1$; d) $v_{{\rm diss}} \propto v_{{\rm rms}}$.

In Sect. 4 we are going to check the validity of these assumptions and confront Eqs. (29) to (34) with numerically obtained values. We expect good agreement as long as $M_{{\rm rms}} \gg 1$, thus dissipation in shocks likely dominates, and as long as $\ell_{{\rm cdl}} \ll {Y}$. The "Euler character'' of the solution should prevail under these conditions. We also determine those quantities that cannot be derived analytically. These are, on the one hand, the coefficients $\eta_{{\rm 1}}$ and $\eta_{{\rm 3}}$, as well as the exponent $\beta_{{\rm 3}}$. On the other hand, there are quantities for which we have no analytical expression at all, like the wiggling of the confining shocks, the associated distribution of the angle $\alpha$, or the Mach-number dependence of the length of the confining shocks.

   
4 Numerical results

We now present our numerical results. After a brief phenomenological description of the solution in Sect. 4.1, we give quantitative results for initial conditions I=0 in Sect. 4.2. Results for initial conditions I=1 and I=2 are given in Sect. 4.3, and asymmetric settings are briefly addressed in Sect. 4.4. Discretization and domain studies are the topic of Sect. 4.5.

   
4.1 Brief phenomenological description

We begin with a brief qualitative description of the CDL. As an example, the density structure of run R22_1.2.2 is shown in Fig. 3 for three different times.

A first characteristic is the local bending of the confining shocks. The spatial scale of these wiggles increases linearly with time, as the CDL accumulates more and more matter and gets more and more extended. The inclination of the wiggles with respect to the direction of the upstream flows decreases with increasing upstream Mach-number (see Sect. 4.2.2). Occasionally, we observe a superimposed "bending mode'' (e.g. bottom panel in Fig. 3), which in appearance is somewhat similar to the bending modes of the NTSI described by Vishniac (1994).

A second characteristic is the patchy appearance of the CDL. The turbulent interior is organized in filaments and patches, regions within which a flow variable remains more or less constant. The spatial extension of these patches increases as well as the CDL accumulates more and more matter. The flow variables clearly mirror the supersonic character of the turbulence: the contrast between high-density filaments and extended patches in Fig. 3 easily reaches two orders of magnitude, the root-mean-square velocity is well above sound, and the mean density is substantially reduced compared to the 1D case. Shocks within the CDL are ubiquitous.

   
4.2 Settings without CDL at t = 0

For symmetric settings, and if there is no CDL at time t=0, we expect to see the self-similar relations we derived in Sect. 3.2. We express the time evolution of the solution in terms of

$$\ell( N ) \equiv N / N_{{\rm0}} = \frac{\rho_{{\rm m}} \ell_{{\rm cdl}}} {\rho_{{\rm u}} {Y}_{{\rm0}}}\cdot$$ (35)

This function monotonically increases at about the same rate as the mean extension of the CDL, since $\rho_{{\rm m}} \approx \eta_{{\rm 1}} \rho_{{\rm u}}$ (Eq. (29)). In fact, $\rho_{{\rm m}} \approx 30 \rho_{{\rm u}}$(Sect. 4.2.1) and thus $\ell ( N ) = 60$ corresponds to $\ell_{{\rm cdl}} \approx 2 {Y}_{{\rm0}}$. For the symmetric case we consider in this paper, $\ell ( N )$ is proportional to the elapsed time. Using Eq. (12) to express N, we can write

$$\ell ( N ) \equiv N / N_{{\rm0}} = \frac{2 \tau \rho_{{\rm u... ...rm u}} {Y}_{{\rm0}}} = \tau \frac{2 v_{{\rm u}}}{Y_{{\rm0}}},$$ (36)

and $\ell ( N ) = 60$ then corresponds to a time $\tau = 30 Y_{{\rm0}} / v_{{\rm u}}$. Or, if we use $v_{{\rm u}} \approx 5 v_{{\rm rms}}$(Sect. 4.2.1) and ${Y}_{{\rm0}} \approx \ell_{{\rm cdl}} / 2$ for $\ell ( N ) = 60$, we obtain $\tau \approx 3 \ell_{{\rm cdl}} / v_{{\rm rms}}$.

Unless otherwise stated, averages and best fits in this section are always taken over the interval $10 \le \ell ( N ) \le 70$ and over all runs without CDL at time t=0. The interval was chosen such that initialization effects have died away and that domain effects do not matter yet (Sect. 4.5.1).

	Figure 3: The interaction zone of run R22_1.2.2, shown in density (logarithmic scale, in units of $\rho _{{\rm u}}$, color bar from 0 to 4), for three different times: $\ell ( N ) \approx 34$ ( top), $\ell ( N ) \approx 54$ ( middle), $\ell ( N ) \approx 74$ ( bottom). The spatial scale of patches, filaments, and wiggling of the confining shocks increases with $\ell ( N )$.
Open with DEXTER

We mention here already that the two most extreme simulations in terms of  $M_{{\rm u}}$, R5_0.2.4 and R87_0.2.4, often differ somewhat from the other simulations. In the case of R5_0.2.4, we ascribe the deviation to the only subsonic turbulence and the correlation of density and velocity ( $M_{{\rm rms}} \approx 0.9$ and corr $(\rho,v) \approx -0.4$, see Sect. 4.2.1). In the case of R87_0.2.4, the shocks become sometimes too strongly inclined with respect to the computational grid to be properly resolved by our numerical grid (Sect. 4.2.2).

   
4.2.1 CDL mean quantities and correlations

We first turn to the correlation of $\rho$ and v and the CDL mean quantities  $\rho_{{\rm m}}$ and  $M_{{\rm rms}}$, Eqs. (29) and (30). One of our basic assumptions in deriving these self-similar relations, namely point b) that the CDL density and velocity are uncorrelated, we find confirmed by our simulations. For nearly all symmetric simulations without initial CDL and for $10 \le \ell ( N ) \le 70$, we have $0.1 \ge {\rm corr}(\rho,v) \ge -0.1$. The only exceptions are the three low Mach-number runs R11_0.2.4, R11_0.2.2, and R5_0.2.4 with correlations of about -0.2, -0.2, and -0.4 respectively. The top panel of Fig. 4 shows the time evolution of corr$(\rho ,v)$ for five selected runs that differ only in their upwind Mach-number, $5 \le M_{{\rm u}} \le 90$.

In the middle and bottom panel of the the same figure, $\rho _{{\rm m}}/\rho _{{\rm u}}$ and $M_{{\rm rms}}/M_{\rm u}$are shown as a function of $\ell ( N )$ for the same runs. Two things are apparent. First, the ratios take similar values for all five runs, indicating that indeed $\beta_{{\rm 1}} \approx 0$ and $\beta_{{\rm 2}} \approx 1$for the exponents in Eqs. (29) and (30). Second, the ratios are not constant with $\ell ( N )$, indicating that the numerical solution is indeed only approximately self-similar. We come back to this point in Sect. 5.

To determine optimum exponents $\beta_{{i}}$, i=1,2, we rewrite Eqs. (29) and (30) as equations for  $\eta_{{\rm 1}}$ and  $\eta_{{\rm 2}}$ and minimize the variance  $\sigma^{2}(\eta_{{i}})$. Considering all data points within $10 \le \ell ( N ) \le 70$ of all runs without a CDL at t=0, we find the smallest variances for $\beta_{{\rm 1}}=0$ and for $\beta_{{\rm 2}}=1$. The corresponding means are $\mu(\eta_{{\rm 1}}) \approx 28$ and $\mu(\eta_{{\rm 2}}) \approx 0.21$. Although clearly identifiable, the minima of $\sigma$ are relatively shallow. Changing  $\beta_{{\rm 1}}$ or  $\beta_{{\rm 2}}$ by $\pm$0.1, or excluding the very low Mach-number case R5_0.2.4 (for which $M_{{\rm rms}} \approx 0.9$) changes $\sigma$ by only about 5%. By repeating the analysis but allowing for a linear dependence of  $\eta_{{i}}$ on $\ell ( N )$, we obtain the same optimum values for  $\beta_{{\rm 1}}$ and $\beta_{{\rm 2}}$ but with considerably smaller variance. As $\ell ( N )$ increases from 10 to 70, $\eta_{{\rm 1}}$ rises by about 25% (from 25 to 31), while $\eta_{{\rm 2}}$ decreases by about 15% (from 0.22 to 0.19).

Part of our assumption a), namely the simple Mach-number dependencies of $\rho_{{\rm m}}$ and $M_{{\rm rms}}$, thus seems justified. With $\eta _{{\rm 2}} = 0.21$, assumption c), $\eta_{{\rm 2}}^{2} M_{{\rm u}}^{2} \gg 1$, is also fulfilled for most of our simulations. An exception is again run R5_0.2.4, for which $\eta_{{\rm 2}}^{2} M_{{\rm u}}^{2} \approx 1$.

In summary, the simulation results, $\rho_{{\rm m}} \approx 28 \rho_{{\rm u}}$ and $M_{{\rm rms}} \approx 0.21 M_{{\rm u}}$, essentially confirm the expected relations, Eqs. (29) and (30). $\eta_{{\rm 1}}^{1/2} \eta_{{\rm 2}}=1$, predicted by Eq. (20), is fulfilled to within 10% at any given time. The mean density is (nearly) independent of  $M_{{\rm u}}$. As expected, the solution is only approximately self-similar, $M_{{\rm rms}}$ decreases by about 15% as $\ell ( N )$increases from 10 to 70.

Figure 4: Time evolution of corr$(\rho ,v)$ ( top), $\rho _{{\rm m}}/\rho _{{\rm u}}$ ( middle), and $M_{{\rm rms}}/M_{\rm u}$ ( bottom) for runs R5_0.2.4 (dotted, dark blue), R11_0.2.4 (dashed, purple), R22_0.2.4 (solid, red), R33_0.2.4 (dash-dotted, orange), R43_0.2.4 (dash-three-dots, green), and R87_0.2.4 (long dashes, pink). For these runs, $\ell ( N ) = 60$ corresponds to $\ell _{{\rm cdl}} \approx {Y}/2$.

Open with DEXTER

   
4.2.2 Confining shocks

The turbulence within the CDL is driven by the upstream flows. The confining shocks of the CDL affect this driving in two ways. The less inclined the shocks are on average with respect to the direction of the upstream flows (smaller angle $\alpha_{{\rm eff}}$ in Eq. (24)), the more kinetic energy survives shock passage and is available for driving the turbulence. The smaller the spatial scale on which the angle $\alpha$ varies, the smaller the scale on which the energy input changes. In the following, we analyze how these shock properties depend on  $M_{{\rm u}}$ and on  $\ell _{{\rm cdl}}$.

For this purpose, we specify the following basic quantities. The discrete x-position of the left and right shocks, $s_{{\rm l}}$ and $s_{{\rm r}}$, defined for each discrete y-position y_j as the two cell boundaries where the Mach-number drops for the first time from its upwind value M_u to 0.8 M_u. We determine the average extension of the CDL, $\ell _{{\rm cdl}}$, as

$$\ell_{{\rm cdl}} = \frac{1}{J} \sum_{j=1}^{J} [ s_{{\rm r}}(y_{j}) - s_{{\rm l}}(y_{j}) ].$$ (37)

The length of the left and right shocks, $\ell_{{\rm sh,l}}$ and $\ell_{{\rm sh,r}}$, are computed as

$$\ell_{{\rm sh},i} = \sum_{j=1}^{J} \left[ (s_{{i}}(y_{j}) - s_{{i}}(y_{j-1}))^{2} + (y_{j} - y_{j-1})^{2}\right]^{1/2},$$ (38)

where J is the number of cells in y-direction, and $i=\rm l,r$. We define the angle $\alpha_{{\rm l,r}}(y_{j})$ as the angle between the x-axis and the tangent to the shock (see Fig. 1). Its numerical computation is described in Appendix A. To obtain a number distribution, we sort the values $\alpha_{{\rm l,r}}(y_{j}) \in [0,\pi/2]$ into 60 bins. Finally, to obtain a measure for the scale on which the shocks are wiggled, we look at the auto-correlation functions $\Gamma_{{\rm l,r}}$,

$$\Gamma_{{i}}(y_{{\rm corr}}) = \frac{\langle[s_{{i}}(y_{j}) -... ...{\rm corr}}) - \bar{s}_{{i}}]\rangle} {\sigma^{2}_{{\rm s}}},$$ (39)

where $\sigma_{{\rm s}}^{2}$ is the variance of the shock position s_i, and  $\langle.\rangle$ denotes the mean over all discrete position y_j. For each time, we determine $y_{{\rm corr_{{\rm0}}}}$such that $\Gamma_{{i}} (y_{{\rm corr_{{\rm0}}}}) = 0.5$. Averaging $y_{{\rm corr_{{\rm0}}}}$ over both shocks gives a mean auto-correlation length $\ell _{{\rm corr}}$,

$$\ell_{{\rm corr}} = \frac{1}{2} \left[ y_{{\rm corr_{{\rm0... ..._{{\rm l}}) + y_{{\rm corr_{{\rm0}}}}(s_{{\rm r}}) \right].$$ (40)

A larger auto-correlation length $\ell _{{\rm corr}}$ then indicates that the shocks are wiggled on a larger spatial scale, but it does not give the scale of the wiggles in absolute units (see below).

Figure 5: Quantities related to the confining shocks: average extension $\ell _{{\rm cdl}}$ of the CDL ( first panel), total normalized shock length $l_{{\rm sh}}/ ({Y} M_{{\rm u}}^{0.8})$ ( second panel), number distribution (60 bins) of obliqueness angle $\alpha$ averaged over $10 \le \ell ( N ) \le 70$ ( third panel), auto-correlation length $\ell _{{\rm corr}} / {Y_{{\rm0}}}$ ( fourth panel), and scaled auto-correlation length $\ell _{{\rm corr}} / (\ell _{{\rm cdl}} M_{{\rm u}}^{-0.6})$, ( fifth panel). Individual curves denote the same runs as in Fig. 4.

Open with DEXTER

Figure 6: Variation of $\Gamma _{{\rm l}}$, color coded, as a function of $y_{{\rm corr}}$ for run R43_0.2.4 ( top panel). To allow for better display the color scale is limited to a range $-0.5 \le \Gamma _{{\rm l}} \le +0.5$. Lower or higher values of $\Gamma _{{\rm l}}$ are uniformly colored in dark blue or red, respectively. For the same run, $\Gamma _{{\rm l}}$ is shown as a function of $y_{{\rm corr}}$ for three selected times ( bottom panel). $\ell ( N ) = 30$ (solid), $\ell ( N ) = 50$ (dotted), $\ell ( N ) = 70$ (dashed).

Open with DEXTER

All four quantities, CDL extension, number distribution of angle $\alpha$, shock length, and correlation length, are shown in Fig. 5. The first panel of Fig. 5 shows the essentially linear growth of the CDL with $\ell ( N )$. The growth rate, however, slowly decreases with increasing $\ell ( N )$. The slope of a linear fit in the range $40 < \ell( N ) < 70$ is roughly 10% flatter than the slope obtained in the range $10 < \ell( N ) < 40$. This fits with the slight increase in $\rho_{{\rm m}}$, observable in the middle panel of Fig. 4. The second panel of Fig. 5 shows that the average shock length $\ell_{\rm sh} = 0.5(\ell_{\rm sh,l} + \ell_{\rm sh,r})$ is fairly constant with respect to $\ell ( N )$ but increases with $M_{{\rm u}}$. Assuming a dependence of the form $\ell_{\rm sh} = \eta_{{\rm sh}} {Y} M_{{\rm u}}^{\beta_{{\rm sh}}}$, the variance $\sigma^{2} (\eta_{{\rm sh}})$ becomes minimal for $\beta_{{\rm sh}}=0.8$. As can be seen, the two runs R5_0.2.4 and R87_0.2.4 again behave somewhat differently. If we neglect these two runs, $\beta_{{\rm sh}}$ remains unchanged but $\sigma$ is reduced by about 40%. The third panel of Fig. 5 shows that larger upwind Mach-numbers lead to less inclined shocks with respect to the direction of the upstream flows (lower values of $\alpha$). Shown is the number distribution of $\alpha$, averaged over $10 \le \ell ( N ) \le 70$. Individual runs show a slight shift towards higher values of $\alpha$ as $\ell ( N )$ increases. This shift is, however, small compared to the effect of $M_{{\rm u}}$. The fourth panel of Fig. 5 shows the auto-correlation length $\ell _{{\rm corr}}$. It not only depends on $M_{{\rm u}}$ but is also proportional to $\ell _{{\rm cdl}}$. The best fit is found to be $\ell_{{\rm corr}} \approx 0.7 \ell_{{\rm cdl}} M_{{\rm u}}^{-0.6}$. The fifth panel of Fig. 5 shows $\ell _{{\rm corr}}$scaled with this best fit. From these scaling properties of $\ell _{{\rm corr}}$, we take that higher values of $M_{{\rm u}}$ lead to smaller scale wiggling of the shocks with respect to $\ell _{{\rm cdl}}$.

The absolute value of $\ell _{{\rm corr}}$ clearly depends on the choice of the threshold value in our definition, $\Gamma ( y_{{\rm corr}} ) = 0.5$. Figure 6 illustrates the variation of $\Gamma _{{\rm l}}$ as a function of  $y_{{\rm corr}}$ at the example of run R43_0.2.4. The top panel of Fig. 6 shows that the initially present sinusoidal wiggling of the confining shocks does not get lost until about $\ell ( N ) = 15$, which is rather late compared to the other runs. Mode-like signatures again appear around $\ell ( N ) \ga 50$. Our data give, however, no clear answer to how typical and persistent such signatures are. A basic problem is that their wave length soon becomes comparable (within a factor of 2 or so) to the domain size in the y-direction, which may affect the signatures. From the bottom panel of Fig. 6, on the other hand, it can be taken that $\Gamma _{{\rm l}}$ essentially decreases linearly from 1 to about 0.2. The other simulations show a similar behavior. Consequently, the above scaling properties of $\ell _{{\rm corr}}$should also be obtained if smaller threshold values are used, down to about $\Gamma ( y_{{\rm corr}} )= 0.2$.

Figures 4 and 5 also allow some insight into why runs R5_0.2.4 and R87_0.2.4 sometimes fit not so well. The third panel of Fig. 5 shows that our spatial resolution is barely sufficient for run R87_0.2.4, the largest upwind Mach-number we have considered. The number distribution here peaks at around $\alpha \approx 0.1$. In terms of discrete positions this means that the shock position changes by about 15 cells in the x-direction as one moves from y_j to y_j+1. Run R5_0.2.4, on the other hand, may deviate just because of its low Mach-number. The turbulence within its CDL is subsonic, $M_{{\rm rms}} \approx 0.9$; and with $\eta_{{\rm 2}}^{2}M_{{\rm u}}^{2} \approx 1.1$ and ${\rm corr}(\rho,v) \approx -0.4$ (Fig. 4, top panel), it violates two of the basic assumptions made when deriving the self-similar scaling laws in Sect. 3.2.

In summary, as $M_{{\rm u}}$ increases, the bounding shocks become less inclined with respect to the direction of the upstream flows (smaller $\alpha$), the fraction of upstream kinetic energy that survives the passage through the bounding shocks increases, and the bounding shocks themselves are wiggled on progressively smaller scales (smaller $\ell_{{\rm corr} }/ \ell_{{\rm cdl}}$).

   
4.2.3 Energy balance

	Figure 7: Driving efficiency ( top panel) and best fit $\eta _{{\rm 3}} = (1 - f_{{\rm eff}}) M_{{\rm u}}^{0.7}$ ( bottom panel). For details see text.
Open with DEXTER

Energy input into the CDL occurs only at its confining interfaces. Energy dissipation, on the other hand, occurs throughout the CDL volume. Nevertheless, according to the analysis in Sect. 3.2 both $\dot{\cal E}_{{\rm drv}}$ and $\dot{\cal E}_{{\rm diss}}$ should be independent of the CDL extension if dissipation is only due to shocks and if $\ell _{{\rm cdl}}$ is small compared to Y. The average distance between shocks must then increase and/or the average strength of the shocks must decrease as the CDL grows.

To determine $\dot{\cal E}_{{\rm drv}}$ we must compute the driving efficiency $f_{{\rm eff}} = \dot{\cal E}_{{\rm drv}} / {\cal F}_{{\rm e_{{\rm kin}},u}}$. The corresponding integral in Eq. (24) is evaluated numerically, and the resulting driving efficiency is shown in the top panel of Fig. 7. As can be seen, larger Mach-numbers lead to more efficient driving, and a smaller part of the upstream kinetic energy is thermalized already at the confining shocks. The driving efficiency  $f_{{\rm eff}}$ increases by about a factor of four between runs R5_0.2.4 and R87_0.2.4. Also noteworthy is that the absolute value of the driving power $\dot{\cal E}_{{\rm drv}}$ differs by more than 4 orders of magnitude between runs R5_0.2.4 and R87_0.2.4. The best fit for the assumed Mach-number dependence (minimization of $\sigma(\eta_{{\rm 3}})$ in Eq. (25)) yields $\beta _{{\rm 3}}= -0.7$. The corresponding values of $\eta _{{\rm 3}} = (1 - f_{{\rm eff}}) M_{{\rm u}}^{0.7}$ are shown in the bottom panel of Fig. 7. From the figure we take that the second part of our assumption a), the simple Mach-number dependence of  $f_{{\rm eff}}$, seems justified. The figure also shows that  $f_{{\rm eff}}$, and thus the driving power $\dot{\cal E}_{{\rm drv}}$, is not strictly independent of $\ell _{{\rm cdl}}$ but decreases with increasing $\ell ( N )$. Repeating the best fit analysis but allowing for a linear dependence of $\eta_{{\rm 3}}$ on $\ell ( N )$ again leads to $\beta _{{\rm 3}}= -0.7$, while $\eta_{{\rm 3}}$ changes from 3.1 to 3.6 as $\ell ( N )$ goes from 10 to 70. The average value of $\eta_{{\rm 3}}$ is 3.3. Omission of the extreme runs R5_0.2.4 and R87_0.2.4 does not change the result.

We determine the dissipated energy as $\dot{\cal E}_{{\rm diss}} = \dot{\cal E}_{{\rm drv}} - \dot{\cal E}_{{\rm kin}}$(Sect. 3.2.2), where $\dot{\cal E}_{{\rm kin}}$ is the change per time of the kinetic energy within an average column of the CDL, and $\dot{\cal E}_{{\rm kin}}$ is directly from our simulation data. Figure 8 shows the numerically obtained value $\dot{\cal E}_{{\rm diss}}$ (top panel) and the theoretically expected value (Eq. (34)) $\dot{\cal E}_{{\rm diss}}^{{\rm th}}$(middle panel), both in units of ${\cal F}_{{\rm e_{{\rm kin}},u}} = \rho _{{\rm u}} v_{{\rm u}}^{3}$, as well as the ratio of the two (bottom panel). For better display, the theoretical value, which must not depend on $\ell ( N )$, is shown as a (constant) function of $\ell ( N )$. For the constants in Eq. (34) we used the numerically obtained average values, $\eta _{{\rm 3}} = 3.3$, $\beta _{{\rm 3}}= -0.7$, and $\eta _{{\rm 2}} = 0.21$. We used $\eta _{{\rm 3}} = 2.75$ only for R5_0.2.4, in accordance with the bottom panel of Fig. 7. The numerically obtained value was smoothed for better display using a running mean with window size $\Delta \ell ( N ) = \pm 1$. The effect of the smoothing is illustrated in Fig. 9 with the example of run R11_0.2.4.

Figure 8: Numerically obtained ( top panel) and theoretically expected ( middle panel) energy dissipation in units of the upstream kinetic energy flux density ${\cal F}_{{\rm e_{{\rm kin}},u}} = \rho _{{\rm u}} v_{{\rm u}}^{3}$. The constants in Eq. (34) were set to the best fit values, $\eta _{{\rm 3}} = 3.3$, $\beta _{{\rm 3}}= -0.7$, and $\eta _{{\rm 2}} = 0.21$. We used $\eta _{{\rm 3}} = 2.75$ for run R5_0.2.4 (for details see text). The bottom panel shows the ratio of the two quantities. Individual curves denote the same runs as in Fig. 4. For better display, $\dot{\cal E}_{{\rm diss}}$ was smoothed using a running mean with time window $\Delta \ell ( N ) = \pm 1$.

Open with DEXTER

	Figure 9: Effect of smoothing $\dot{\cal E}_{{\rm diss}}$ with a running mean and window $\Delta \ell ( N ) = \pm 1$, illustrated by run R33_0.2.4. Shown is $\dot{\cal E}_{{\rm diss}}$ in units of ${\cal F}_{{\rm e_{{\rm kin}},u}} = \rho _{{\rm u}} v_{{\rm u}}^{3}$, before (dashed, black) and after (solid, red) smoothing, in units of erg cm-3 s-1.
Open with DEXTER

Looking at the data of $\dot{\cal E}_{{\rm diss}}$ and $\dot{\cal E}_{{\rm drv}}$, three points may be stressed. First, $\dot{\cal E}_{{\rm diss}}$ (Fig. 8, top panel) mirrors $\dot{\cal E}_{{\rm drv}} = {\cal F}_{{\rm e_{{\rm kin}},u}} f_{{\rm eff}}$(Fig. 7, top panel), and the values usually differ by less than 10%. This is not surprising. It implies, however, that for larger upstream Mach-numbers, a larger fraction of the upstream kinetic energy is thermalized only within the volume of CDL and not already at its confining shocks. For $M_{{\rm u}} \ga 20$, the energy dissipated within the CDL exceeds 50% of the upstream kinetic energy (Fig. 8, top panel).

Second, the bottom panel of Fig. 8 shows that $\dot{\cal E}_{{\rm diss}}^{{\rm th}}$ and $\dot{\cal E}_{{\rm diss}}$ agree to within 10% most of the time. Given the wide range covered (5 orders of magnitude in $\dot{\cal E}_{{\rm diss}}$, a factor of 20 in  $M_{{\rm u}}$, and an increase by a factor of 7 in  $\ell ( N )$), we conclude that the self-similar solution gives a good estimate.

Third, from the same figure it can be seen that $\dot{\cal E}_{{\rm diss}}$ generally decreases, except for run R5_0.2.4. Excluding R5_0.2.4, a linear fit to  $\dot{\cal E}_{{\rm diss}} / \dot{\cal E}_{{\rm diss}}^{{\rm th}}$ yields a decrease of 10% as $\ell ( N )$ increases from 10 to 70. A similar fit to $\dot{\cal E}_{{\rm drv}} / \dot{\cal E}_{{\rm drv}}^{{\rm th}}$ with $\dot{\cal E}_{{\rm drv}}^{{\rm th}} = \rho_{{\rm u}} v_{{\rm u}}^{3} ( 1 - 3.3 M_{{\rm u}}^{-0.7})$ yields an even slightly larger decrease of 13%. The net dissipation, $\dot{\cal E}_{{\rm diss}} / \dot{\cal E}_{{\rm drv}}$, in fact increases by 3%. Thus, as the CDL size increases, the absolute dissipation within an average column decreases while the net dissipation increases.

In summary, the predicted scaling laws, Eqs. (32) to (34), are - within the range of applicability - essentially confirmed by the simulations. The fraction of upstream kinetic energy dissipated only within the CDL, and not at the confining shocks, thus increases with  $M_{{\rm u}}$. Best-fit analysis for the numerical constants yields  $f_{{\rm eff}} = 1 - 3.3 \; M_{{\rm u}}^{-0.7}$. Both $\dot{\cal E}_{{\rm drv}}$ and $\dot{\cal E}_{{\rm diss}}$ decrease slightly with increasing $\ell _{{\rm cdl}}$. The net dissipation rate $\dot{\cal E}_{{\rm diss}} / \dot{\cal E}_{{\rm drv}}$ increases, but only slightly (3% increase as $\ell ( N )$ goes from 10 to 70).

   
4.2.4 Length scales of the turbulence

In Sect. 4.2.2 we looked at the scaling properties of the confining shocks and pointed out that shorter auto-correlation lengths  $\ell _{{\rm corr}}$ imply smaller-scale wiggling, thus smaller scale changes of the kinetic energy entering the CDL. In the following, we show that the interface based quantity  $\ell _{{\rm corr}}$ is proportional to the length scale derived from the volume properties of the turbulence. We take this as evidence of the tight coupling between volume and interface properties, between the turbulence and its driving.

On dimensional grounds, we can define two length scales based on volume properties of the turbulence,

  

$$\ell_{{\rm e_{{\rm kin}}}} \equiv \frac{ N ^{-1/2} {\cal E}_{{\rm kin}}^{3/2}}{\dot{\cal E}_{{\rm diss}}},$$ (41)

$$\ell_{{\rm v_{{\rm rms}}}} \equiv \frac{ N v_{{\rm rms}}^{3}}{\dot{\cal E}_{{\rm diss}}},$$ (42)

where ${\cal E}_{{\rm kin}} = \frac{\ell_{{\rm cdl}}}{{2V}}\int_{{V}} \rho v^{2}$ is the average column integrated kinetic energy density. Here V is again the 2D volume of the CDL, introduced in Sect. 2.3.1. The two scales are equal up to a numerical constant if the density and velocity are uncorrelated, in which case we can replace the average over the product $\rho v^{2}$ by the product of the averages of $\rho$ and v², ${\cal E}_{{\rm kin}} = \ell_{{\rm cdl}} \rho_{{\rm m}} v_{{\rm rms}}^{2} = N v_{{\rm rms}}^{2}$. As this is the case in most of our simulations we look at only one of the above quantities in the following, $\ell _{{\rm e_{{\rm kin}}}}$, shown in the top panel of Fig. 10. For better display, as $\ell _{{\rm e_{{\rm kin}}}}$ inherits the large time variability of $\dot{\cal E}_{{\rm diss}}$, it is smoothed in the same way as $\dot{\cal E}_{{\rm diss}}$ in the bottom panel of Fig. 8.

Figure 10: Characteristic length $\ell _{{\rm e_{{\rm kin}}}}$ of the turbulence ( top), in units of $\ell _{{\rm corr}}$ ( middle), and scaled with best-fit $\ell _{{\rm cdl}} M_{{\rm u}}^{0.6}$ ( bottom) as functions of $\ell ( N )$. Individual curves denote the same runs as in Fig. 4. For better display, $\ell _{{\rm ekin}}$ was smoothed by a running mean with window $\Delta \ell ( N ) = \pm 1$.

Open with DEXTER

Figure 11: Plots of ${\rm div}(\vec{v})$ for two runs that are identical except for their upstream Mach-number. Larger upstream Mach-numbers lead, on average, to finer structure within the CDL and smaller scale wiggling of the confining shocks. Shown are runs R33_0.2.4 ( top) and R11_0.2.4 ( bottom), both at a time when $\ell_{{\rm cdl}} \approx 2 ~ {Y}_{{\rm 0}} = {Y}/2$. Blue (dark lines) indicates convergence, red (dark patches) divergence.

Open with DEXTER

Assuming a relation of the form $\ell_{{\rm e_{{\rm kin}}}} = \alpha_{{\rm e_{{\rm kin}}}} \ell_{{\rm corr}}$, we obtain optimal fits (minimum of $\sigma^{2}(\alpha_{{\rm e_{{\rm kin}}}})$) for $\alpha_{{\rm e_{{\rm kin}}}} \approx 1.3$. The fits become only slightly better if a weak linear dependence of $\alpha_{{\rm e_{{\rm kin}}}}$ on $\ell ( N )$ is allowed (13% change as $\ell ( N )$ goes from 10 to 70). $\ell_{{\rm e_{{\rm kin}}}} / \ell_{{\rm corr}}$ is shown in the middle panel of Fig. 10. Looking directly at the dependence of $\ell _{{\rm e_{{\rm kin}}}}$ on $\ell _{{\rm cdl}}$ and $M_{{\rm u}}$, we find $\ell_{{\rm e_{{\rm kin}}}} \propto \ell_{{\rm cdl}} M_{{\rm u}}^{-0.6}$. This is the same dependence we found for $\ell _{{\rm corr}}$ in Sect. 4.2.2, $\ell _{{\rm e_{{\rm kin}}}}$ scaled with this best fit is shown in the bottom panel of Fig. 10.

With increasing upstream Mach-number the characteristic length scale $\ell _{{\rm e_{{\rm kin}}}}$ thus decreases with respect to the CDL extension. This is consistent with our observation that for the same  $\ell _{{\rm cdl}}$ the interior of the CDL shows finer structuring (patches, filaments) for higher values of  $M_{{\rm u}}$. Figure 11 illustrates this observation with the example of runs R11_0.2.4 and R33_0.2.4. Shown in the figure is  ${\rm div}(\vec{v})$, as the flow patterns, especially shocks, are better visible in this quantity than in density.

In summary, our simulations show that the inherent length scale of the turbulence is proportional to the auto-correlation length of the confining shocks, independent of  $M_{{\rm u}}$ and  $\ell _{{\rm cdl}}$. With increasing $M_{{\rm u}}$, both length scales decrease relative to the CDL extension, $\ell_{{\rm e_{{\rm kin}}}} / \ell_{{\rm cdl}} \propto M_{{\rm u}}^{-0.6}$. The appearance of the CDL, the size of its patches and filaments, behaves similarly.

   
4.3 Settings with CDL at t = 0

Figure 12: Comparing runs with and without an initial CDL. Shown are $\rho _{{\rm m}}/\rho _{{\rm u}}$ (first), $M_{{\rm rms}}/M_{\rm u}$ (second), the scaled driving efficiency $(1-f_{{\rm eff}}) M_{{\rm u}}^{0.7}$ (third), and the scaled characteristic length of the turbulence $\ell _{{\rm e_{{\rm kin}}}}/\ell _{{\rm cdl}} \cdot M_{{\rm u}}^{0.6}$ for all symmetric runs. Line styles and colors denote initial conditions, 0 (solid line, blue), 1 (dashed line, red), and 2 (dash-dotted line, orange).

Open with DEXTER

Figure 13: Time evolution of angle distribution for run R22_2.2.2. Shown is the average angle distribution for $10 < \ell ( N ) < 70$ (dashed, blue), $70 < \ell ( N ) < 130$ (dash-dotted, green), $130 < \ell ( N ) < 190$ (dash-three-dots, orange), $190 < \ell ( N ) < 250$ (long dashes, purple), $250 < \ell ( N ) < 310$ (solid, red). Also shown are the distributions for run R5_0.2.4 (black dots, right line) and for run R11_0.2.4 (black dots, left line), both averaged over $10 < \ell ( N ) < 70$.

Open with DEXTER

We performed additional runs to study the influence of an initially present CDL. Figure 12 illustrates the results for some selected quantities. Shown are all the runs we performed with initial condition I=0 (no CDL at t=0), I=1 (moderate CDL at t=0), and I=2 (massive CDL at t=0).

Comparison of the I=1 and I=0 curves in Fig. 12 shows that an initially present CDL of moderate column density ( $N = 14 ~ N_{{\rm0}}$) soon develops characteristics similar to those found in simulations without initial CDL. A quasi-steady state is reached for $\ell( N ) \ga 40$. The I=1 and I=0 curves then agree to within about a factor of two for volume quantities like $\rho_{{\rm m}}$ and  $M_{{\rm rms}}$ (first two panels in Fig. 12). Agreement seems slightly better for interface related quantities. For $(1-f_{{\rm eff}})~M^{0.7}$, shown in the third panel of Fig. 12, the I=1 and I=0 curves lie more or less on top of each other. The same is true for $\ell_{{\rm e_{{\rm kin}}}}/ \ell_{{\rm cdl}} ~M_{{\rm u}}^{0.6}$, shown in the bottom panel of Fig. 12.

The situation is slightly different for runs with an initially rather massive CDL (I=2, with initially $N = 56~ N_{{\rm0}}$). Also in these simulations the CDL gets more and more turbulent. For all the quantities shown in Fig. 12, the I=2 curves approach the I=1 and I=0 curves. However, it takes these runs much longer to saturate. Only for $\ell( N ) > 240$ the curves finally seem to saturate, at similar values as the I=0 and I=1curves. That saturation does indeed occur around that time is also supported by Fig. 13. As can be seen, the average angle distribution of the confining shocks for run R22_2.2.2 first shifts to higher and higher values as $\ell ( N )$ increases. It then stagnates for the last two averaging periods, $190 < \ell ( N ) < 250$ and $250 < \ell ( N ) < 310$.

In summary, we conclude that our symmetric simulations all end up in a similar quasi-steady final state. An initially present CDL only delays the development. The incoming flows also manage to generate (and sustain) a similar level of turbulence also within an initially massive CDL.

Figure 14: Average $f_{{\rm eff}}$ as a function of $M_{{\rm rms}}$ for all our symmetric simulations (triangles). In addition, we included data from our asymmetric runs (asterisks), for which $1.6~M_{{\rm r}} \le M_{{\rm l}} \le 64~M_{{\rm r}}$ and which initially have no CDL. Averages were taken over $10 \le \ell ( N ) \le 70$ for simulations without initial CDL (blue triangles and green asterisks), over $40 \le \ell ( N ) \le 70$ for runs with a moderate initial CDL (red triangles), and over $70 \le \ell ( N ) \le 140$ for runs with a massive initial CDL (orange triangles). Lines show $f_{{\rm eff}} = 1 - M_{{\rm rms}}^{\xi }$ with $\xi =-0.6$ (dashed) and $\xi =-0.6 \pm 0.1$ (dotted).

Open with DEXTER

   
4.4 Asymmetric cases

We also computed a few asymmetric cases, where the two upwind Mach-numbers are different, $M_{{\rm l}} \ne M_{{\rm r}}$. For the same reason as given in Sect. 3, we expect the solution to only depend on $M_{{\rm l}}$ and $M_{{\rm r}}$. These dependencies are more complicated than those assumed in Sect. 3 as we now have two different upwind Mach-numbers. The simple dependencies of Sect. 3 should, however, be recovered in the limit $M_{{\rm l}} \rightarrow M_{{\rm r}}$.

The basic physical reason for the more complicated dependencies on the upwind Mach-numbers lies in the strong back coupling between the turbulence within the CDL and the driving of this turbulence by the upwind flows. Our asymmetric simulations demonstrate clearly (much more clearly than the symmetric simulations) that the turbulence crucially affects the driving: although  $M_{{\rm l}}$ and  $M_{{\rm r}}$ are strongly different, the corresponding driving efficiencies are about equal, $f_{{\rm eff,l}} \approx f_{{\rm eff,r}}$. Thus the efficiency does not depend primarily on the upwind flow. In fact, Fig. 14 shows that for both symmetric and asymmetric runs  $f_{{\rm eff}}$ (averaged now over both shocks) can be described well by

$$f_{{\rm eff}} = 1 - M_{{\rm rms}}^{-0.6}.$$ (43)

The angle distribution of the two shocks behaves accordingly in that it is similar for both shocks and determined by $M_{{\rm rms}}$rather than by either $M_{{\rm l}}$ or $M_{{\rm r}}$.

A more detailed analysis of the asymmetric case, including an approximate analytical solution, will be presented in a subsequent paper.

   
4.5 Grid and domain studies

The numerical results presented in Sect. 4.2 were all based on simulations with a domain ${Y}=4 {Y}_{{\rm0}}$ and a discretization of $1.5 \times 10^{-3} {Y}_{{\rm0}}$ (R=2) or 2560 cells in the y-direction. Here we want to check whether these choices have any systematical effect on the numerical results of Sect. 4.2.

   
4.5.1 Different y-extent

To check whether the size of the computational domain has any systematic effect on the results of Sect. 4.2, we performed some of the simulations again, but this time on smaller domains of ${Y}=2 {Y}_{{\rm0}}$ and ${Y}={Y}_{{\rm0}}$. We also performed one simulation on a larger domain ${Y}=6 {Y}_{{\rm0}}$.

Figure 15 illustrates our findings for simulations on domains ${Y}=2 {Y}_{{\rm0}}$ and  ${Y}=4 {Y}_{{\rm0}}$. $M_{{\rm rms}}$ shows no systematic effect and is, as such, representative of other volume-related quantities (Fig. 15, top panel). As a typical representative for interface-related quantities, $f_{{\rm eff}}$ also shows no clear overall effect of the domain size (Fig. 15, middle panel). The quantity for which we find the most clear effect is the auto-correlation length  $\ell _{{\rm corr}}$ (Fig. 15, bottom panel). However, even for  $\ell _{{\rm corr}}$ the effect sets in only for two of the four runs and only for $\ell( N ) \ga 30$, i.e. once the CDL extension reaches about half the size of the smaller domain. For the numerical results in Sect. 4.2, $\ell _{{\rm cdl}} \approx {Y}/2$ corresponds to $\ell ( N ) = 60$. We conclude that the y-extent of the computational domain has no apparent systematic effect on these results up to  $\ell( N ) \la 30$ and probably even up to  $\ell( N ) \la 60$.

A systematic effect of the computational domain on the numerical solution does become apparent if the simulations are carried on much longer. One pair of runs, R22_0.2.2 and R22_0.2.6, were carried on much longer, till $\ell ( N ) \approx 200$. For this pair of runs, Fig. 16 shows the evolution of $M_{{\rm rms}}$ for each run, as well as their ratio, $M_{{\rm rms,3y}}/M_{{\rm rms,y}}$. The run on the smaller domain apparently shows a faster decay in $M_{{\rm rms}}$ after $\ell ( N ) \approx 100$. From Fig. 17 we take that the behavior of this one pair of runs is most likely the rule, and not the exception. The top panel of Fig. 17 shows $\ell _{{\rm corr}}$, scaled, for all the symmetric runs we have performed and whose domain has a y-extent $\le$ $2 {Y}_{{\rm0}}$. The bottom panel of Fig. 17 gives the same quantity for all the runs with a domain extention $\ge$ $4{Y}_{{\rm0}}$. Comparison of the two figures shows that runs on a domain $\le$ $2 {Y}_{{\rm0}}$ saturate around $\ell_{{\rm corr}} M_{{\rm u}}^{0.6} \approx 1.6 {Y}_{{\rm0}}$. For runs on a domain $\ge$ $4{Y}_{{\rm0}}$, $\ell _{{\rm corr}}$ reaches much higher values.

Figure 15: Comparing runs that differ only in the y-extent of the domain ( ${Y} = 2~{Y}_{{\rm 0}}$ and ${Y} = 4~{Y}_{{\rm 0}}$). Shown are $M_{{\rm rms,2}Y} / M_{{\rm rms,}Y}$ ( top), $f_{{\rm eff,2}Y} / f_{{\rm eff,}Y}$ ( middle), and $\ell _{{\rm corr,2}Y} / \ell _{{\rm corr,}Y}$ ( bottom). Individual curves denote runs R11_0.2.* (dashed, purple), R22_0.2.* (solid, red), R33_0.2.* (dash-dotted, orange), R43_0.2.* (dash-three-dots, green).

Open with DEXTER

	Figure 16: Comparison of runs R22_0.2.2 and R22_0.2.6, illustrating the effect of a three-times larger y-extent of the computational domain on long time scales. Shown are $M_{{\rm rms}}/M_{\rm u}$ ( top) for R22_0.2.2 (solid, light blue) and R22_0.2.6 (dashed, dark red) and the ratio $M_{{\rm rms,3}Y} / M_{{\rm rms,}Y}$ ( bottom).
Open with DEXTER

	Figure 17: Scaled auto-correlation lengths of all symmetric simulations on domains with a y-extent less or equal to $2~{Y}_{{\rm 0}}$ ( top) and a y-extent greater or equal to $4 ~ {Y}_{{\rm0}}$ ( bottom).
Open with DEXTER

   
4.5.2 Different discretization

The results presented in Sect. 4.2 were all based on simulations with a discretization of $1.5 \times 10^{-3} {Y}_{{\rm0}}$ (R=2) or 2560 cells in the y-direction. To check the effect of the discretization on our results, we repeated several simulations with coarser and/or finer discretization. These simulations indeed reveal a systematic effect of the discretization on the values of average quantities. Nevertheless, the general properties of the solution, its approximate self-similarity and Mach-number dependences, remain unaltered. Only the numerical constants $\eta_{{i}}$ are affected. The changes are, however, small when compared to the differences between the 1D and 2D solution (for example, $\rho_{{\rm m}} = \eta_{{\rm 1}} \rho_{{\rm u}}$ in 2D, while $\rho_{{\rm m}} = M_{{\rm u}}^{2} \rho_{{\rm u}}$ in 1D).

We find that finer discretization generally leads to reduced turbulence. Using finer meshes we obtained larger mean densities and lower values of  $M_{{\rm rms}}$, as shown in Fig. 18. The driving efficiency gets lower and the shocks become more inclined with respect to the direction of the upstream flows, and the angle distribution is shifted to lower values. The characteristic length scale $\ell _{{\rm e_{{\rm kin}}}}$ remains about constant if taken in units of $\ell _{{\rm cdl}}$.

A possible explanation for the reduction of turbulence (smaller  $M_{{\rm rms}}$) on finer grids could be the dominance of shocks for the energy dissipation in the CDL. On a coarser grid, the network of shocks within the CDL is less dense. The divergence plots shown in Fig. 19 illustrate this effect. A closer analysis of this idea is, however, beyond the scope of the present paper.

	Figure 18: Comparison of $M_{{\rm rms}}$ for runs whose spatial resolution differs by a factor of 2 (subscript c = coarse, f = fine). Shown are (giving only the name of the finer run) runs R22_0.2.4 (solid, red), R22_0.4.4 (dash-three-dots, blue), R43_0.2.4 (long dashes, purple), and R11_0.2.4 (dash-dotted, orange).
Open with DEXTER

Figure 19: Plots of ${\rm div}(\vec{v})$ for two runs that are identical to run R11_0.2.4, shown in Fig. 11, except for their discretization. The runs shown here were computed with two times lower ( top) and at four times lower ( bottom) resolution. Blue (dark lines) indicates convergence, red (dark patches) divergence. As can be seen, the number of convergent regions within an average CDL column decreases with decreasing resolution.

Open with DEXTER

We stress that, so far, no convergence has been reached in our discretization studies. Looking at the comparison of the three runs R22_0.1.4, R22_0.2.4, and R22_0.4.4 in Fig.18 shows that each reduction of the cell size by a factor of two leads to a reduction of about 20% in $M_{{\rm rms}}$. This indicates that the resolution of 2560 cells in y-direction in our standard runs (R*_0.2.4) and of 5120 cells in the y-direction in the refined runs is still not sufficient. This should be kept in mind when interpreting these results or any results on shock bound turbulent structures in 2D, let alone 3D.

Also, no clear picture emerges regarding the deviation of $M_{{\rm rms}}$from the constant value predicted by Eq. (15). A linear fit to $M_{{\rm rms}}$ for  $10 \le \ell ( N ) \le 70$ yields -12% for run R22_0.2.4 and -23% for the two times coarser run R22_0.1.4. For runs R43_0.*.4, the grid dependence is the other way round: R43_0.2.4 shows a decrease of -25%, the twice coarser run R43_0.1.4 decreases by only -15%.

   
5 Discussion

We want to address four points in this section. First, we sketch possible reasons for the slight difference between the numerical solution and the relations we derived in Sect. 3. Second, we look once more at the driving of the turbulence and, in particular, the back-coupling between interface and volume properties. Third, we briefly consider our results in an astrophysical context, in particular with regard to molecular clouds. Finally, based on preliminary numerical results, we sketch the effect of some additional physics.

5.1 The numerical solution versus the analytical solution

In Sect. 3.2 we suggested that a self-similar solution to our 2D model problem may still exist for the limiting case where the system approaches infinity. The relations derived in that section give a reasonable estimate for the numerical results of Sect. 4. However, while $M_{{\rm rms}}$ is constant in Sect. 3.2, the numerical simulations show a gradual decrease in $M_{{\rm rms}}$ already for small CDLs, $\ell_{{\rm cdl}} \la {Y}/2$ (15% decrease of  $M_{{\rm rms}}$ as $\ell ( N )$ increases from 10 to 70, Sect. 4.2). We have no firm explanation for this difference. We sketch three possible effects in the following, but stress that the available data do not allow us to clearly distinguish between them.

A first, obvious reason could be the finite y-extent of the computational domain, Y. It sets an upper limit on the total energy input into the CDL, thus on the amount of mass within the CDL that can be driven. Once the CDL has accumulated too much mass, the driving per unit mass weakens and the turbulence starts to weaken. The spatial growth of the CDL slows down while the average density increases. The following considerations on time scales may illustrate this point further.

An upper limit to the time at which Y starts to affect the solution is given by the time t_y at which $\ell_{{\rm cdl}} = {Y}$. At later times structures may still grow in the x-direction (up to $\ell _{{\rm cdl}}$ at most) but cannot grow any more in the y-direction (where Y sets an upper limit). For the runs in Sect. 4.2, $\ell_{{\rm cdl}} = {Y}$corresponds to $\ell ( N ) \approx 120$ or $t_{{y}} = 12 {Y}_{{\rm0}}/v_{{\rm rms}}$. A lower limit for the decay time scale of the turbulence may be obtained as follows. For the case of uniformly driven isothermal hydrodynamic turbulence in a 3D periodic box, Mac Low (1999) has shown that the typical decay time once the driving is turned off, $t_{{\rm0}}$, and the initial driving wave length, $\lambda_{{\rm drv}}$, are related by $t_{{\rm0}} \approx \lambda_{{\rm drv}} / v_{{\rm rms}}$. Assuming that this result also holds for our slab, that $\lambda_{{\rm drv}} = {Y}$, and that driving is turned off completely, it follows that $t_{{\rm0}} \approx {Y} / v_{{\rm rms}}$, or $t_{{\rm0}} \approx 4 {Y_{{\rm0}}} / v_{{\rm rms}}$ for the runs in Sect. 4.2. However, driving continues in our simulations and so the effective decay time scale of the turbulence is likely to be much longer than $t_{{\rm0}}$. Finally, for the runs in Sect. 4.2, and a typical integration time of $\ell ( N ) = 60$ corresponds to about $\tau = 6 {Y}_{{\rm0}} / v_{{\rm rms}}$, a typical turbulent crossing time at $\ell ( N ) = 60$) is $\tau_{{\rm cross}} = \ell_{{\rm cdl}} / v_{{\rm rms}} \approx 2 {Y}_{{\rm0}}/ v_{{\rm rms}}$. Comparing these different time scales makes it seem likely that at $\ell ( N ) = 60$, turbulence in the center of the CDL is still essentially driven, not essentially decaying.

Our simulation data do not allow us to either clearly confirm or reject the hypothesis that the finite y-extent of the computational domain is responsible for the slight decrease in $M_{{\rm rms}}$ that we observe at early times, $\ell ( N ) \la 70$. If the finite domain size were responsible, $M_{{\rm rms}}$ should decay differently on different domains. Comparison of simulations on different domains up to $\ell ( N ) \approx 70$(Sect. 4.5.1), however, gives no clear picture. The data are rather noisy, and simulations on domains $2 {Y}_{{\rm0}}$ and $4{Y}_{{\rm0}}$ show no systematic differences as long as $\ell( N ) \la 30$ ( $\ell_{{\rm cdl}} < {Y}/2$ on the smaller domain). Only for much later times, $\ell ( N ) \gg 70$, well beyond the range for the results in Sect. 4.2, does Y have a clear effect and $M_{{\rm rms}}$ decreases faster on smaller domains (Fig. 16).

A second, more speculative, reason might be numerical dissipation, provided that its effect were to increase with $\ell _{{\rm cdl}}$. While we have no evidence that the latter is really the case, it may also be hasty to discard this possibility right away. Porter & Woodward (1994) found, by observing how simple 2D hydrodynamical flows (shear flows and sound waves of definite wave number, their Sect. 3.3) damp with time, that the decay rate due to numerical dissipation alone is a non-linear function of the wave number. Their results are certainly not directly applicable to the present case. But in view of these results, and given the change in structure size with $\ell _{{\rm cdl}}$ as suggested by Fig. 3, it might be possible that the effect of numerical dissipation indeed changes with $\ell _{{\rm cdl}}$. Note that this would also imply that the MILES approach, outlined in Sect. 2.1, were not strictly valid for the problem we consider. The currently available data do not allow us to clearly reject or confirm the effect.

A third reason, or rather an amplifying mechanism, could be back-coupling between $M_{{\rm rms}}$ and the driving efficiency. Once the turbulence within the CDL is slightly reduced (for whatever reason), the reduction is further amplified by the back-coupling between turbulence and driving, $f_{{\rm eff}} = (1 - M_{{\rm rms}}^{-0.6})$. The decrease in $M_{{\rm rms}}$ results in larger inclination of the shocks with respect to the direction of the upstream flows, more energy is dissipated at the confining shocks of the CDL, and less driving energy enters the CDL. For the observed 15% reduction of $M_{{\rm rms}}$, the reduced driving may, in fact, play a dominant role: as $\ell ( N )$ increases from 10 to 70, $\dot{\cal E}_{{\rm drv}} / \dot{\cal E}_{{\rm drv}}^{{\rm th}}$decreases by 13% (Sect. 4.2.3). But to really estimate the relative importance of the three effects just sketched, further studies are certainly necessary.

Two more points seem noteworthy to us in this section. One concerns the near independence of $\dot{\cal E}_{{\rm diss}}$ on $\ell _{{\rm cdl}}$. From Fig. 3 (increase in structure size with increasing $\ell _{{\rm cdl}}$), we take that it is rather the increasing average distance between shocks that allows $\dot{\cal E}_{{\rm diss}}$ to be essentially independent of $\ell _{{\rm cdl}}$ and not so much the, on average, decreasing strength of shocks (Sect. 3.2.3). Whether this is indeed true, only a closer analysis of the structure within the CDL along the lines of Mac Low & Ossenkopf (2000) can tell, which is, however, beyond the scope of the present paper. Such an analysis could also shed light on whether (or in which sense) $\ell _{{\rm e_{{\rm kin}}}}$(see Sect. 4.2.3) is indeed a measure of the average distance between shocks. It would also allow us to quantify our impression that small scale structures are preferably located close to the confining interfaces. If true, this would fit with the result by Smith et al. (2000) that the high-frequency part of the shock spectrum is lost most efficiently.

The other point concerns run R5_0.2.4. With corr $(\rho,v) \approx -0.4$ $M_{{\rm rms}} \approx 0.9$, it violates two of the basic assumptions we made in Sect. 3.2. Its mean density is close to the isothermal value for strong shocks, $\rho_{{\rm m}} \approx 22 \rho_{{\rm u}} \approx 0.9 \rho_{{\rm u}} M_{{\rm u}}^{2}$. Both $\dot{\cal E}_{{\rm diss}}$ and $\dot{\cal E}_{{\rm drv}}$ increase with $\ell _{{\rm cdl}}$. With these characteristics, R5_0.2.4 may mark the transition from compressible supersonic turbulence, the topic of this paper, to compressible subsonic turbulence.

5.2 CDL and confining shocks: a coupled system

The turbulence within the CDL is "naturally driven'' in the sense that we control neither what fraction of the total upstream kinetic energy, $\rho_{{\rm u}} M_{{\rm u}}^{2}$, really enters the CDL nor the spatial scale on which this energy input varies. Both are directly determined by the confining shocks instead and indirectly depend on the system as a whole. The driving efficiency at each confining shock scales with $M_{{\rm rms}}$, even for situations where $M_{{\rm l}} \ne M_{{\rm r}}$ (see Sect. 4.4). The auto-correlation length of the confining shocks and the characteristic length scale of the turbulence within the CDL are proportional to each other, both scaling as $\ell_{{\rm cdl}} M_{{\rm u}}^{-0.6}$. We take these facts as evidence that the CDL as a whole, its interface and volume properties, forms a tightly coupled, quasi-stationary, and self-regulating system. Back coupling between post shock flow and shock is also described in other contexts, for example by Foglizzo (2002) for the case of Bondi-Hoyle accretion.

An aspect that remained elusive in Sect. 4 is the spatial scale on which the energy input varies, the energy injection scale. To really tackle this issue, it would be necessary to analyze the energy spectrum of the CDL. This task requires, however, some caution because of the highly irregular boundary of the CDL, and we postpone it for the moment. Nevertheless, we would like to present a few thoughts on the subject.

A first question is whether it is justified to speak at all of only one injection scale, of monochromatic driving. The homogeneous upstream flow is modulated by the confining shocks. These are wiggled on a variety of spatial scales at any given moment. This strongly suggests that the kinetic energy input into the CDL is most likely not monochromatic but occurs at a whole spectral range instead. Consequences of such non-monochromatic driving have been studied, for example, by Norman & Ferrara (1996).

It also seems worthwhile to briefly look at monochromatically-driven turbulence, in particular at the numerical simulations by Mac Low (1999). For the case of artificially, monochromatically driven hydrodynamic turbulence in a 3D box with periodic boundaries, he found that the characteristic length of the turbulence is proportional to the driving wave length, independent of the Mach-number: $\lambda / \ell^{{\rm 3d}}_{{\rm e_{{\rm kin}}}} = 1.42$, where $\lambda$ is the (known) driving wave length and $\ell^{{\rm 3d}}_{{\rm e_{{\rm kin}}}}$ is the 3D analogon of $\ell _{{\rm e_{{\rm kin}}}}$ in Eq. (41). In addition, Mac Low (1999) observed that $\ell^{{\rm 3d}}_{{\rm e_{{\rm kin}}}}$ increases with $\lambda$, which is mirrored in the apparent increase in the structure size (patches, filaments).

Although our setting clearly differs from that of Mac Low (1999), two thoughts come to mind. The first is an actual observation, namely that we also observe an increase in structure size with $\ell _{{\rm e_{{\rm kin}}}}$. The second thought is more of a question or speculation. Mac Low (1999) determines the proportionality constant between the characteristic scale of the turbulence and the monochromatic driving wave length. One may wonder about the implications of this finding if not one driving wave length is present but a whole spectrum. How will the characteristic length scale of the turbulence, which can still be determined following Eq. (41), depend on this spectrum? And, given our finding that $\ell_{{\rm e_{{\rm kin}}}} \propto \ell_{{\rm corr}}$, what does $\ell _{{\rm corr}}$ tell us about this spectrum? Both questions should become tractable once the energy spectrum of the CDL is determined.

5.3 A glimpse at astrophysics

With regard to astrophysics, the presented work basically suggests that, within the frame of isothermal hydrodynamics and a roughly plane parallel setting, larger Mach-numbers of the colliding flows results in a finer and finer network of higher and higher density contrast within the interaction zone. In different types of wind-driven structures, this connection between Mach-number and structure may be directly observable.

For the clumping of line-driven hot-star winds, our results suggest that the sheets or clumps formed by the instability of the line-driving are not homogeneous but possess fine-scale substructure with high density contrast.

Concerning molecular clouds, we first mention that recent arguments support the idea, originally brought forward by Hunter (1979) and Larson (1981), that molecular clouds result from the collision of large-scale flows in the ISM. Basu & Murali (2001) make the point that small-scale driving ($\approx$0.1-1 pc) of molecular clouds is incompatible with observed total luminosities, unless the energy dissipation rates derived from MHD simulations are seriously overestimated. Using a principal component analysis of $^{{\rm 12}}$CO (J=1-0) emission, Brunt (2003) identifies large-scale flows of atomic material in which the globally turbulent molecular clouds are embedded. Similar observational results were reported by Ballesteros-Paredes et al. (1999a).

Driven supersonic turbulence as a structuring agent for the interior of molecular clouds was examined by many authors (Audit & Hennebelle 2005; Elmegreen 1993; Ballesteros-Paredes et al. 1999a; Mac Low 1999; Mac Low & Klessen 2004; Hartmann et al. 2001; Vazquez-Semadent et al. 1995; Burkert & Hartmann 2004; Hunter et al. 1986; Vázquez-Semadeni et al. 2006; Joung & Mac Low 2004; Heitsch et al. 2005; Ballesteros-Paredes et al. 2006; Kim & Ryu 2005; Mac Low et al. 1998; Ballesteros-Paredes et al. 1999b). The driving wave length of the turbulence, and thus the largest structure size (Mac Low 1999; Ballesteros-Paredes & Mac Low 2002), is usually a free parameter. Our results show instead that, at least for the case of an isothermal, shock compressed, supersonically turbulent 2D slab, the structure size rather depends on the size of the slab or cloud.

5.4 Additional physics: an outlook

The model presented in this paper covers only some very basic physics. To obtain results with a more direct relation to reality, additional physics must be included in the future, among these the following. Strongly asymmetric flows, where  $M_{{\rm l}} \ne M_{{\rm r}}$, lead to more complicated dependences, as we will demonstrate in a forthcoming paper. Inclusion of radiative cooling, instead of assuming isothermal conditions, can affect the problem in different ways. Thermal instability can lead to additional dynamical effects (Audit & Hennebelle 2005; Chevalier & Imamura 1982; Walder & Folini 1996; Hennebelle & Pérault 2000,1999; Mignone 2005; Koyama & Inutsuka 2002; Heitsch et al. 2005; Pittard et al. 2005; Vázquez-Semadeni et al. 2000; Strickland & Blondin 1995; Gaetz et al. 1988). Extended cooling layers, on the other hand, tend to act as a cushion. Simulations by Walder & Folini (1999) and Walder & Folini (2000b), which include radiative cooling but have otherwise similar parameters as some of the simulations presented here, show comparatively more small scale structure and even roll-ups at the interfaces confining the CDL. The CDL as a whole evolves less violently, and mean densities are about a factor of four to eight higher that what we found here for the isothermal case. Strongly asymmetric flows, where $M_{{\rm l}} \ne M_{{\rm r}}$, lead to a qualitatively different solution if radiative cooling is included (Walder & Folini 1998) and to more complicated dependences on the upwind Mach-numbers in the isothermal case, as we will demonstrate in a forthcoming paper. The role of thermal conduction has only been considered by relatively few publications so far (Koyama & Inutsuka 2004; Begelman & McKee 1990; Myasnikov & Zhekov 1998). Global bending of the interaction zone affects the stability properties of the interaction zone as a whole and thus probably also its interior properties. In colliding wind binaries, for example, matter is transported out of the central part of the system and diluted in the outer part. Simulations of bow shocks and colliding winds in binaries show strong traveling waves, together with a systematic change of the mean properties in the flow off from the stagnation point (Stevens et al. 1992; Walder & Folini 1995; Blondin & Koerwer 1998).

   
6 Summary and conclusions

We looked at symmetric, supersonic ( $5 \la M_{{\rm u}} \la 90$), isothermal, plane-parallel, colliding flows in 2D. The resulting shock-confined interaction zone (CDL) is supersonically turbulent ( $1 \la M_{{\rm rms}} \la 10$). We investigated the CDL and its interplay with the upstream flows by dimensional analysis and numerical simulations. The latter we generally stopped when $\ell _{{\rm cdl}} \approx {Y}/2$. The results are interesting not only with regard to flow collisions, but also shed new light on the properties of supersonic turbulence in general.

The numerical simulations show that the CDL has an irregular shape and a patchy, supersonically turbulent interior. The driving of the turbulence is natural in that it depends on the shape of the confining shocks. The dimensional analysis is based on isothermal Euler equations in infinite space. Within this frame, a self-similar solution may exist that would depend on  $M_{{\rm u}}$ but must not depend on $\ell _{{\rm cdl}}$. Relations for average quantities are obtained under some further simplifying assumptions (Sect. 3.3).

Based on both the analytical and numerical results, we arrive at the following conclusions.

1.

Comparison of the numerical and the self-similar solution shows generally good agreement if $M_{{\rm rms}} \ga 1$. The modest deviation between the numerical and the self-similar solutions increases with $\ell _{{\rm cdl}}$. We suggest some explanations for the deviation, but our data do not allow any clear conclusions on the issue. For $M_{{\rm rms}} \la 1$, we have but one simulation. It shows clear differences to the other runs and may be more characteristic of compressible subsonic turbulence than of supersonic turbulence;

2.

The CDL is characterized by $M_{{\rm rms}} \approx \eta_{{\rm 1}}^{-1/2} M_{{\rm u}}$ and $\rho_{{\rm m}} \approx \eta_{{\rm 1}} \rho_{{\rm u}}$. The average compression ratio of the CDL is thus independent of $M_{{\rm u}}$. This is in sharp contrast to the 1D case, where $\rho_{{\rm m,1d}} = M_{{\rm u}}^{2}\rho_{{\rm u}}$. From the numerical simulations, we find $\eta_{{\rm 1}} \approx 30$;

3.

The turbulence within the CDL and the driving efficiency are related by $f_{{\rm eff}}= 1 - M_{{\rm rms}}^{-0.6}$. The relation also holds for asymmetric settings, where $M_{{\rm l}} \ne M_{{\rm r}}$, emphasizing the mutual coupling between volume and interface properties. For larger upstream Mach-numbers, the shocks confining the interaction zone are less strongly inclined with respect to the direction of the upstream flows. The driving is more efficient, a larger fraction of the upstream kinetic energy is dissipated only within the CDL and not already at the confining shocks;

4.

The characteristic length scale of the turbulence, $\ell _{{\rm e_{{\rm kin}}}}$, and the auto-correlation length of the confining shocks, $\ell _{{\rm corr}}$, are proportional to each other. Both scale as $\ell_{{\rm cdl}} M_{{\rm u}}^{-0.6}$, this although the former is based on volume quantities while the latter is derived from interface properties;

5.

The separation of filaments and the size of patches within the CDL both get larger as $\ell _{{\rm cdl}}$ increases and/or $M_{{\rm u}}$ decreases.

For increasing upstream Mach-numbers in summary we thus expect a faster expanding CDL with less strongly inclined confining interfaces with respect to the direction of the upstream flows, similar mean density, finer interior structure relative to the CDL size, and a gradual shift of the energy dissipation from the confining shocks to internal shocks within the CDL. We expect to observe these general dependencies in real objects where shock-confined slabs play a role, like molecular clouds, wind driven structures, supernova remnants, or $\gamma$-ray bursts.

Acknowledgements

The authors wish to thank the crew running the Cray SV1 at ETH Zürich, where the simulations were performed, the system administrator of the institute for astronomy, ETH Zürich, P. Steiner, for steady support, and J. Favre from the Swiss Center of Scientific Computing CSCS/SCSC, Manno, for graphics support. The authors also would like to thank the referee, E. Vazquez-Semadeni, for the detailed and engaged report.

   
Appendix A: numerical computation of obliqueness angle

While shocks are smeared over approximately 3 grid cells in our simulations, the confining shocks in our analysis are specified as a series of discrete x,y-coordinate pairs only (see Sect. 4.2.2). This information is sufficient to compute most shock-related quantities to good accuracy, for example the shock length  $\ell_{{\rm sh}}$. The only quantity that requires a more careful proceeding is the obliqueness angle $\alpha$. If it were computed directly from the discrete shock positions, only discrete values would be obtained, for example 0$^{\circ}$, 45$^{\circ}$, 63.4$^{\circ}$ etc. for one-sided differences.

To compute the obliqueness angle $\alpha_{{i}}(y_{{j}})$ (see Fig. 1 and Sect. 4.2.2) at each position $y_{{\rm j}}$, $1\le j \le J$, of the left and right shock ( $s_{{\rm l}}$ and  $s_{{\rm r}}$), we proceed as follows. In a first step, we use spline interpolation to double the number of points in the y-direction along the shock front. Next, we smooth the shock front slightly, using a running mean with an averaging window of $\pm$5 points (this corresponds to an averaging window of $\pm$2.5 points in the original data. Then we compute the derivative at each point of this smoothed shock front, using a 3-point Lagrangian interpolation. To avoid abrupt changes in the derivative from one point to the next, we smooth it again by a running mean with averaging window $\pm$5 points. We finally obtain the obliqueness angle $\alpha_{{i}}(y_{{j}})$, $1 \le j \le 2J$, as the arctan of the derivative.

We checked that the size of the averaging window ($\pm$3 points or $\pm$7 points) has only a marginal effect on the angle distribution and the driving efficiency. For the latter, which is an integral over both shocks, tests show that $\alpha$ can even be computed directly from the discrete positions.

   
Appendix B: list of runs, their parameters, and naming schemes

Table B.1: List of performed simulations.

References

• Anninos, P., & Norman, M. L. 1996, ApJ, 460, 556 [NASA ADS] [CrossRef]
• Audit, E., & Hennebelle, P. 2005, A&A, 433, 1 [NASA ADS] [CrossRef] [EDP Sciences]
• Ballesteros-Paredes, J., Hartmann, L., & Vazquez-Semadeni, E. 1999a, ApJ, 527, 285 [NASA ADS] [CrossRef]
• Ballesteros-Paredes, J., Vázquez-Semadeni, E., & Scalo, J. 1999b, ApJ, 515, 286 [NASA ADS] [CrossRef]
• Ballesteros-Paredes, J., & Mac Low, M. 2002, ApJ, 570, 734 [NASA ADS] [CrossRef]
• Ballesteros-Paredes, J., Gazol, A., Kim, J., et al. 2006, ApJ, 637, 384 [NASA ADS] [CrossRef]
• Basu, S., & Murali, C. 2001, ApJ, 551, 743 [NASA ADS] [CrossRef] (In the text)
• Begelman, M. C., & McKee, C. F. 1990, ApJ, 358, 375 [NASA ADS] [CrossRef]
• Berger, M. J. 1985, Lect. Appl. Math., 22, 31 (In the text)
• Blondin, J. M., & Marks, B. S. 1996, New Astron., 1, 235 [NASA ADS] (In the text)
• Blondin, J. M., & Koerwer, J. F. 1998, New Astron., 3, 571 [NASA ADS]
• Boris, J., Grinstein, F., Oran, E., & Kolbe, R. 1992, Fluid. Dynam. Res, 10, 199 (In the text)
• Brunt, C. M. 2003, ApJ, 280 (In the text)
• Burkert, A., & Hartmann, L. 2004, ApJ, 616, 288 [NASA ADS] [CrossRef] (In the text)
• Chevalier, R. A., & Imamura, J. N. 1982, ApJ, 261, 543 [NASA ADS] [CrossRef]
• Colella, P. 1990, J. Comput. Phys., 87, 171 [NASA ADS] [CrossRef] (In the text)
• Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef]
• Corcoran, M. F., Pittard, J. M., Stevens, I. R., Henley, D. B., & Pollock, A. M. T. 2005, in X-Ray and Radio Connections, ed. L. O. Sjouwerman, & K. K. Dyer, Published electronically by NRAO, http://www.aoc.nrao.edu/events/xraydio Held 3-6 February 2004 in Santa Fe, New Mexico, USA (E2.02) 13 pages
• Elmegreen, B. G. 1993, ApJ, 419, L29 [NASA ADS] [CrossRef]
• Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211 [NASA ADS] (In the text)
• Fan, Y. Z., & Wei, D. M. 2004, ApJ, 615, L69 [NASA ADS] [CrossRef]
• Feldmeier, A., & Owocki, S. 1998, Ap&SS, 260, 113 [NASA ADS]
• Feldmeier, A., Puls, J., & Pauldrach, A. W. A. 1997, A&A, 322, 878 [NASA ADS]
• Foglizzo, T. 2002, A&A, 392, 353 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Folini, D., & Walder, R. 2000, in Thermal and Ionization Aspects of Flows from Hot Stars, ASP Conf. Ser., 204, 267
• Gaetz, T. J., Edgar, R. J., & Chevalier, R. A. 1988, ApJ, 329, 927 [NASA ADS] [CrossRef]
• Gammie, C. F., & Ostriker, E. C. 1996, ApJ, 466, 814 [NASA ADS] [CrossRef] (In the text)
• Garnier, E., Mossi, M., Sagaut, P., Comte, P., & Deville, M. 1999, J. Comput. Phys., 153, 273 [NASA ADS] [CrossRef]
• Ghisellini, G., Celotti, A., & Costamante, L. 2002, A&A, 386, 833 [NASA ADS] [CrossRef] [EDP Sciences]
• Hartmann, L., Ballesteros-Paredes, J., & Bergin, E. A. 2001, ApJ, 562, 852 [NASA ADS] [CrossRef]
• Heitsch, F., Burkert, A., Hartmann, L. W., Slyz, A. D., & Devriendt, J. E. G. 2005, ApJ, 633, L113 [NASA ADS] [CrossRef]
• Hennebelle, P., & Pérault, M. 1999, A&A, 351, 309 [NASA ADS]
• Hennebelle, P., & Pérault, M. 2000, A&A, 359, 1124 [NASA ADS]
• Heyer, M. H., & Brunt, C. M. 2004, ApJ, 615, L45 [NASA ADS] [CrossRef]
• Hirsch, C. 1995, Numerical computation of internal and external flows, Vol. 1, Fundamentals of discretization (John Wiley and Sons)
• Hueckstaedt, R. M. 2003, New Astron., 8, 295 [NASA ADS]
• Hunter, J. H. 1979, ApJ, 233, 946 [NASA ADS] [CrossRef] (In the text)
• Hunter, J. H., Sandford, M. T., Whitaker, R. W., & Klein, R. I. 1986, ApJ, 305, 309 [NASA ADS] [CrossRef]
• Jasak, H., & Weller, H. 1995, internal Report, CFD research group, Imperial College, London
• Joung, M. K. R., & Mac Low, M.-M. 2004, Am. Astron. Soc. Meet. A., 204
• Kaiser, C. R., Sunyaev, R., & Spruit, H. C. 2000, A&A, 356, 975 [NASA ADS] (In the text)
• Kang, H., Ryu, D., Cen, R., & Song, D. 2005, ApJ, 620, 21 [NASA ADS] [CrossRef]
• Kim, J., & Ryu, D. 2005, ApJ, 630, L45 [NASA ADS] [CrossRef]
• Klein, R. I., Woods, D. T., & Tod, D. 1998, ApJ, 497, 777 [NASA ADS] [CrossRef] (In the text)
• Koyama, H., & Inutsuka, S.-i. 2002, ApJ, 564, L97 [NASA ADS] [CrossRef]
• Koyama, H., & Inutsuka, S.-i. 2004, ApJ, 602, L25 [NASA ADS] [CrossRef]
• Larson, R. B. 1981, MNRAS, 194, 809 [NASA ADS] (In the text)
• Lele, S. 1992, J. Comput. Phys., 103, 16 [CrossRef] (In the text)
• Lesieur, M. 1999, Turbulence in fluids (fluid mechanics and its applications) (Springer) (In the text)
• LeVeque, R. 2002, Finite volume methods for hyperbolic problems (Cambridge University Press)
• Mac Low, M.-M. 1999, ApJ, 524, 169 [NASA ADS] [CrossRef] (In the text)
• Mac Low, M., & Ossenkopf, V. 2000, A&A, 353, 339 [NASA ADS] (In the text)
• Mac Low, M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125 [NASA ADS] [CrossRef] (In the text)
• Mac Low, M.-M., Klessen, R., & Burkert, A. 1998, Phys. Rev. Lett., 80, 2754 [NASA ADS] [CrossRef]
• Marchenko, S. V., Moffat, A. F. J., Ballereau, D., et al. 2003, ApJ, 596, 1295 [NASA ADS] [CrossRef]
• Matzner, C. D., & McKee, C. F. 1999, ApJ, 526, L109 [NASA ADS] [CrossRef] (In the text)
• Mignone, A. 2005, ApJ, 626, 373 [NASA ADS] [CrossRef]
• Mimica, P., Aloy, M. A., Müller, E., & Brinkmann, W. 2004, A&A, 418, 947 [NASA ADS] [CrossRef] [EDP Sciences]
• Myasnikov, A. V., & Zhekov, S. A. 1998, MNRAS, 300, 686 [NASA ADS]
• Norman, C. A.,, & Ferrara, A. 1996, ApJ, 467, 280 [NASA ADS] [CrossRef] (In the text)
• Nussbaumer, H., & Walder, R. 1993, A&A, 278, 209 [NASA ADS]
• Oskinova, L. M., Feldmeier, A., & Hamann, W.-R. 2004, A&A, 422, 675 [NASA ADS] [CrossRef] [EDP Sciences]
• Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef]
• Panaitescu, A., Spada, M., & Mészáros, P. 1999, ApJ, 522, L105 [NASA ADS] [CrossRef]
• Piran, T. 2004, Rev. Mod. Phys., 76, 1143 [NASA ADS] [CrossRef]
• Pittard, J. M., Dobson, M. S., Durisen, R. H., et al. 2005, A&A, 438, 11 [NASA ADS] [CrossRef] [EDP Sciences]
• Porter, D. H., & Woodward, P. R. 1994, ApJS, 93, 309 [NASA ADS] [CrossRef] (In the text)
• Porter, D., Pouquet, A., & Woodward, P. 1992, Theor. Comput. Fluid Dyn., 4, 13 [CrossRef]
• Porter, D., Pouquet, A., & Woodward, P. 1994, Phys. Fluids, 6, 2133 [NASA ADS] [CrossRef]
• Quirk, J. J. 1994, Internat. J. Numer. Methods Fluids, 18, 555
• Rees, M. J.,& Meszaros, P. 1994, ApJ, 430, L93 [NASA ADS] [CrossRef]
• Robinet, J.-C., Gressier, J., Casalis, G., & Moschetta, J.-M. 2000, J. Fluid Mech., 417, 237 [NASA ADS] [CrossRef] (In the text)
• Scalo, J., & Elmegreen, B. G. 2004, ARA&A, 42, 275 [NASA ADS] (In the text)
• Smith, M. D., Mac Low, M.-M., & Zuev, J. M. 2000, A&A, 356, 287 [NASA ADS] (In the text)
• Stevens, I. R., Blondin, J. M., & Pollock, A. M. T. 1992, ApJ, 386, 265 [NASA ADS] [CrossRef]
• Strickland, D., & Blondin, J. M. 1995, ApJ, 449, 727 [NASA ADS] [CrossRef] (In the text)
• Vázquez-Semadeni, E. 2004, Ap&SS, 292, 187 [NASA ADS]
• Vazquez-Semadent, E., Passot, T., & Pouquet, A. 1995, ApJ, 441, 702 [NASA ADS] [CrossRef]
• Vázquez-Semadeni, E., Gazol, A., & Scalo, J. 2000, ApJ, 540, 271 [NASA ADS] [CrossRef]
• Vázquez-Semadeni, E., Ryu, D., Passot, T., González, R. F., & Gazol, A. 2006, ApJ, 643, 245 [NASA ADS] [CrossRef]
• Vishniac, E. T. 1994, ApJ, 428, 186 [NASA ADS] [CrossRef] (In the text)
• Walder, R., & Folini, D. 1995, in IAU Symp., 163, 525
• Walder, R., & Folini, D. 1996, A&A, 315, 265 [NASA ADS]
• Walder, R., & Folini, D. 1998, A&A, 330, L21 [NASA ADS] (In the text)
• Walder, R., & Folini, D. 1999, in Hyperbolic Problems: Theory, Numerics, Applications, ed. M. Fey, & R.Jeltsch, 973 (In the text)
• Walder, R., & Folini, D. 2000a, in Thermal and Ionization Aspects of Flows from Hot Stars: Observations and Theory, ed. H. J. G. L. M. Lamers, & A. Sapar, ASP Conf. Ser., 281
• Walder, R., & Folini, D. 2000b, Ap&SS, 274, 343 [NASA ADS] (In the text)