Bow shock fragmentation driven by a thermal instability in laboratory astrophysics experiments

BOW SHOCK FRAGMENTATION DRIVEN BY A THERMAL INSTABILITY IN LABORATORY ASTROPHYSICS EXPERIMENTS F. Suzuki-Vidal 1 , S. V. Lebedev 1 , A. Ciardi 2,3 , L. A. Pickworth 1,7 , R. Rodriguez 4 , J. M. Gil 4 , G. Espinosa 4 , P. Hartigan 5 , G. F. Swadling 1,7 , J. Skidmore 1,8 , G. N. Hall 1,7 , M. Bennett 1 , S. N. Bland 1 , G. Burdiak 1 , P. de Grouchy 1,9 , J. Music 1 , L. Suttle 1 , E. Hansen 6 , and A. Frank 6 1 Blackett Laboratory, Imperial College London, Prince Consort Road, London SW7 2BW, UK; [email protected] 2 Sorbonne Universités, UPMC Univ. Paris 6, UMR 8112, LERMA, F-75005, Paris, France 3 LERMA, Observatoire de Paris, PSL Research University, CNRS, UMR 8112, F-75014, Paris, France 4 Departamento de Fisica de la Universidad de Las Palmas de Gran Canaria, E-35017 Las Palmas de Gran Canaria, Spain 5 Department of Physics and Astronomy, Rice University, 6100 S. Main, Houston, TX 77521-1892, USA 6 Department of Physics and Astronomy, University of Rochester, Rochester, NY 14627, USA Received 2015 September 14; accepted 2015 November 5; published 2015 December 14 ABSTRACT The role of radiative cooling during the evolution of a bow shock was studied in laboratory-astrophysics experiments that are scalable to bow shocks present in jets from young stellar objects. The laboratory bow shock is formed during the collision of two counterstreaming, supersonic plasma jets produced by an opposing pair of radial foil Z-pinches driven by the current pulse from the MAGPIE pulsed-power generator. The jets have different flow velocities in the laboratory frame, and the experiments are driven over many times the characteristic cooling timescale. The initially smooth bow shock rapidly develops small-scale nonuniformities over temporal and spatial scales that are consistent with a thermal instability triggered by strong radiative cooling in the shock. The growth of these perturbations eventually results in a global fragmentation of the bow shock front. The formation of a thermal instability is supported by analysis of the plasma cooling function calculated for the experimental conditions with the radiative packages ABAKO/RAPCAL. Key words: Herbig–Haro objects –instabilities –ISM: jets and outflows –methods: laboratory: atomic –plasmas – shock waves 1. INTRODUCTION One of the characteristic features of protostellar jets is the presence of shocks. They can be seen as large-scale terminal bow shocks or working surfaces (known as Herbig–Haro or HH objects)that form as the jet interacts with previous jet ejections. Smaller-scale, internal shocks are also present, which are driven by highly variable flow velocity in the jet. These are formed as the flow moves away from the protostar and encounters and overtakes slower material from previous ejections, with the process repeating along the jet beam. Proper-motion measurements of the jets in HH 34 (Reipurth et al. 2002)and HH 111 (Hartigan et al. 2001)show that the jet flow can reach peak velocities of ∼200–300 km s −1 with a typical velocity variation along the jet of ∼40 km s −1 . Shocks from protostellar jets exhibit complex dynamics in which different effects such as shear, hydrodynamic instabilities, and radiative cooling can be present simultaneously (Hartigan 2003). In this work we are particularly interested in the effect of radiative cooling, as it can drastically modify the shock morphology. Radiative losses are strongly dependent on the opacity (Drake et al. 2006), and if the shock region is optically thin, then radiation can escape the shock, leading to an increase in the postshock density. As the shock cools down, it can be prone to the growth of thermal instabilities (Field 1965; Hunter 1970),which can fragment and ultimately break up the shock. This effect has previously been studied mostly through numerical simulations (see, e.g., Blondin et al. 1989,1990; de Gouveia dal Pino & Benz 1993; Stone & Norman 1993; Frank et al. 1998;Teşileanu et al. 2008 and Asahina et al. 2014). In this paper we describe laboratory experiments that provide a complementary approach to study the effects of radiative cooling on the structure of the bow shocks. The similarity in the key dimensionless parameters characterizing our experiments, such as the Mach number and the cooling parameter, mean that the results are scalable to the internal shocks observed in young stellar object (YSO)jets. We observe the formation of a bow shock in the flow and its subsequent fragmentation, consistent with the onset of thermal instabilities. The bow shock in the experiments is formed from the interaction of two counterstreaming flows, which is equivalent to observations of internal shocks in YSO jets from a reference frame moving with the shock. An overall approach to modeling astrophysical phenomena in laboratory experiments is presented in the review by Remington et al. (2006), while a recent review of the synergy between observations, theory, and experiments relevant to the studies of YSO jets can be found in Frank et al. (2014). 2. EXPERIMENTAL SETUP The experimental setup is represented schematically in Figure 1. Two counterstreaming, supersonic plasma outflows are produced using two coaxial and oppositely facing radial foil Z-pinches (Suzuki-Vidal et al. 2009). Each radial foil is a metallic disk(40 mm diameter, 14 μm thick aluminum), subjected to a fast-rising electrical current pulse from the MAGPIE generator (Mitchell et al. 1996), which for the present experiments used a peak current of ∼1MA in ∼330 ns. The The Astrophysical Journal, 815:96 (9pp), 2015 December 20 doi:10.1088/0004-637X/815/2/96 © 2015. The American Astronomical Society. All rights reserved. 7 Present address: Lawrence Livermore National Laboratories, CA, USA. 8 Present address: AWE Aldermaston, UK. 9 Present address: Cornell University, NY, USA. 1 current is driven to the foils through 6.35 mm diameter stainless steel tubes touching each foil at its center, with the same current going through both of the foils via vertical posts (shown schematically in Figure 1(a)). The distance between the two foilsurfaces was ∼30 mm. The plasma is produced by continuous ablation of the surfaces of the foils as they are heated by the current, and the ablated plasma is accelerated by the axial pressure gradient produced by the current-induced azimuthal magnetic field as it diffuses through the foils. Each of the foils produces a supersonic plasma outflow that consists of a dense central jet surrounded by lower-density ambient plasma. Both components propagate with the same axial velocity of ∼50–100 km s −1 . Details on the formation of the outflows by a single radial foil Z-pinch can be found in Suzuki-Vidal et al. (2009), Ciardi et al. (2009), Gourdain et al. (2010), and SuzukiVidal et al. (2012). In the experiments presented here, the interaction of the jets was diagnosed from the side-on (radial) direction using an optical framing camera (Invisible Vision UHSi 12/24)that imaged the optical self-emission from the plasma. This camera is capable of taking up to 12 images per experiment, with 5 ns exposure and 30 ns interframe separation. We also used simultaneous optical laser shadowgraphy and interferometry (λ=532 nm, pulse duration 0.3 ns). The latter diagnostic was used to measure the electron density distribution of the different plasma features present, i.e., jets, surrounding plasma, and bow shock (Swadling et al. 2014). 2.1. Scaling Jets from YSOs are generally well described by ideal magnetohydrodynamics (MHD), and the experiments presented here are designed to produce flows in a similar regime. It is evident that these evolve over hugely different length scales and timescalesand have different physical characteristics in terms of density, temperature, and chemical composition. Nevertheless, the invariance properties of the ideal MHD equations provide a framework that allows meaningful scaling of the jet dynamics over many orders of magnitude (see, e.g., Ryutov et al. 1999,2000,2001), provided that certain constraints between the flow variables are satisfied. It is important to note that such invariance is also applicable to flows with shocks and to radiative flows under a more restrictive set of constraints (Falize et al. 2011). In order for a fluid description to be applicable, we require the flows to have a localization parameter δ=1. This is equivalent to an ion mean free path much less than the characteristic spatial scale of the system. Using the parameters from Table1 in Hartigan et al. (2009), we estimate that jets from YSOs have δ∼10 −5 . Applicability of an MHD description also requires that the transport of momentum, magnetic field, and thermal energy occurs predominantly through advection with the flow, i.e., negligible dissipation through viscosity (Reynolds number Re ?1), magnetic diffusivity (magnetic Reynolds number Re M ?1), and heat conduction (Peclet number Pe ?1), respectively. Owing to their large spatial scales, YSO jets are characterized by Re, Re M ,Pe10 7 . Our experiments are characterized by δ∼10 −4 , Re ∼10 5 ,Pe∼10 3 , and Re M ∼10 3 (Suzuki-Vidal et al. 2012), and thus we expect a similar overall physical behavior. Additionally, the interaction of radiative jets with the interstellar medium can be broadly classified by three dimensionless parameters (Blondin et al. 1990): the Mach number M(the ratio of flow speed to sound speed), the density contrast η(the ratio of density between the flow and the ambient medium where it propagates), and the cooling parameter χ cool (the ratio of the cooling time τ cool to the characteristic hydrodynamical time of the flow τ hydro ), which quantifies the effect of radiative losses in the plasma. The specific values of M,η, and χ cool essentially determine the overall morphology of the flow. Internal bow shocks in YSO jets are characterized by a relative velocity of ∼40 km s −1 and M∼10, and we can expect η∼1 and χ cool 1(Hartigan et al. 2009). As will be discussed later in the paper, the values of these three parameters in our experiments are close to those in YSO jets. Overall, the similarity of the dimensionless parameters allows applying the Eulerian scaling relations described in detail in Ryutov et al. (1999,2000,2001). Flows with identical Mach numbers will evolve with identical morphology, but on different hydrodynamic temporal (τ hydro )and spatial (r jet ) scales, related via the corresponding flow velocities (V flow )as τ hydro =r jet /V flow . Taking the jet radius as a characteristic spatial scale, young stellar jets typically have r YSO ∼50 AU ∼10 15 mm and flow velocities of V YSO =40 km s −1 (i.e., typical velocity variability in YSO jets). For the experiments here r exp ∼1.5 mm and V exp =140 km s −1 (i.e., the relative velocity of a single jet in our experiments;see explanation in the next section). Thus, the characteristic temporal scales are τ YSO =1.9 ×10 8 s(∼6yr) and τ exp =10 ns, respectively. The total time interval over which the evolution of the flows is followed in the experiments of ∼300 ns (≈30τ exp )thus corresponds to ∼180 yrof evolution for the astrophysical counterpart, which exceeds the typical timescale for multi-epoch observations of YSOs with Figure 1. Schematic experimental configuration to study the formation of a bow shock from the interaction between two counterstreaming jets with different relative axial velocities, represented as opposite vertical arrows on axis. The schematic depicts a side-on (radial), cut view of the system, which has azimuthal symmetry. The dashed (red)arrows represent the path of the current that drives the two plasma flows. The (blue)arrows pointing into and out of the page correspond to the azimuthal magnetic field generated by the current, which provides the driving force for the two outflows. The jets are surrounded bylower-density plasma (yellow regions), which moves with the same axial velocity as the jets. Smaller arrows in these regions represent the plasma flow direction. The images depict the two counterstreaming outflows (a)before their collision and(b)after they collide, triggering the formation of a bow shock moving toward the bottom foil. 2 The Astrophysical Journal, 815:96 (9pp), 2015 December 20 Suzuki-Vidal et al. the Hubble Space Telescope(HST)of ∼10 yr(see, e.g., Hartigan et al. 2011). 3. EXPERIMENTAL RESULTS AND DISCUSSION The overall evolution of the interaction of the two counterstreaming plasma jets can be seen in Figure 2, which shows side-on optical emission images obtained in the same experiment at different times after the start of the current pulse driving the jets. Images in the top row of Figure 2show the entire region between the two foils and correspond to the early times of the interaction (310–430 ns). The images show the formation of two well-collimated jets on the axis of the system propagating toward each other. Although the top and bottom foils are nominally the same and are driven by the same current, an asymmetry between the top and bottom flows is evident from the emission images. The top jet is formed first compared to the bottom jet, which is seen by looking at the positions of the tips of the jets indicated by horizontal arrows on the first three frames. This asymmetry is caused by oppositepolarity current drive for each foil (radially outwardon the bottom foil compared to radially inwardon the top foil)and is reproducible from experimenttoexperiment. The asymmetry affects the details of the plasma ablation and initial acceleration of the flow at the foil (Gourdain & Seyler 2013). This results in different flow parameters (e.g., ram pressures), which, as the jet collide, lead to the formation of a bow shock. This is first evidenced at 370 ns as a highly emitting region aligned with the head-on collision between the two jet tips. The bow shock is fully formed at 400 ns and is seen to move downward, i.e., toward the bottom foil. At this time the postshock region is seen as a highly compressed, highly collimated column, while the bow shock begins to develop small-scale structures from here onward. The long-term evolution of the bow shock and of these structures is highlighted in the subsequent images, corresponding to times 460–640 ns. The first collision between the two jets seen in the image at 370 ns occurs ∼2 mm below the midplane, which suggests that the two jets have slightly different velocities. This is confirmed Figure 2. Counterstreaming jet interaction results from optical self-emission of the plasma obtained from the same experiment. The arrows in the first three frames indicate the position of the tip of both jets (visibility dependent on image contrast levels), with their collision highlighted at 370 ns. The last six frames are focused on the bow shock region, which is seen to fragment most evidently in the last three times. 3 The Astrophysical Journal, 815:96 (9pp), 2015 December 20 Suzuki-Vidal et al. by following the positions of the visible tips of the jets in time before their collision, resulting in tip velocities of V topjet ∼80 ±10 km s −1 and V bottomjet ∼60 ±10 km s −1 . The higher velocity in the top jet, combined with a larger density inferred from the stronger self-emission and from laser interferometry, impliesan imbalanced ram pressure (P ram =ρV 2 )in the collision, which explains the shape and orientation of the bow shock and its downward propagation. The velocity of the leading edge of the bow shock was measured as V bow ∼40 ±10 km s −1 . The velocity of the bow shock measured in the experiments can be compared with the velocity expected from a onedimensional, momentum flux conservation argument (Norman et al. 1983; Hartigan 1989; Blondin et al. 1990; de Gouveia dal Pino & Benz 1994; de Gouveia Dal Pino 2005; Nicolaï et al. 2008). This is convenient to do in the reference frame of a single jet driving a shock through a stationary medium ahead of it. The velocity of a bow shock V bow is related to the velocity of the jet V jet and the density contrast ηbetween the jet and the preshock external medium by V bow ≈V jet (1+η −1/2 ) −1 . The measured experimental jet tip velocities of V jettop ∼80 km s −1 and V jetbottom ∼60 km s −1 for the top and bottom jets, respectively, can be translated into a single jet with a relative tip velocity of V jetrel =140 km s −1 , which in turns forms a bow shock with a relative velocity of V bowrel =100 km s −1 . Using these relativevelocities in the equation above results in a density contrast η∼6, i.e., the experiment is equivalent to the interaction of a single jet that is six times denser than the ambient medium where the shock is generated; for this particular counterstreaming geometry, the ambient medium is the opposite (bottom)jet. We can calculate the internal Mach number of the jet working surface by taking the ratio of the jet relative velocity (V jetrel =140 km s −1 )and the typical ion acoustic speed in the jet flow of c s ∼15 km s −1 (Suzuki-Vidal et al. 2012), resulting in an internal Mach number of M∼10. Therefore, the values of ηand Min the experiments are similar to those in YSO jets (Hartigan et al. 2009), and thus the experiments should evolve with a similar jet/shock morphology as discussed earlier. In addition to the bow shock, the images in Figure 2also show the collision between the two flows off-axis, i.e., the plasma surrounding the jets, which leads to the formation of a double-shock structure that extends at large radii. This feature is especially evident at the top of the images from ∼460 ns onward. These standing shocks remain approximately symmetric with respect to the midplane between the two foils. Preliminary MHD simulations indicate that these features arise as a result of the presence of toroidal magnetic field advected by the counterstreaming flows (Suzuki-Vidal et al. 2014).A more detailed discussion of the pileup of the magnetic fields will be published separately. 3.1. Bow Shock and Jet Working Surface: Evolution and Small-scale Structures The emission images in Figure 2show in detail the evolution of the bow shock from its formation at 400 ns through its rapid fragmentation. At 400 ns the surface of the bow shock and the postshock region just behind it are smooth. The formation of small-scale structures occurs on timescales that are shorter than the interframe separation for this particular diagnostic. This gives an upper limit for their development timescale of ∼30 ns, which is consistent with the appearance of new structures between frames at later times. The characteristic spatial scale of the structures seen in the images ranges from ∼200 μmto ∼2 mm, and at the final stages of the evolution (the last two panels at 610 and 640 ns)the bow shock has fragmented into a large number of emitting clumps. The smallest observed size (∼200 μm)seen on the images obtained with the optical framing camera is comparable to the spatial resolution of this diagnosticand is also comparable to the motional blurring due to the temporal resolution of the diagnostic (5ns×40 km s −1 =200 μm). Further details of the small-scale structures present in the bow shock and postshock region were obtained with laser probing (0.3 ns pulse duration), which has significantly better spatial resolution (50 μm)and reduced motional blurring (∼10 μm). The interferometry channel of this diagnostic provides measurements of the spatial distribution of the plasma electron density, and the shadowgraphy channel, which is sensitive to spatial gradients of the electron density, provides information on the characteristic spatial scales of the nonuniformities. Figure 3shows laser probing images obtained at 400 ns in the same experiment as the optical emission images shown in Figure 2, i.e., the timing of the laser probing images corresponds to the fourth panel in Figure 2. Figure 3(a)shows the raw image from the interferometry channel. The apparent distortion of the interference fringes, which were initially horizontal and uniformly spaced, is caused by changes in the interference state that are induced by a spatially varying phase delay imparted on the probe beam by the plasma. This phase delay can be extracted from the image (Swadling et al. 2014) and is proportional to the electron column density of the plasma, n e L(in cm −2 ). Here Lis the length of plasma along the probing beam, which changes as a function of position in the plane of the image. The sufficiently good axial symmetry of the object allows applying Abel inversion (Hutchinson 2005)to the electron column density, resulting in the axisymmetric (radial) distribution of electron density n e (r)(in cm −3 )shown in Figure 3(b). The highest electron density near the axis, reaching n e ∼10 19 cm −3 (up to 75% uncertainty due to errors in the choice of the central axis and left-right asymmetries),is observed in the narrow, compressed region formed above the bow shock (the postshock region), while the electron density in the flow below the bow shock (the preshock region)is n e ∼(3±2)×10 18 cm −3 (i.e., ∼60% uncertainty). We note that, although these errors seem large, for measurements of electron density off-axis these errors decrease rapidly and can reach values 20%. The higher electron density present in the jet driving the bow shock is qualitatively consistent with the density contrast (heavy jet propagating through ambient medium)inferred from application of the one-dimensional momentum flux conservation argument introduced in the previous section. The laser shadowgraphy diagnostic in Figure 3(c)shows several interesting features. First, there are two dark regions in the top part of the image positioned on both sides of the axis, and just below them a dark vertical region on axis. These features are produced by large density gradients, and their shape is consistent with compression of the jet to a diameter smaller than the initial jet diameter due to converging flows driven by large postshock pressures. Formation of such shocks was observed in numerical simulations of high Mach number jets propagating through ambient media (e.g., Norman et al. 1982), and the structure we see in the experiment is 4 The Astrophysical Journal, 815:96 (9pp), 2015 December 20 Suzuki-Vidal et al. very similar to that observed in simulations of radiatively cooled jets with a similar density contrast, e.g., presented in Figure7(c)of Blondin et al. (1990)for a density contrast of η=3, which is of the same order to the one inferred for these experiments of η=6. The simulations by Blondin et al. (1990) also show the presence of a second standing conical shock, positioned between the bow shock and the region of maximum compression that opens up toward the bow shock. The image in Figure 3(c)indicates that this conical shock is also present in the experiment, though it is less pronounced, indicating a smaller density gradient than that in the downstream conical shock. We also note the presence of a horizontal dark region connecting the dark features at the top of the image, which is consistent with the expected shape of a Mach disk (see, e.g., Figure1(a)in Hartigan 1989 and Figure7(b)in Blondin et al. 1990). Finally, the image in Figure 3(c)shows the presence of small-scale structures already developing behind the bow shock at this early stage of the evolution. The perturbations appear to be elongated, predominantly in the direction normal to the jet axis, and the smallest detectable spatial scale is ∼120–170 μm, which is a factor of ∼3 larger than the spatial resolution of the diagnostic (∼50 μm in the diffraction-limited case).We interpret the observed development of small-scale perturbations in the bow shock leading to its complete fragmentation later in time as being a result of strong radiative cooling in the postshock plasma. The fragmentation of the bow shock and the development of dense clumps in this region are associated with a thermal instability (Field 1965), which is discussed thoroughly in the next section. The fragmentation of shocks in protostellar jets due to thermal instabilities has been previously studied by numerical simulations (see, e.g., Blondin et al. 1989,1990; Frank et al. 1998);however, the effect of this instability in YSO jets and shocks is still unclear. 3.2. Radiative Cooling and Thermal Instabilities in the Bow Shock The most striking result in the experiments is the rapid development of small-scale spatial structures in the bow shock and postshock region. Our interpretation is that the fragmentation of the bow shock is related to a dynamic, local thermal instability that leads to the condensation of density perturbations by radiative cooling (see, e.g., Field 1965; Mathews & Bregman 1978; Fall & Rees 1985; Balbus 1986 and Blondin & Cioffi1989). The instability develops over the characteristic radiative cooling time, τ cool , which for optically thin plasmas, such as the ones in our experiments (Espinosa et al. 2015),is given by the ratio of thermal energy density, U, to the radiated power per unit volume P rad =n e n i Λ(n i ,T e )as τ cool =U/n e n i Λ (n i ,T e ), where Λ(n i ,T e )is the normalized cooling function (in erg cm 3 s −1 )and n i is the ion density. This expression can also be written as (Ryutov et al. 1999) ZT Zn n T s2.410 1eV cm , ,1 cool 12 i3ie () () [] ¯[] ¯() ⎡ ⎣⎤ ⎦ t=´ + L - - where Z ¯is the average ionization in the plasma. When radiative cooling occurs on timescales that are long compared to the sound speed crossing time, perturbations over a region of size λ iso =c s τ cool (where c s is the ion acoustic speed)tend to maintain a common pressure (isobaric). In this case nT ∼constant and the cooling rate then scales as TT. cool 12 ()t G =µL -In the isobaric regime, if the cooling rate increases for a decreasing temperature, i.e., dΓ/dT <0, then radiation losses will be even more efficient in removing energy from the plasma and further reducing its temperature. To maintain pressure equilibrium with its surroundings,the density increases, thus further increasing the radiated power losses and potentially leading to a runaway condensation of the initial density perturbations. Detailed analysis of the instability Figure 3. Optical laser probing of the bow shock at 400 ns (same experiment as Figure 2).(a)Raw data from laser interferometry;(b)axisymmetric electron density (n e )from analysis of (a);(c)laser shadowgraphy, with the field of view shown schematically in the dashed inset in (a)and (b). The fringing on this image is an artifact present in this particular diagnostic. 5 The Astrophysical Journal, 815:96 (9pp), 2015 December 20 Suzuki-Vidal et al. was done by Field (1965), Hunter (1970), and Balbus (1986). For isobaric modes, the condition for stability is given by dnT dT log , log 2. 2 i () () () () L < While λ iso sets an upper limit for the the lengthscale of isobaric condensation, thermal conduction will suppress shortwavelength perturbations that are of the order of the so-called Fieldʼs length ( 21 Field th cool 12 [( ) ] l pg ct=- ), where χ th =κ/n e is the thermal diffusivity and γis the ratio of heat capacities. The most unstable wavelength, valid in the regime of large wavenumbers (but still smaller than a critical wavenumber k 2 crit Field 1 pl=-), is given by the geometric mean (Field 1965)max Field iso 12 () l ll=so that density and temperature perturbations will be unstable in the wavelength range .3 Field iso ()lll<< Although cooling in our experiments is different from that in astrophysics owing to the differences in elements and physical conditions of temperature and densities, for optically thin plasmas the relevant parameters for comparison are the dimensionless cooling parameter χ cool and the dependence of the cooling function Λ(n i ,T)on temperature. To estimate the effects of radiative cooling in the postshock aluminum plasma, we use new cooling rates calculated with the computational packages ABAKO/RAPCAL (Rodriguez et al. 2008; Florido et al. 2009). The codes calculate the plasma level populations and average ionizations by solving the set of rate equations of the collisional-radiative model implemented in ABAKO (assuming that the plasma is optically thin and in steadystate). The model is capable of accounting for coronal equilibrium andlocaland nonlocal thermodynamic equilibrium regimes. The atomic data required were obtained using the FAC code (Gu 2008)in the relativistic detailed configuration accounting approach, with the spin–orbit split arrayformalism (BaucheArnoult et al. 1985)and including configuration interaction within the same nonrelativistic atomic configurations. The databases of cooling functions and average ionizations obtained with ABAKO/RAPCAL were subsequently parameterized as a function of the plasma density and temperature using the PARPRA code (Rodriguez et al. 2014). Results from these numerical calculations are presented in Figures 4(a)–(b), showing the variation of the cooling function Λ(n i ,T e )and the average ionization Z ¯for an aluminum plasma as a function of electron temperature (in the range of ∼10 4 –10 6 K, Figure 4. Calculated (a)cooling function and (b)average ionization for aluminum as a function of electron temperature for different ion densities relevant to our experiments (the legend shown in (a)applies to the entire figure).(c)Cooling time calculated using Equation (1).(d)Thermal instability analysis applied to the cooling functions shown in (a). The regions below each horizontal line indicate the thresholds for the onset of isobaric (<2), isochoric (<1), and adiabatic (<−0.5)thermal instabilities. 6 The Astrophysical Journal, 815:96 (9pp), 2015 December 20 Suzuki-Vidal et al. ∼1–85 eV)for different ion densities typical of our experiments (n i =10 17 –10 20 cm −3 ). Figure 4(a)shows the rate of cooling increases with decreasing ion density, varying up to two orders of magnitude. The average ionization in Figure 4(b) shows an increase with temperature, with little variation as a function of ion density. With the cooling function and the average ionization we can calculate the cooling time τ cool (Equation (1)), and this is plotted for different aluminum ion densities in Figure 4(c). The plot shows that overall the cooling time decreases with increasing ion density, as τ cool ∝ 1/Λ(n i ,T e ), and longer cooling times can be expected around an electron temperature of T e ∼10 5 K. With the cooling time we can calculate the cooling parameter χ cool =τ cool /τ hydro , which quantifies the importance of radiative losses in the plasma. Taking the estimated hydrodynamical time of the experiments of τ exp ∼10 ns (see Section 2.1), we see that, independent of ion density,τ cool >τ exp in the region around T e ∼10 5 K, which implies χ cool 1 and thus thatradiative losses are not important. This region correlates with a dip in the cooling functions at this temperature. Besides this temperature region, however, the cooling time is overall of the order of the hydrodynamical time, and thus the cooling parameter χ cool ∼1. The hydrodynamical time can also be estimated as the time it takes for the shock to form and become fully disrupted, which from Figure 2can be taken as τ exp ∼550−370 ns =180 ns, resulting in χ cool ∼(3–30 ns)/180 ns <1, and thus the plasma in the shock is expected to be radiatively cooled. Knowledge of the dependence of the cooling function as a function of temperature allows performing a thermal instability analysis. The stability condition for isobaric modes presented in Equation (2)(dlogΛ(n i ,T e )/dlogT e <2)is plotted in Figure 4(d)using the cooling functions for different ion densities in Figure 4(a). The threshold for isobaric instabilities shows that, almost independently of the ion density, the plasma is expected to become unstable at temperatures of T e ∼(1.5–9)×10 4 K(∼5eV)and at T e 2.5 ×10 5 K (20 eV). Following the analysis in Shchekinov (1978),we also plot other instability thresholds, namely, isochoric and adiabatic modes, which correspond to dlogΛ(n i ,T e )/dlogT e <1 and <−0.5, respectively. Furthermore, the onset of thermal instabilities in the postshock plasma depends on the thermal evolution of the ion and electron plasma components. Because of the large difference between the ion and electron masses, there is a large temperature difference between the ion and electron components immediately behind the bow shock. In the strong-shock approximation (Mach number M?1), the postshock ion temperature is TAmVk21 110K, Bi,PS p bow rel 226 () ()gg=- +~ where Ais the atomic weight (A=27 for aluminum),m p is the proton mass, V bowrel is the shock velocity in the reference frame of a stationary preshock medium (i.e., V bowrel =100 km s −1 ),and γ=5/3 assuming an ideal gas. Compared to the ions, the electrons are compressed adiabatically across the shock and their temperature increases only by a factor of 4 γ−1 ∼2.5 to T e, PS ∼(2–3)×10 5 K. Although the plasma can be unstable at these temperatures, the existence of a well-defined unstable range of wavelengths, as given in Equation (3), is not satisfied. The postshock is characterized by a relaxation region where the energy exchange between ions and electrons competes with radiative cooling. Immediately after the shock, the electrons are thermally decoupled from the ions and their cooling timescale (1ns)is much shorter than the ion–electron energy equilibration timescale (Z20 ns ie ei  tt=~ ). Thus, over timescales of a nanosecond, the electron temperature rapidly drops while the ion temperature remains essentially constant. The cooling rate also decreases considerably, and an equilibrium is immediately reached behind the shock, where electron heating due to the energy received by the ions balances their radiative energy losses TT n k nT1,. 4 B ie ei iie () () ()  tg -»- L In this regime 15 20 ns ie cool  tt»~- and the electrons are approximately isothermal. Detailed calculations of temperature equilibration show that the equilibrium electron temperature in the postshock relaxation layer is T e ∼(1–1.5)×10 5 K;thus, the plasma is expected to be thermally stable. After ∼40 ns the electron and ion plasma components are fully equilibrated and their common temperature decreases below 7×10 4 K. This places the plasma in conditions corresponding to a thermally unstable regime. The condition on the wavelengths is satisfied, and density and temperature perturbations with wavelengths λ∼30–100 μm grow over very short timescales of the order of 1 ns. 4. CONCLUSIONS We presented a new experimental configuration that aims at studying experimentally the formation of bow shocks relevant to those present in young stellar jets. In our experiments we produce two counterpropagating, supersonic plasma jets from plasma ablation of two radial foil Z-pinches. A bow shock is driven by the head-on collision between the two jets, and the bow shock properties are determined by their relative velocities and densities. A key result is that initially the bow shock is smoothbut then quickly, within timescales of ∼30 ns, develops small-scale spatial features behind the shock front. This timescale is consistent with the expected cooling times in the experiments, which were estimated using the cooling functions and average ionization calculated with the radiative packages ABAKO/RAPCAL. This allowed performing a thermal instability analysis for isobaric modes, resulting in expected temperature ranges at which the plasma should become thermally unstable. Detailed analytical calculations of the thermal evolution of the ion and electron components predict, for the experimental shock conditions, the development of an isobaric thermal instability in timescales of ∼40 nsand with typical spatial scales of ∼30–100 μm. This is in very good agreement with our experimental results that show an initially smooth bow shock that, within ∼30 ns, fragments into smallscale spatial features with typical sizes of ∼120 μm. Thermal instabilities may play a role in HH shocks provided that the shock velocities are high enough to raise the temperature in the postshock gas significantly above the peak of the cooling curve around 2 ×10 5 K. The critical shock velocity for the onset of cooling instabilities has been estimated to be ∼200 km s −1 in 1D simulations (Smith 1989), though Sutherland et al. (2003)found that thermal instabilities in their 2D simulations generated filaments and voids in the postshock gas for shock velocities as low as 120 km s −1 . However, Innes (1992)demonstrated that even a fairly weak preshock magnetic field stabilized shocks up to 175 km s −1 . In HH jets, most shocks have velocities 100 km s −1 and should not be affected by thermal instabilities. However, the strongest bow shocks, such as HH1 and HH2, have high ionization lines (Boehm 7 The Astrophysical Journal, 815:96 (9pp), 2015 December 20 Suzuki-Vidal et al. et al. 1993)and broad line profiles (Hartigan et al. 1987) indicative of shock velocities ∼200 km s −1 , well above the criteria for thermal instabilities. As shown in Figure 5, bright knots do appear along the bow shock in HH1 on timescales of decades. Using the scaling presented at the end of Section 2.1, the laboratory temporal scale of 30ns would correspond to an astrophysical timescale of ∼17 yr, consistent with the astronomical observations. The spatial scales associated with the astronomical knots also scale well with the nonuniformities in the experiment. For example, the minimum size of the nonuniformities in the experiments of 0.12 mm scales to 4 AU in the astronomical observations, where the spatial resolution is ∼20 AU. The larger nonuniformities present in the optical self-emission (e.g., the last three panels of Figure 2) have a typical size of ∼1 mm, which scales to 30 AU in Figure 5, consistent with the size of the new knots in HH1. Driven from the same source on the other side of the outflow, HH2 also has a high shock velocity and shows small knots that appear and merge along the strongest shock fronts (Hartigan et al. 2011). Of course, there are other ways to generate clumps along bow shocks, such as Kelvin–Helmholtz instabilities or even a clumpy preshock density, and we cannot rule out these possibilities without further observational data. More generally, the thermal instability analysis could be applied to cooling curves used for numerical simulations of protostellar jets and shocks (see, e.g., Dalgarno & McCray 1972; Kafatos 1973; Sutherland & Dopita 1993). Although in these cases we expect the temperature ranges for the onset of thermal instabilities to differ from those expected in the experiments owing to the different elements and abundances that characterize them (e.g., H in simulations compared to Al in the experiments), in both cases the instability should develop at the appropriate slope of the cooling curve by the condition given in Equation (2). The experiments presented here are important as they shed light on the onset and nonlinear evolution of the thermal instability. Numerical simulations of thermally unstable flows can be challenging as thermal conduction must be explicitly included to suppress the growth of small-scale perturbations. In particular, the Fieldʼs length has to be resolved by at least a few computational cells (Koyama & Inutsuka 2004)to avoid the growth of perturbation and fragmentation at the grid scale. Furthermore, the presence of magnetic fields further increases the complexity of numerical calculations by making thermal conduction anisotropic. In that direction, similar experiments to the ones presented here can be designed to increase sufficiently the magnetic field in the flow to allow studying such a regime in the future. This could be relevant to previous theoretical studies that predict the suppression of thermal instabilities by magnetic fields (e.g., Innes 1992; Lesaffre et al. 2004). This work was supported in part by the Royal Society through a University Research Fellowship, by EPSRC Grant No. EP/G001324/1, by DOE cooperative agreements No. DEF03-02NA00057 and No. DE-SC-0001063, and by the LABEX Plas@Par project. This work received financial state aid managed by the Agence Nationale de la Recherche, as part of the program “Investissements d’avenir”under ref. ANR-11IDEX-0004-02. REFERENCES Asahina, Y., Ogawa, T., Kawashima, T., et al. 2014, ApJ,789, 79 Balbus, S. A. 1986, ApJL,303, L79 Bauche-Arnoult, C., Bauche, J., & Klapisch, M. 1985, PhRvA,31, 2248 Blondin, J. M., & Cioffi, D. F. 1989, ApJ,345, 853 Blondin, J. M., Fryxell, B. A., & Konigl, A. 1990, ApJ,360, 370 Blondin, J. M., Konigl, A., & Fryxell, B. A. 1989, ApJL,337, L37 Boehm, K.-H., Noriega-Crespo, A., & Solf, J. 1993, ApJ,416, 647 Ciardi, A., Lebedev, S. V., Frank, A., et al. 2009, ApJL,691, L147 Dalgarno, A., & McCray, R. A. 1972, ARA&A,10, 375 de Gouveia Dal Pino, E. M. 2005, AdSpR,35, 908 de Gouveia dal Pino, E. M., & Benz, W. 1993, ApJ,410, 686 de Gouveia dal Pino, E. M., & Benz, W. 1994, ApJ,435, 261 Drake, R. P., Davison, L., & Horie, Y. 2006, High-Energy-Density Physics: Fundamentals, Inertial Fusion, and Experimental Astrophysics (Berlin: Springer) Espinosa, G., Gil, J., Rodriguez, R., et al. 2015, HEDP,17, 74 Falize, É., Michaut, C., & Bouquet, S. 2011, ApJ,730, 96 Fall, S. M., & Rees, M. J. 1985, ApJ,298, 18 Field, G. B. 1965, ApJ,142, 531 Florido, R., Rodríguez, R., Gil, J. M., et al. 2009, PhRvE,80, 056402 Frank, A., Ryu, D., Jones, T. W., & Noriega-Crespo, A. 1998, ApJL,494, L79 Frank, A., Ray, T. P., Cabrit, S., et al. 2014, Protostars and Planets VI (Tuscon, AZ: Univ. Arizona Press) Gourdain, P.-A., Blesener, I. C., Greenly, J. B., et al. 2010, PhPl,17, 012706 Gourdain, P.-A., & Seyler, C. E. 2013, PhRvL,110, 015002 Gu, M. F. 2008, CaJPh,86, 675 Hartigan, P. 1989, ApJ,339, 987 Hartigan, P. 2003, Ap&SS,287, 111 Hartigan, P., Foster, J. M., Wilde, B. H., et al. 2009, ApJ,705, 1073 Hartigan, P., Frank, A., Foster, J. M., et al. 2011, ApJ,736, 29 Hartigan, P., Morse, J. A., Reipurth, B., Heathcote, S., & Bally, J. 2001, ApJL, 559, L157 Hartigan, P., Raymond, J., & Hartmann, L. 1987, ApJ,316, 323 Hunter, J. H., Jr. 1970, ApJ,161, 451 Hutchinson, I. H. 2005, Principles of Plasma Diagnostics (Cambridge: Cambridge Univ. Press) Figure 5. Multiepoch observations of HH1 showing the formation of knots in the bow shock (adapted from Figure2 in Hartigan et al. 2011). The image shows two out of three epochs of HST images, with Hαin green, [SII]in red, and yellow denotingemission in both filters. 8 The Astrophysical Journal, 815:96 (9pp), 2015 December 20 Suzuki-Vidal et al. Innes, D. E. 1992, A&A, 256, 660 Kafatos, M. 1973, ApJ,182, 433 Koyama, H., & Inutsuka, S.-i 2004, ApJL,602, L25 Lesaffre, P., Chièze, J.-P., Cabrit, S., & Pineau des Forêts, G. 2004, A&A, 427, 147 Mathews, W. G., & Bregman, J. N. 1978, ApJ,224, 308 Mitchell, I. H., Bayley, J. M., Chittenden, J. P., et al. 1996, RScI,67, 1533 Nicolaï, P., Stenz, C., Kasperczuk, A., et al. 2008, PhPl,15, 082701 Norman, M. L., Winkler, K.-H., & Smarr, L. 1983, BAAS, 15, 962 Norman, M. L., Winkler, K.-H. A., Smarr, L., & Smith, M. D. 1982, A&A, 113, 285 Reipurth, B., Heathcote, S., Morse, J., Hartigan, P., & Bally, J. 2002, AJ, 123, 362 Remington, B. A., Drake, R. P., & Ryutov, D. D. 2006, RvMP,78, 755 Rodriguez, R., Espinosa, G., Gil, J., et al. 2014, CCoPh, 16, 612 Rodriguez, R., Florido, R., Gil, J., et al. 2008, LPB, 26, 433 Ryutov, D., Drake, R. P., Kane, J., et al. 1999, ApJ,518, 821 Ryutov, D. D., Drake, R. P., & Remington, B. A. 2000, ApJS,127, 465 Ryutov, D. D., Remington, B. A., Robey, H. F., & Drake, R. P. 2001, PhPl, 8, 1804 Shchekinov, I. A. 1978, SvA, 22, 182 Smith, M. D. 1989, MNRAS,238, 235 Stone, J. M., & Norman, M. L. 1993, ApJ,413, 198 Sutherland, R. S., Bicknell, G. V., & Dopita, M. A. 2003, ApJ,591, 238 Sutherland, R. S., & Dopita, M. A. 1993, ApJS,88, 253 Suzuki-Vidal, F., Bocchi, M., Lebedev, S. V., et al. 2012, PhPl,19, 022708 Suzuki-Vidal, F., Lebedev, S. V., Ciardi, A., et al. 2009, Ap&SS,322, 19 Suzuki-Vidal, F., Lebedev, S., Pickworth, L. A., et al. 2014, in APS Meeting Abstracts, Session J05: HEDP Laboratory Physics (Ridge, NY: APS), 00005 Swadling, G. F., Lebedev, S. V., Hall, G. N., et al. 2014, RScI, 85, 110000 Teşileanu, O., Mignone, A., Massaglia, S., et al. 2008, A&A,488, 429 9 The Astrophysical Journal, 815:96 (9pp), 2015 December 20 Suzuki-Vidal et al.