Cooling and Instabilities in Colliding Radiative Flows with Toroidal Magnetic Fields

We report on the results of a simulation based study of colliding magnetized plasma flows. Our set-up mimics pulsed power laboratory astrophysical experiments but, with an appropriate frame change, are relevant to astrophysical jets with internal velocity variations. We track the evolution of the interaction region where the two flows collide. Cooling via radiative loses are included in the calculation. We systematically vary plasma beta ($\beta_m$) in the flows, the strength of the cooling ($\Lambda_0$) and the exponent ($\alpha$) of temperature-dependence of the cooling function. We find that for strong magnetic fields a counter-propagating jet called a"spine"is driven by pressure from shocked toroidal fields. The spines eventually become unstable and break apart. We demonstrate how formation and evolution of the spines depends on initial flow parameters and provide a simple analytic model that captures the basic features of the flow.


INTRODUCTION
Shocks arise in hypersonic plasma flows when faster moving material collides with slower material ahead of it.Radiative shocks, where energy is lost to optically thin photons, can arise in a wide variety of physical settings, such as protostellar jets (Raga & Kofman 1992;Frank et al. 2014), supernova explosions (Wheeler et al. 2002), gamma ray bursts (Piran 2004), active galactic nuclei (Rees 1978), and High Energy Density Plasma (HEDP).Hydrodynamic and MHD evolution in collision zones is expected to lead to high degrees of heterogenity or "clumpiness" (Hansen et al. 2017).In the context of protostellar jets, for example, flow collisions correspond observationally to chains of knots (Hartigan et al. 2011).
Supersonic jets are structurally composed of a supersonic beam, a cocoon of shocked jet gas, a region of shocked ambient gas, and a bow shock (Blondin et al. 1990), as shown in figure 1.When two supersonic flows collide, an interaction region forms, bounded by a pair of jet shocks also known as working surfaces.If radiative losses are significant, the interaction region will include a cooling region, in which shocked gas cools from its initial post-shock temperature Ts.Further behind the shock, gas reaches its final post-cooling temperature T f and is deposited onto the cold slab, where densities are highest.Some shocked material gets ejected laterally by the high pressures throughout the interaction region, producing shocked lateral outflows (SLOs) (Falle & Raga 1993;Markwick et al. 2022).
When magnetic fields are present, additional structures can form within the interaction region.One-dimensional simulations of pulsed jets conducted by Hartigan et al. (2007) found that colliding magnetized blobs can essentially 'bounce' off of one another as magnetic pressure increases in cooling zones behind the shocks.Hansen et al. (2015) investigated the lateral structure of this phenomenon with an axisymmetric model and found that while adiabatic shocks generate a net outward radial force in the jet column from high postshock pressures, the radial forces are inward for strongly-cooled isothermal shocks because the increased postshock magnetic pressure dominates the thermal pressure in that case.This inward motion produces a two-component structure: a disk and a spine.The disk spans the width of the interaction region in a manner similar to the cold slab in the hydrodynamic case (as discussed in Markwick et al. 2022), while the spine is the result of magnetic forces drawing cooled material towards the axis.The spine then expands along the z-axis and becomes a secondary jet that flows against the converging flows of the colliding columns.
Figure 1 contains a diagram highlighting the various features of hydrodynamic and MHD collisions.
Efforts have been made to study the MHD of plasma flow collisions in a laboratory setting, with differences in physical scale accounted for by the use of dimensionless parameters (Ryutov et al. 2000(Ryutov et al. , 2001;;Falize et al. 2011).Jets have been produced in the laboratory using both pulsedpower (e.g.Lebedev et al. 2005;Ciardi et al. 2009;Gourdain et al. 2010) and laser-driven (e.g.Albertazzi et al. 2014;Lu et al. 2019;Gao et al. 2019) setups.Of particular relevence to this paper, Suzuki-Vidal et al. (2015) studied pairs of supersonic jets in collision.Among the results of this experiment was the emergence of small-scale structures within the interaction region.
In Markwick et al. (2021), we began a series of simulations to study colliding radiative flows like those of Suzuki-Vidal et al. (2015).The goal of this work has been to both recover the behavior seen in the experiments and understand the more general plasma dynamics of converging, magnetized radiative flows.Our work began with a simplified model, featuring hydrodynamic flows and an analytic form of radiative cooling.These simplifications allowed us to focus our attention on instabilities in the interaction region (specifically the cold slab).In this way they provided a starting point for explaining the origin of small scale structures seen in the experiments.Long-term evolution of the simulations was found to be dominated by bending modes characteristic of the nonlinear thin shell instability (Vishniac 1994), which could be triggered either by sufficiently short cooling lengths or by oscillations resulting from the radiative shock instability (Langer et al. 1981).Meanwhile, we found no evidence for the Field (1965) instability in these simulations.
In another paper (Markwick et al. 2022), we continued our work by focusing on the larger-scale effects resulting from variation of flow parameters, namely density, velocity, and jet radius.While our simulations remained in the hydrodynamic limit, we did begin to use a detailed cooling function Λ(n, T ) which was computed to mimic the cooling behaviour of Aluminium (the material used in the laboratory experiments).When densities or velocities of the two flows are not identical, a net motion of the interaction region can result.Meanwhile a difference in cross-sectional area (i.e.jet radius) results in a deflection of the angle at which shocked lateral outflow emerges from the interaction region.
In this paper, we continue our study of colliding radiative flows by examining the effects of magnetic fields.We limit initial magnetic fields to toroidal geometry.In order to better focus on instabilities, we return to the analytic cooling model used in Markwick et al. (2021).Cooling parameters Λ0 and α are varied, as is the strength of the magnetic fields.Variation of both of these cooling parameters allows us to cover cases of both shorter and longer cooling lengths as well as the presence and absence of the radiative shock instability.
This paper is organized as follows: In Section 2 we discuss the model system and simulation parameters.In section 3, we present the results of the simulations.Section 4 will include a discussion of spine growth and instabilities of the interaction region.In sections 5 and 6, we relate our findings to the laboratory and astrophysical settings, respectively.

METHODS AND MODEL
The simulations in this study were conducted using As-troBEAR1 (Cunningham et al. 2009;Carroll-Nellenback et al. 2013), which is a massively parallelized adaptive mesh refinement (AMR) code that includes a variety of multiphysics solvers, such as radiative transport, self-gravity, heat conduction, and ionization dynamics.Our study uses the magnetohydrodynamic solvers with an energy source term associated with the radiative cooling.Simulations were conducted in three dimensions in cartesian coordinates.Our governing equations are where ρ is the mass density, n is the number density of nuclei, v is the fluid velocity, p = p th + 1 8π B2 is the combined thermal (p th ) and magnetic pressure, I is the identity matrix, and E = 1 γ−1 p th + 1 2 ρv 2 + 1 8π B 2 is the combined internal and bulk kinetic energies.In all runs an average particle mass of 1 amu was used.n 2 Λ(n, T ) is the cooling funciton, which gives the rate of radiative energy loss per unit time per unit volume.It may be worth noting that for this work we ignore magnetic resistivity; this will be considered in a future work.
As in Markwick et al. (2021), we use a power-law cooling function 2 of the form Λ(T ) = Λ0 T T 0 α .Cooling is applied only at temperatures above the floor temperature of the simulation in order to safeguard against runaway cooling.Since realistic cooling curves tend to vanish at low temperatures, this behaviour is justifiable on physical grounds.Reference temperature T0 was fixed at 2.25 × 10 4 K (1.94 eV) across all runs.The simulations were conducted in a space of 64 by 64 by 64 computational units, with one computational unit corresponding to 0.1 mm.Within a radius of 16 computational units (centred on the jet axis), up to three levels of refinement were permitted for a maximum resolution of 0.0125 mm; outside this radius, refinement was limited to two levels for a maximum resolution of 0.025 mm.Extrapolated boundary conditions (specifically, Neumann boundary conditions with a derivative of zero) were used in all directions.
Our simulations feature collisions of two cylindrical jets driven from the top and bottom z-boundaries with speed vi = 70 km s −1 .Jets were initialised as cylindrical regions with height 0.5 computational units, in which the density and velocity were held constant.In an effort to match the experiments of Suzuki-Vidal et al. (2015), the jets were of density 3×10 18 particles per cm 3 , while the ambient medium was set to a density of 5 × 10 17 particles per cm 3 and a temperature of 4320 K.
Initial conditions are such that hydromagnetic pressure equilibirum is established across the jet.To achieve this, we used a profile previously used by Hansen et al. (2015) and others (e.g.Lind et al. 1989;Frank et al. 1998).The initial toroidal magnetic field and pressure inside the jet are given by Here Rm = 0.6 mm is the radius at which the magnetic field is maximized at B = Bm, Rj = 1.5 mm is the jet radius, and βm = 8πp(Rm) is the plasma-beta.Bm is a parameter defining the maximum magnetic field, computed such that p(Rj) is equal to the thermal pressure of the ambient medium.In the hydrodynamic limit, this gives a jet temperature of 710 K, though jet temperatures can be as low as 596 K in the βm = 0.84 case.Figure 2 shows these initial conditions.
Our study varies three parameters: one which describe magnetic fields and two which describe cooling.All parameters are summarised in table 1.First, we vary the plasmabeta between 0.84 and 9.84.We also include two cases of βm = 10 5 , which corresponds to the hydrodynamic limit.The magnetic fields of the two jets are aligned in the same direction; this is done to correspond with both Suzuki-Vidal et al. (2015) as well as protostellor jets, as well as to limit the need to consider the effects of reconnection.
We selected our cooling parameters in a way which corresponds to different hydrodynamic instabilities as studied in Markwick et al. (2021).First, cooling power α is varied between α = −1 and α = +1.Based on results from Strick-  m /(8π) .We also note the expected sonic M = mach numbers, computed post-collision using γ = 5 3 .
land & Blondin (1995) and Markwick et al. (2021), the former case is likely to exhibit the radiative shock instability (Langer et al. 1981), while the latter is unlikely to do so.In addition to the power law exponent, we also vary cooling strength Λ0 between low and high strengths.In Markwick et al. (2021), we found that momentum conservation for a γ = 5 3 gas gives a cooling length of approximately . We adjust the values of Λ0 with α such that the cooling length remains fixed: Our 'low' strengths are defined as Λ0 = 1.0 × 10 −23 erg cm 3 s −1 for α = −1 and Λ0 = 2.0 × 10 −23 Low βm, low Λ 0 , α = +1 Figure 3. Run 1 from table 1, which features a relatively strong magnetic field, along with weak cooling which weakens as temperature decreases.A prominent feature is the formation of a spine, driven by magnetic hoop stresses (∇B 2 ϕ ), which eventually collapses as a result of an instability.In this figure and all figures within seciton 3, images are taken as midplane slices of density, and tick marks are spaced at 1mm intervals.erg cm3 s −1 for α = +1, while the 'high' strengths are equal to ten times the corresponding low strength.

Weak Cooling
We begin our results with the four cases in which Λ0 is of order 10 −23 erg cm 3 s −1 .In all four of these cases we observe prominent formation of the spine, though specific details vary with both the magnetic field strength (characterized by βm) and cooling behaviour (determined by exponent α) are varied.
We begin with the case of strong magnetic field (βm = 0.84) and a positive cooling exponent (α = +1.0).This run  1 or figure 3).After an initial surge, the growth speed remains roguhly constant until around 500 ns.After this time, the spine growth decelerates at an approximately constant rate.
can be seen in figure 3 which shows density slices through the 3-D simulations.At the edges of the jet we see material ejected laterally into the surrounding medium.Such shocked lateral outflows (SLO) where seen in the hydrodynamic simulations Markwick et al. (2021Markwick et al. ( , 2022)).These outflows are driven by the difference between high pressure within the interaction region and the lower pressure of the ambient medium.The initial lateral pressure discontinuity smooths into a gradient of thermal pressure spanning the outer portion (r > Rm) of the interaction region 3 .
Almost immediately after the collision between the two counter-flowing jets, magnetic forces in regions of the cold slab closer to the axis become stronger than thermal pressure gradients in those regions.The dominance of magnetic forces in the inner region causes material to be drawn towards the axis, increasing the density (Hansen et al. 2015;De Colle et al. 2008).Soon afterwards a spine -an axial bipolar flow propagating back into both jets -begins to emerge along the axis.The spine quickly grows to a length comparable to the diameter of the counter propagating jets.As the spine propagates, it pushes through jet material, deforming it.This back-reaction modifies the structure of the interaction region and the two "working surfaces" which bound it.After an initial surge and up through around 500 ns, the growth speed of the spine is approximately constant, as seen in figure 4.
At later times, the spine becomes subject to what appears to be a form of a kink instability.Once this occurs, the disrupted spine is unable to support itself against the ram pressure of the jets and begins to decelerate back towards the centre of the interaction region.This can be seen in figure 4, in which the growth speed of the spine undergoes approximately constant deceleration until it eventually begins to collapse inwards.In the process of collapse, the entire interaction region is strongly distorted.
The creation of the spine leads to a significant change in the evolution of the interaction region compared with the hydrodynamic cases (Markwick et al. 2021), in which instabilities of the cooling region and the disk dominate longterm evolution of the interaction system.
Next we examine Run 2 which is similar to the Run 1, but uses a different power law exponent.Run 1 had α = +1.0,indicating that cooling is stronger at high temperatures and weakens as the material cools.Run 2 uses α = −1.0,indicating that cooling is weaker at high temperatures but strengthens as the material cools.The initial postshock cooling strength is adjusted so that both runs feature similar cooling lengths, as can be seen by comparing figures 3 and 5.
Unlike in the α = +1.0case, the cold slab in the α = −1.0retains a noticeable amount of cooled material which we refer to as the disk.This is likely the result of the nature of density gradients within the cooling region: for the α = +1.0case, the gradient is appreciable throughout the cooling region; for cases with α < 0, since the cooling strengthens at lower temperatures there is a sharper transition to the cold slab (see figure 3 of Markwick et al. (2021)).
After some time, the cold slab experiences significant perturbations.Oscillations drive cold material into the cooling regions above and below.As these perturbations grow Low (left) vs high (right) βm.Low Λ 0 , α = +1 Figure 6.Runs 1 (left) and 3 (right) from table 1.Both cases feature identical cooling (weak cooling which weakens as temperature decreases), but run 3 has a weaker magnetic field.As a result, the spine forms more slowly and doesn't grow as large.
they lead to interactions with the spine, which excites further instabilities.We note that Run 2's perturbed cold slab leads to the kink occurring in the spine long before it was seen in Run 1.
Our next case, Run 3, is shown on the right side of figure 6. Run 3 features identical cooling to Run 1 (low Λ0, α = +1.0),but has a weaker magnetic field (βm = 9.84).
The most significant differences in the results of Runs 1 and 3 (which in setup differ only in the value of βm) are found in the spine.To understand these differences, however, it is first helpful to focus on the evolution of the disk.
In the βm = 9.84 case the disk more prominent than in its βm = 0.84 counterpart, though it is not as prominent as it is in the α = −1.0(even-numbered runs) or high Λ0 (Runs 5-10) cases.This is likely a result of magneticallydriven collimation forces being weaker which allows more of the cooled material to remain in the disk rather than being drawn into the spine.Since the disks is partly responsible "feeding" the spine, such differences should be expected to effect the latter's evolution.These expectations are born out in Run 3 as the initial formation of the spine takes nearly twice as long to begin as in Run 1.After spine formation begins subsequent growth is also significantly slower in the βm = 9.84 case (Run 3) than its βm = 0.84 counterpart (Run 1).This is likely the result of weaker magnetic fields driving the collimation process, especially in cases such as these where a limited amount of cold slab material is available for collimation.As with the βm = 0.84 case, the spine is eventually disrupted by an effect which appears to be similar to the kink instability.This disruption begins at a similar duration after the initial spine growth begins, which results in a shorter spine length at the time of collapse.
We now examine a case with α = −1.0 and βm = 9.84 (figure 7), which combines some of the features of Run 2 and Run 3. As with the other α = −1.0case (Run 1), the disk is more prominent when the cooling exponent is negative.
Also like the other α = −1.0run, the strong cooling at lower temperature drives strong NTSI in the disk which then interacts with the spine.Once again we see these instabilities leading to a collapse of the spine before the spine instability has time to fully develop.
As with the other weaker field βm = 9.84 case (Run 2), spine growth is somewhat delayed compared to the stronger field βm = 0.84 cases (see figure 15 in section 4).Here however the delay is far less than Run 3, the α = +1.0,βm = 9.84 case.The more rapid onset of spine formation is likely the result of the stronger cooling at low temperature.More disk material is available in this run to feed spine growth.Once the spine begins to form however its propagation speed is the lowest of any of our βm ⩽ 10 cases (see table 2 in section 4), though it is still comparable to the other βm = 9.84 cases.

Strong Cooling
For our remaining simulations, we now increase the cooling strength Λ0 by a factor of 10, resulting in cooling of order 10 −22 erg cm 3 s −1 .With stronger cooling, the radiative shock and nonlinear thin shell instabilities of the cooling region and disk contribute significantly to the long-term evolution of the system, in some cases before the spine is able to properly form.
We begin examining the effects of this change by comparing low and high cooling strengths for the models with βm = 0.84 and α = +1.0,(i.e. for strong field, positive cooling exponent).Results are shown in figure 8.In the case of stronger overall cooling due to higher Λ0, the spine emerges as earlier compared with the low Λ0 runs.We note that the spine has a somewhat larger radius.A prominent disk is also present.Both of these effects can be traced to the shorter cooling length expected for higher Λ0: a smaller cooling region limits opportunity for cooling matter to escape via lateral outflows, so a greater portion of the initial flow (extending into r > Rm) is destined for the cold slab.
When we compare the length of the spines, we find that over a sufficient length of time cooling strength Λ0 does not appear to have an appreciable effect.Like the α = +1.0cases with lower cooling, the instability of the growing spine still appears to be the dominant effect in its evolution (as opposed to any instabilities which may arise within the disk).
In figure 9, we again consider strong cooling (high Λ0) with βm = 0.84, but we we add the effects of α = −1.0.
Here we begin to see some of the hydrodynamic instabilities studied in Markwick et al. (2021) play a role in shaping the evolution interaction region, though magnetic fields continue to have effects as well.
Not long after collision occurs, the radiative shock instability promotes oscillation in size of the cooling region.The timescale of these oscillations is inversely proportional to the cooling length, making them more prominent in the cases with stronger cooling.As shown in the first five panels of figure 9 and discussed in Markwick et al. (2021), these oscillations imprint variation within the disk, allowing instabilities which disrupt the entire interaction region to occur earlier than in the corresponding low Λ0 case.
We now consider pairing strong cooling (high Λ0) with weaker fields (βm = 9.84).If α = +1.0(figure 10), the result is an interaction region which exhibits similarities to several  1.Both feature relatively strong magnetic fields and positive α, but the absolute strength of cooling is higher in run 5.While the spines are of comparable length, stronger cooling provides more cold material, resulting in a wider spine.
other cases before being disrupted by multiple simultaneous instabilities.
The disk in this case is similar to that of the corresponding βm = 0.84 case with the same cooling parameters.We see the disk is thickened by the stronger total cooling while lacking the disturbances of the radiative shock instability.
Low βm, high Λ 0 , α = −1 Figure 9. Runs 6 from table 1, which features strong fields and strong cooling which strengthens as temperature decreases.The first five images highlight the oscillations of the cooling region at an earlier time window.At later times we zoom out to show spine growth; stronger cooling provides more cold material, resulting in a wider spine.
The spine is also of a similar thickness to the other strong cooling cases.Finally, the spine growth rate is comparable to other βm = 9.84 cases, but the onset of its growth is closer to the α = −1.0,low Λ0 case than it is to the α = +1.0,low Λ0 case.
As expected, with stronger overall cooling, bending modes characteristic of the nonlinear thin shell instability appear within the disk before the spine begins to develop instabilities.As these bending modes grow, the NTSI alters High βm, high Λ 0 , α = +1 Figure 10.Run 7 from table 1, which features weaker magnetic fields and strong cooling which weakens as temperature decreases.Stronger cooling in the disk allows the spine to grow larger, though instabilities in the disk also contribute to the spine's eventual collapse.
disk at a similar time to when the spine undergoes its own instability.
If we combine weaker cooling, α = −1.0,and βm = 9.84 (shown in figure 11), the hydrodynamic instabilities occur at an earlier time and prevent the formation of a spine.As was the case in Markwick et al. (2021), the radiative shock instability imprints variation onto the cold slab, allowing the NTSI to begin at an earlier timescale than it would otherwise.At the same time, weaker magnetic fields (βm = 9.84) result in a delay in the formation of the spine (We do note however that the magnetic fields do still cause cold slab material to congregate near the axis).The combination of these effects places spine formation at a later time than the disruption of the disk by the NTSI, therefore preventing the former altogether.
Lastly, we ran a pair of cases with βm = 10 5 , again using strong cooling.This should approximate the hydrodynamic limit we studied in Markwick et al. (2021).As seen in figure 12, we do recover the results from that paper.Without a significant magnetic field, a spine does not form.Instead the cold slab forms a disk which, over time, becomes unstable.It is the NTSI which drives bending modes which eventually disrupt the interaction regions.As expected, this occurs sooner in the case where the cooling law exponent is negative, since the radiative shock instability imprints a seed for the NTSI.  1, which features which features weaker magnetic fields and strong cooling which strenghtens as temperature decreases.The combination of radiative shock and nonlinear thin shell instability serve as the primary drivers of this run's evolution, the presence of a nontrivial magnetic field still results in cold slab material moving towards the centre.

Magnetic fields in the Cold Slab
We now present a simplified model for the formation and propagation of the spines seen in our simulations.We do this by considering unperturbed conditions for the spine formation.We must first find expressions for the density ρ d , pressure p d , and magnetic field B d of the disk in terms of the corresponding preshock values ρi, pi, Bi, and jet velocity vi.Also relevant to this calculation is sound speed cs = γp i ρ i .Radial variation of the preshock magnetic fields will result in some of these quantities being dependent on radius.
Approximating the cooling region as a stationary isothermal shock (i.e. the quantity p ρ returns to its preshock value when it reaches the cold slab), conservation of mass and momentum can be expressed as where βr is value of 8πp i B 2 i at a radius r.In the βr >> γM 2 limit, we recover the ρ d = γM 2 ρi result found in Markwick We note that M varies as p and βm varies as pi, so equation 7 varies inversely with B and is independent of pi.Using equation 2 for our initial B(r), the density and pressure jumps are therefore given as where pi(r) is given by equation 3. Flux freezing tells us that the magnetic fields grow proportional to density; the magnetic fields in the disk are approximately given as B d = M √ γβmBm, which is independent of position.
Figure 13 shows density, pressure, and magnetic fields in the disk, both using the βm << γM 2 approximation (equation 7) and the full isothermal approximation (equation 6).We find that the two approximations are reasonably close to each other, especially for stronger magnetic fields (lower βm).The βm << γM 2 approximation however is inaccurate at locations very close to the axis, since the magnetic fields are weak and we instead approach the hydrodynamic jump conditions.

Growth of the spine
We now develop a simple model for the growth of the spine in which magnetic forces in the disk draw material towards the axis.The build-up of magnetic and gas pressure there then create a collimated bipolar flow pushing back into both jets.
We begin by assuming the spine has uniform density ρs and grows with velocity vs. Let 2hs and rs be the height4 and radius of the spine.The growth of the spine is driven by the high magnetic pressure of the cold slab, so we define the total pressure of the disk to be p d + B 2 d 8π .It is also fed momentum flux −ρiv 2 i from above.We have assumed the ram pressure from the disk itself to be negligible, as is the preshock pressure5 and magnetic field6 .
Consider a cylindrical region of radius rs and height hr centered on the spine's axis.Let one end of the cylinder be located in the preshock region, while the other end is at the disk.Such a region is pictured in figure 14.The momentum of the spine within this region is given by (ρsvs)(πr 2 s hs), while the momentum in the rest of the region is given by (−ρivi)(πr 2 s [hr − hs]).Here we denote vi as being positive towards the disk, and vs being positive away from the disk.Since ρi, vi, rs, and hr are all constant, vertical momentum conservation is given by Assuming density and velocity to be constant, we have Consider now a cylinder similar to the one we have been examining, but with the lower boundary located at some Away from the jet axis, the two approximations agree reasonably well, especially for lower βm.distance h− above the disk (see figure 14).Mass conservation for this new cylinder is given by where ρc is the density of the cooling region immediately ahead of the spine head7 .We can solve this equation for ∂hs ∂t : Plugging this into equation 11 gives us Assuming ρc << ρs gives us Since ρi ⩽ ρc we must also have ρi << ρs, and if we assume8 ρsvs ∼ ρivi we can also infer vs << vi.We therefore can eliminate the factor in square brackets as being close to unity: Let us strengthen our assumption to ρsvs = ρivi exactly.Dividing by ρs gives us Solving the quadratic gives which can be approximated as Given that the cooling region is approximately an isothermal shock, p d ρs is equal to p i ρ i , which is in turn To replace the remaining terms B d and ρs with their Table 2.The observed speeds for the growth of spine size hs in the simulations, compared to the estimated speeds for the corresponding value of βm.Times are given in ns, ∆hs is given in mm, and speeds are given in km s −1 .
preshock counterparts, we consider how alfven speed vA = B √ 4πρ scales across a shock.The density increases by a factor of M √ γβm, but so does the magnetic field as a result of flux freezing.The alfven speed therefore increases by a net factor of M 1 2 (γβm) 1 4 , leaving our overall equation as , this is equal to Finally, going back to equation 13 tells us that the growth speed is approximately twice this, or . In terms of βm, this can be expressed as To compare our estimate to our simulations, we selected a starting and stopping time for each simulation and measured the change in hs over that time interval.The results are shown in table 2, and a time plot of spine growth is shown in figure 15.
We find that the model predicts growth speeds which compare well with the simulations given the model's simplicity, with differences ranging between approximately 10% to 50% For runs 1 and 5 (which feature βm = 0.94 with α = −1), we find that the model predicts growth speeds which overestimate the speeds by only about 12%.Runs 2 and 6, which are the corresponding α = −1 cases, are slower than their α = +1 counterparts despite this not being accounted for in the estimation, while the βm = 9.84 cases (runs 3, 4, 7, and 8) show faster observed growth than is estimated.In general, runs with positive cooling exponent and lower beta (stronger field) produce the best correspondence to the model.

Timescale Analysis of Disk Instabilities
In addition to unstable behaviour which causes the spine to break apart, there are two additional hydrodynamic instabilities occurring within the interaction region: the NTSI in the disk, and the radiative shock instability in the cooling region.The latter affects the system by imprinting variation onto the disk and promoting the NTSI.
Before examining these two instabilities further, we briefly comment on the breakup of the spine.The nature of the dynamics in this case, in terms of a dominant mode, is unclear as the β value in the spine is intermediate between hydrodynamic (in which instabilities such as Kelvin Helmholtz might apply) and magnetically dominated (in which instabilities such as kink or sausage might apply) regimes.We observe timescales for the spine to break-up on the order of a few hundred nanoseconds after the onset of its growth.We leave a more detailed analysis of its dynamics for a future project and focus on the the jet interaction region.

Nonlinear Thin Shell Instability
The nonlinear thin shell instability promotes growth of bending modes within the disk.We define the amplitude of a perturbation as ψ, and its wavelength as λ.Its growth speed in the hydrodynamic case given as ψ ∼ ψ λ 3 2 cs.
(25) (Vishniac 1994).An appropriate timescale for growth would therefore be There is a constraint that ψ and λ must both exceed the thickness of the system (i.e. of order d cool ).With this in mind, we approximate our initial timescale as We therefore find that for our strong cooling cases we have an estimated timescale of order 170 ns; this appears to be in agreement with the evolution of the α = +1.0hydrodynamic case presented in figure 12.
We must consider the impact of magnetic fields on the NTSI.Heitsch et al. (2007) found that magnetic fields perpendicular to the shock can inhibit the NTSI, but that the instability will often still develop.For toroidal magnetic fields, we expect magnetic tension to only resists perturbations with azimuthal variation; radial variation would be perpendicular to the magnetic field, yielding magnetic tension equal to the unperturbed case.
Assuming a purely azimuthal perturbation of the form ψ ∼ e imϕ and a sufficient distance from the axis, the Euler equation can be expressed as where ψhydro is the growth in velocity the hydrodynamic case, which per equation 25 is governed by . Meanwhile flux freezing should bend the field into a sinusoidal shape to match our perturbation, giving us Bz of order Since β ′ m in the disk is smaller than global βm by a factor of order M √ γβm, we have .
Annular Structure for the Nonlinear Thin Shell Instability Figure 16.An isosurface plot of density in the α = −1, strong cooling, weak field case.Note that perturbations with azimuthal variation are suppressed, with the most prominent perturbations featuring annular structure.
The above calculation suggests that magnetic tension will prevent growth of toroidal perturbations smaller than ψ ∼ 16π 2 3 γM 2 βm 1 2 λ, which is quite large compared to the wavelength.However, as noted earlier, radial perturbations (such as those seeded by the radiative shock instability) should be unaffected by magnetic tension.Heitsch et al. (2007) also found that other conditions such as sufficient turbulent energy may also allow magnetic tension to be overcome.Examining our simulations (e.g.figure 16) we find that perturbations do in fact adopt an annular structure; the little azimuthal variation that is present appears to be most prominent along the grid axes, suggesting numerical artefact.

Radiative Shock Instability
The radiative shock instability produces oscillations in the size of the cooling region; these oscillations can disrupt the rest of the interaction region by imprinting spatial variation on the cold slab, which excites the NTSI.An appropriate time scale for the radiative shock instability is that of its oscillation period, which is given by (Chevalier & Imamura 1982): where f (α) = π 2δ I for δI given in table 1 of Chevalier & Imamura (1982).For α = −1, f (α) = 5.97 for the fundamental mode and f (α) = 1.65 for the first overtone.
In figure 9 shows this case undergoing an oscillation period over the course of 40 seconds or so.In Markwick et al. (2021) we found the cooling time to be approximately where ns is the immediate postshock density.Using this estimate, we find that the first overtone has an estimated period of 57 ns.
Comparing timescales for the radiative shock and nonlinear thin shell instabilities offers insight into the observed behaviour of the disk.At early times, the radiative shock instability -which for high mach number has a shorter initial timescale than the NTSI -causes growth to occur sooner than it does as a result of the NTSI alone, as seen in figure 12 and discussed in Markwick et al. (2021).As perturbations grow, the timescale for the NTSI is reduced, so this instability has a tendency to dominate at later times if present.

COMPARISON TO EXPERIMENTAL RESULTS
One of the primary goals of our simulations is to better understand the experiments performed by Suzuki-Vidal et al. (2015).In those experiments, a region of compressed material is observed behind the leading bow shock.As seen in figure 3 of that paper, the height of the compressed region is more than twice width of the jets, while none of our previous hydrodynamic simulations (Markwick et al. 2021(Markwick et al. , 2022) ) showed an interaction region which was taller than the width of the jets.In our present work however we see that that magnetic fields allow for the formation of a spine, which can grow to such a height.We also see that the NTSI and radiative shock instabilities continue to be present in magnetized jets; together with the instability of the spine, these instabilities provide a mechanism for the observed fragmentation of the interaction region.Best exhibiting the combination of these effects are cases 6 (figure 9) and 7 (figure 10); both of these feature cooling which is strong enough for hydrodynamic instabilities to contribute, but not so strong (relative to the magnetic field) that the instabilities disrupt the system before a spine is able to fully form.While our latest set of simulations has brought us closer to the experiments of Suzuki-Vidal et al. (2015), it is worth noting that key differences remain.First, we note that the simulations presented in this paper used an analytic cooling model, while cooling in the experiments has a more complicated dependence on temperature and density.A more accurate model was discussed in Suzuki-Vidal et al. (2015) and was used in Markwick et al. (2022).Second, the experiments do not show a spine in the same way that our simulations do: a spine-like structure is observed on one side of the interaction region, but unlike our simulations it only is observed one one side of the cold slab.This may be a result of asymmetry: the experiments feature jets of differing speeds and densities.We previously studied the effects of collisions between non-identical jets in Markwick et al. (2022), but those simulations did not include magnetic fields and thus lacked a spine.In a future work we hope to combine the use of a more sophisticated cooling function and non-identical jets studied with the presence of magnetic fields.

CONCLUSION
In this paper we have extended the results of Markwick et al. (2021) by examining the structure and evolution of colliding flows with toroidal magnetic fields.The hoop stress of the toroidal fields provides a collimating force which results in the emergence of a dense spine along the axis.The growth speed of the spine is determined primarily by the strength of the magnetic fields.After some time, one or more instabilities will contribute to the disruption of the interaction region.For flows with sufficient magnetic fields for spine formation, the spine eventually becomes unstable and collapses after growing for a sufficient amount of time.For flows dominated by radiative cooling, the hydrodynamic instabilities studied in Markwick et al. (2021) will disrupt the cold slab, which in turn disrupts the spine at an earlier time than would occur as a result of the spine instability alone.
The results of colliding flow experiments (both those performed in the laboratory as well as those performed via simulation) can be used to gain insight into astronomical phenomena such as YSO jets.While astrophysical jets do not ordinary collide in the same manner as in experiments, the formation of shocks produced by flow collisions can be connected (via reference frame transformation) to those formed by jet pulsation in which a slower region is overtaken by a faster region of the jet behind it (e.g.Gardiner et al. 2000).Spine growth can be observed in simulations such as Hansen et al. (2015), but in that case the axisymmetric nature of the simulation may have inhibited instabilities; our work here offers some insight into the lateral structure of the spine and its eventual breakup.Conversely, the 3D simulations of Hansen et al. (2017) exhibit the nonlinear thin shell instability, but do not exhibit spine growth or instability as they lack magnetic fields.Real astrophysical jets are produced by magnetic effects (e.g.Blandford & Payne 1982;Lynden-Bell 1996), and so the interaction of magnetic fields with instabilities must be considered.

ACKNOWLEDGMENTS
This work used the computational and visualization resources in the Center for Integrated Research Computing (CIRC) at the University of Rochester.Financial support for this project was provided by the Department of Energy grants GR523126, DE-SC0001063, and DE-SC0020432, the National Science Foundation grant GR506177, and the Space Telescope Science Institute grant GR528562.E.G. Blackman acknowledges working in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452.

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.

Figure 1 .
Figure 1.A diagram showing the structure of colliding MHD columns (jets).The nature of structures such as the disk and spine are articulated in the text.For comparison with a purely radiative hydrodynamic case see Markwick et al. (2021).

Figure 2 .
Figure 2. Thermal Pressure (red) and magnetic fields (Blue) in the initial (pre-collision) jets

Figure 4 .
Figure 4.The size (red) and growth speed (blue) of the spine in Run 1 (see table1or figure3).After an initial surge, the growth speed remains roguhly constant until around 500 ns.After this time, the spine growth decelerates at an approximately constant rate.

Figure 5 .
Figure5.Run 2 from table 1, which features a relatively strong magnetic field along with weak cooling.Unlike the previous case, cooling strengthens as temperature decreases.Additional instabilities result in an earlier collapse of the spine.

Figure 7 .
Figure7.Run 4 from table 1 (weak magnetic field with weak cooling which strenghtens as temperature decreases).Once again, instability of the cold slab contributes to an earlier collapse of the spine.

Figure 8 .
Figure 8. Runs 1 (left) and 5 (right) from table1.Both feature relatively strong magnetic fields and positive α, but the absolute strength of cooling is higher in run 5.While the spines are of comparable length, stronger cooling provides more cold material, resulting in a wider spine.

Figure 11 .
Figure11.Run 8 from table 1, which features which features weaker magnetic fields and strong cooling which strenghtens as temperature decreases.The combination of radiative shock and nonlinear thin shell instability serve as the primary drivers of this run's evolution, the presence of a nontrivial magnetic field still results in cold slab material moving towards the centre.

Hydrodynamic limit with high Λ 0 Figure 12 .
Figure 12.Runs 9 (left) and 10 (right)from table 1.These runs reproduce the NTSI-driven hydrodynamic case seen in Markwick et al. (2021), with the radiative shock instability enhancing the NTSI in run 10

Figure 13 .
Figure 13.Estimated Density, Pressure, and Magnetic Field in the disk.The solid lines are plotted using equation 6, while the dashed lines show the β << γM 2 approximation (equation 7).Away from the jet axis, the two approximations agree reasonably well, especially for lower βm.

Figure 14 .
Figure 14.A diagram showing the regions for which mass and momentum conservation are to be considered for the spine.

Figure 15 .
Figure 15.A) plot of the size of the spine vs time.B) a time derivative of the spine size in the α = +1.0cases, averaged over intervals of 50ns.In both plots, the dotted black lines match the estimated growth speeds for βm = 0.84 (upper line) and βm = 9.84 (lower line).