Magnetic Fields in Multiphase Turbulence: Impacts on Dynamics and Structure

Both multiphase gas and magnetic fields are ubiquitous in astrophysics. However, the influence of magnetic fields on mixing of the different phases is still largely unexplored. In this study, we use both turbulent radiative mixing layer (TRML) and turbulent box simulations to examine the effects of magnetic fields on cold gas growth rates, survival, and the morphology of the multiphase gas. Our findings indicate that, in general, magnetic fields suppress mixing in TRMLs while turbulent box simulations show comparatively marginal differences in growth rates and survival of the cold gas. We reconcile these two seemingly contrasting results by demonstrating that similar turbulent properties result in comparable mixing -- regardless of the presence or absence of magnetic fields. We, furthermore, find the cold gas clump size distribution to be independent of the magnetic fields but the clumps are more filamentary in the MHD case. Synthetic MgII absorption lines support this picture being marginally different with and without magnetic fields; both cases aligning well with observations. We also examine the magnetic field strength and structure in turbulent boxes. We generally observe a higher mean magnetic field in the cold gas phase due to flux freezing and reveal fractal-like magnetic field lines in a turbulent environment.


INTRODUCTION
The fact that most of the astrophysical media are multiphase in nature is a well-established one, from observational (e.g.Tumlinson et al. 2017;Veilleux et al. 2020), numerical, and theoretical investigations (McKee & Ostriker 1977;Donahue & Voit 2022;Faucher-Giguere & Oh 2023).The multiphase nature of the interstellar (ISM), circumgalactic (CGM) and intracluster (ICM) medium is also expected to play an important role in the overall evolution of the associated systems, from the general baryon cycle to feedback processes (Veilleux et al. 2005;Péroux & Howk 2020).However, there are many aspects of multiphase media, like survival and characteristic size of cold media, that are still in question and are an active field of research.
Many forays towards understanding the multiphase gas use an idealised version of the medium.There are studies which focus on the development of the multiphase gas by condensation from an initially static hot ambient medium via thermal instability in 1D (Sharma et al. 2012;Waters & Proga 2019), 2D (McCourt et al. 2018) and 3D simulations (Gronke & Oh 2020b).Such studies are a good way to isolate and investigate the role of different factors like magnetic fields (Sharma et al. 2010;Ji et al. 2018), metallicity (Das et al. 2021), gravity (Mccourt et al. 2012), density fluctuations (Choudhury et al. 2019), rotation (Sobacchi & Sormani 2019) or cosmic rays (Butsky et al. 2020).
It is also well known that the astrophysical media are highly tur-★ E-mail: hitesh@mpa-garching.mpg.debulent, due to their high Reynolds number.This has been shown by many observational (Elmegreen & Scalo 2004;Falceta-Gonçalves et al. 2014;Vidal-García et al. 2021;Li et al. 2022) and numerical studies (Brandenburg & Nordlund 2011;Federrath 2013;Burkhart et al. 2020).Hence, many studies like Mohapatra & Sharma (2019); Gronke et al. (2022); Mohapatra et al. (2022b,c) investigate the evolution of the multiphase gas in the presence of a turbulent astrophysical media.Turbulence both amplifies and destroys multiphase gas.The density and temperature perturbations in a turbulent medium can enhance the creation of multiphase gas, while the same turbulent motions can mix the existing multiphase gas, which might further cool or mix away depending on the cooling timescale.Not just astrophysics, multiphase turbulence is also a very relevant topic at more terrestrial scales and is also an active field of research in general fluid dynamics circles, as there are many applications like combustion dynamics, smoke transport and meteorology where multiphase interactions play a crucial role.One seminal result in the field was by Damköhler (1940), where they found that the behaviour of a flame front in a turbulent medium differs depending on the ratio of the reaction and turbulent timescales, i.e. the Damköhler number.Tan et al. (2021) further explore this parallel in the context of hydrodynamic turbulent radiative mixing layers (TRMLs) in astrophysical media.
Generally, there are three stages in the evolution of a turbulent multiphase medium, with many studies examining each.First, the production or presence of seed multiphase gas.The exact mechanism of this can vary from medium to medium.For the CGM, this can either be in the form of multiphase ISM transported into the CGM by feedback mechanism, or via condensations from the hot medium due to thermal instability.The second stage is the growth of one of the phases in the multiphase gas.And, the final stage is the equilibrium or steady state.
In order to understand the second stage of the multiphase gas evolution, Gronke et al. (2022) study the growth of cold gas in a thermally stable, ambient turbulent medium.They found a critical radius for the size of the seed cold gas cloud in a given turbulent hot ambient medium (akin to the survival criterion previously found for laminar flows Gronke & Peng Oh 2018;Li et al. 2020;Kanjilal et al. 2021).Their results also agreed with the expectations from the previous hydrodynamics TRML results from Tan et al. (2021), indicating TRMLs might be the principle mechanism for mixing in a multiphase medium.
But, apart from being multiphase and turbulent, the astrophysical media are also known to possess substantial magnetic fields as seed primordial magnetic fields are amplified due to structure formation and other baryonic dynamics (Dimopoulos & Davis 1997;Subramanian 2015).There are many observational evidences for ubiquitousness of magnetic fields using techniques like Faraday rotation (Dreher et al. 1987;Kim et al. 1990;Taylor & Perley 1993;Clarke et al. 2001), dust alignment (Ade et al. 2015), and others (Lopez-Rodriguez et al. 2021).And, many numerical studies also point to a similar prevalence of magnetic fields (e.g., Pakmor et al. 2020;van de Voort et al. 2021).
The presence of magnetic fields can be disruptive to mixing via TRMLs.It is well-known that linearly, the Kelvin-Helmholtz instability is suppressed for specific magnetic field orientation (Chandrasekhar 1961), while Ji et al. (2019) show that the nonlinear evolution of the instability with radiative cooling is suppressed for all orientations of the initial magnetic field. 1 Hence, the inclusion of magnetic fields may change the overall evolution of multiphase gas, resulting in different survival criteria and cold gas growth rates.In summary, while it has been shown in recent work that mixing and subsequent cooling can lead to the survival and even the production of cold gas, and thus explain the ubiquitous presence of multiphase gas in turbulent systems -where this cold gas should be destroyed rapidly -magnetic fields might ruin this picture by preventing mixing and hence hindering cooling.
In this paper, we investigate the influence of magnetic fields on the general phenomenon of mixing between the phases in a multiphase gas.For that purpose, we use two kinds of simulations, TRMLs and turbulent boxes, with and without magnetic fields.First, we expand on the parameter space for TRMLs explored in previous studies, to confirm the suppression of mixing for different cooling strengths (and hence different Damköhler numbers).Second, we check for the effects of including magnetic fields in turbulent box simulations similar to Gronke et al. (2022).
The structure of the paper is as follows.We explain the numerical setups for both TRMLs and turbulent boxes in section 2. We present the results from the TRML simulations in section 2.1, and the turbulent boxes in section 4.Then, we discuss our results in section 5 and conclude in section 6.The visualisations related to this study can be found here.2multiphase.html

NUMERICAL SETUP
For our simulations, we use the ATHENA++ code (Stone et al. 2020).We use the default HLLC solver for our hydrodynamic (HD) simulations and the default HLLD solver for our magnetohydrodynamic (MHD) simulations, with Piecewise Linear Method (PLM) applied to primitive variables, second-order Runge-Kutta time integrator, adiabatic EOS and a cartesian geometry.Similar to Gronke et al. (2022), we implemented the Townsend radiative cooling algorithm (Townsend 2009) for computing the radiative losses, using a cooling curve at solar metallicity fitted using 40 segments of power laws.We also enforce a temperature floor  floor = 4 × 10 4 K.

Turbulent Radiative Mixing Layer (TRML)
Turbulent, radiative mixing layers in the astrophysical context have been investigated in the past (e.g., Begelman & Fabian 1990;Slavin et al. 1993;Kwak & Shelton 2010;Hillier et al. 2023;Fielding et al. 2020).We use the same numerical setup as the one used in Tan et al. (2021) and Ji et al. (2019) to investigate the Turbulent Radiative Mixing Layer (TRML), with a small difference in our coordinates convention.The shear velocity in our simulations is along x and the cold/hot interface is normal to ẑ.The different density and velocity profiles are, () =  shear 2 1 + tanh (2) where   and   are set to 2/ box,x and 2/ box,y respectively.We use  shear = 100km/s, 50km/s as the shear velocity (corresponding to sonic Mach numbers of M ≡  shear / s,hot ∼ 0.33, 0.16),  = 0.01 shear ,  =  box,z /20 for the interface thickness, and M A = 1, 10 as initial Alfvénic Mach number.Furthermore,  hot = 1.6×10 −4 cm −3 is the density of hot medium,  cold = 1.6×10 −2 cm −3 is the density of cold medium, and B0 is the initial magnetic field direction.We use a floor temperature of 4 × 10 4 K, and we stop cooling at temperatures above 0.5 hot to emulate the effect of heating, where  hot is the hot medium temperature.We initialise the cold medium at the floor temperature (4 × 10 4 K), and use a fixed pressure over the whole box to ensure pressure balance in the initial conditions.This corresponds to a hot medium at  hot =  cold  with  ≡  cold / hot = 100 being the overdensity.We impose an outflow boundary condition along the normal to the cold/hot interface ( ẑ), and periodic boundary conditions in all other directions.Note that due to the self-similarity of the solution, the chosen numerical values for  and  shear are unimportant (as long as the critical dimensionless quantities are kept constant).
We use a resolution of 64 × 64 × 640 in x, ŷ and ẑ directions respectively, and use different box sizes to vary parameters in our simulations, but keep the ratios of box lengths in different dimension fixed,  box,z = 10 box,x = 10 box,y .We vary the Damkohler number (Da =  turb / cool =  box,x or y /( turb  cool )) in a range of ∼ 10 −4 − 10 4 by changing  box .
In such TRML simulations, the mixing layer tends to move into the hot medium as it is consumed and more cold gas is created.This can cause the mixing layer to go out of the computational domain, especially for high Da cases.To counter that, we add a velocity to the whole box in the opposite direction to keep the mixing layer inside the computational domain.This velocity is calculated using the difference between the current cold/hot interface position and the original cold/hot interface position ( = 0).We verify that this does not affect the mixing rates, and only increases the time the mixing layer spends inside the computational box.We denote the different TRML simulations as Ma(A)_Bx(B), where A is the Alfvénic Mach number and B is the initial magnetic field orientation.

Driven turbulence boxes
We use a separate simulation setup, similar to the one used in Gronke et al. (2022), to study the behaviour of cold gas in fully-developed turbulence with (MHD) and without (HD) magnetic fields.We start with a box filled with isobaric gas at uniform density and temperature ( hot = 4 × 10 6 K), with solar metallicity and H-abundance.In our MHD simulations, we initialize the box with a uniform magnetic field.We use the Ornstein-Uhlenbeck (OU) process (Eswaran & Pope 1988;Schmidt et al. 2006) to drive the turbulence at the largest scale ( = 2/ box ), i.e. the box size.We use driving timescale of 0.001  eddy , correlation timescale ∼  eddy and solenoidal to compressive fraction  sol = 0.3.We also maintain a  box / cloud = 40 for all our simulations.
We drive the turbulence for 7  eddy with the cooling turned off, which gives the setup enough time to reach a steady-state with equilibrium kinetic energy and magnetic energy (when included).We restart the simulation after introducing a dense cloud, with an overdensity  = 100 and radius  cl in the centre of the box while conserving the kinetic, thermal and magnetic energy density.This re-sults in an isobaric, cold, dense cloud with density  cold =  hot and temperature  cold = 4 × 10 4 K =  floor .As we use an adiabatic equation-of-state, the average temperature can increase significantly by the end of the turbulent driving phase, due to the dissipation of the turbulent energy.Hence, to bring the ambient temperature back to the required value, before introducing the cloud, we also rescale the internal energy of the whole box by a fixed constant.We also verify that this abrupt rescaling does not have any significant effect on the velocity distribution.
The input parameters for the driven turbulence are the kinetic energy injection rate ( ), the size of the simulation box ( box ) and the density of the medium ().Given a box size and gas density, we calculate the required  for a required turbulent velocity (see, e.g., Lemaster & Stone 2009), i.e., and assuming equipartition this yields where,  = 1 for hydrodynamic simulations and  = 2 for MHD simulations .Following the convention from previous studies, we also use cloud radii normalised by  shatter = min ( s  cool ).This corresponds to the  s  cool of the cold, dense medium in our simulations, i.e.  s  cool ( cold ,  floor ). .This shows the extent of amplification possible in the different cases.Even with a higher initial , turbulent motions can amplify the magnetic fields to lower .For cases with lower initial  strong cooling in the mixing layer can also lead to amplification.

RESULTS: TURBULENT RADIATIVE MIXING LAYER
Turbulent Radiative Mixing Layers (TRMLs) are mixing layer simulations that also include radiative cooling for the mixed gas.These have long been studied as an idealised small-scale setup for the mixing between different phases in a multiphase gas (e.g.Esquivel et al. 2006;Ji et al. 2019;Fielding et al. 2020;Tan et al. 2021;Yang & Ji 2023).In this section, we study the evolution of TRMLs for different cooling strengths and look for differences caused by the presence of magnetic fields.

Gas & magnetic field morphology
Previous studies of linear and non-linear evolution of mixing layers have shown suppression of Kelvin-Helmholtz (KH) instabilityinduced mixing, in the presence of magnetic fields.In the linear regime, the KH instability is suppressed for cases with magnetic field along the shear direction if Chandrasekhar (1961); Chapter XI, Eq. 205).While, in non-linear regime, the KH instability is suppressed for all magnetic field orientations (Ji et al. 2019).In this subsection, we reproduce these results and extend them by varying the cooling and magnetic field strength.
Fig 1 shows the temperature slices of the different TRML simulations, along with the magnetic field morphology, for different Alfvénic mach number (M A =  shear /  , where   is the Alfvén wave speed), and initial magnetic field orientation ( B0 ).We control the Da =  turb / cool,mix (where  turb = / turb and  cool,mix is the cooling time evaluated at  = 2 × 10 5 K, and  turb is calculated using scaling relations from Tan et al. 2021; we will investigate the role of  turb more below) by varying the box sizes, and M A by changing the initial magnetic field strength, for a given sonic Mach number.We have simulations with M shear,hot fixed to 0.16 and 0.33, corresponding to a  shear of 50 and 100 km s −1 .We find that the amplification in the magnetic fields is very different depending on the initial magnetic field orientation as we discuss below.
Fig. 2 shows the (=  thermal / magnetic ) and has layout similar to that of Fig. 1.The extent of amplification in the different cases is much clearer in Fig. 2, where the darker regions correspond to a higher magnetic field strength and lighter regions to a lower magnetic field strength.
The upper and lower row of Fig. 1 shows the fast and slow cooling cases, respectively.On one hand, in the fast cooling (Da > 1) a sharp temperature edge between  < 10 5 K and  ∼ 10 6 K gas is visible, i.e., a true multiphase structure exists, while on the other hand, a large amount of this "intermediate temperature" gas is visible in the slow cooling (Da < 1) case -with the exact amount depending on the suppression of mixing caused by magnetic fields.Unsurprisingly, in the HD case most mixed gas exists and generally in the MHD case with M A ≈ 1 the least.What is maybe a bit more surprising is the effect (and the evolution of) the magnetic field topology: in the B ∥ x case the suppression of mixing is easy to understand and expected from linear theory (Chandrasekhar 1961).However, we also find in all the other cases a (strong) suppression.B ∥ ŷ also has a strong  8)), Bottom Density profile for the different cases shown above.This shows the difference in the stability of the mixing layers for different cases.The cases with higher magnetic field strength along the shear (initial or amplified) are more stable.
The density profile shows the extent of the mixing layers.
effect, particularly in the M A ≈ 1 simulations.This is due to the amplification of the magnetic fields in the direction of the shear, as discussed below.On the other hand, for B ∥ ẑ initially, one can note two distinct effects depending on the value of M A .For M A > 1 the flow bends the magnetic field lines, resulting in a similar situation as in the B ∥ x case; in fact, even larger suppression since the bending of the field lines leads to a   >  initial .An artefact of this bending can be seen in magnetic field topology in the bottom panel for M A > 1 with B ∥ ẑ.We find a kink in the magnetic field moving downwards at the Alfvén wave speed.For M A ≲ 1, however, the magnetic field lines are so stiff that a bending by 90 degrees is not possible.Instead, we end up with diagonal field lines which, nevertheless, substantially suppress the mixing.
To better understand the exact order of amplification, we first consider the cases where the shear is super-Alfvénic (M A ∼ 10) in the hot medium, in the three central (2 nd -4 th ) columns of Fig. 1 & 2, to explain the extent of amplification in the different cases.
•  initial || ẑ || ninterface (Ma10_Bz): The amplification is the highest for this case.The Alfvén wave velocity in the dense medium is lower by a factor of  1/2 , hence the field lines are more "anchored" in the cold gas, compared to the hot gas.This causes field lines to bunch up near the interface, and result in high amplification of magnetic fields in the direction parallel to the shear.This amplification is so high that the magnetic field strength can get much higher than the initialised magnetic field strength.
• For  initial || x || ì  shear (Ma10_Bx): As the Kelvin-Helmholtz (KH) instability grows, it gives rise to vortices around ŷ. Fig. 1 shows that vortices can stretch and bend the magnetic fields, leading to their amplification.These vortices can further become turbulent and cause more amplification due to local dynamo effect.All of these put together, result in the second-highest amplification of magnetic fields along the shear direction.
• Lastly, for  initial || ŷ ⊥ ì  shear ⊥ ninterface (Ma10_By): The amplification is the lowest as the only process for amplification of magnetic fields is due dynamo effect from the turbulent motions generated in the mixing layer due to non-linear evolution of the KH instability.
This results in a general order for the magnetic field strength along the shear direction in Super-Alfvénic TRML simulations as: Ma10_Bz > Ma10_Bx > Ma10_By.
Similarly, for the cases where the shear is sub/trans-Alfvénic (M A ∼ 1) in the hot medium, in the three rightmost columns in Fig. 1, we again check for magnetic field strength along the shear direction.
• For  initial || ẑ || ninterface (Ma1_Bz): As the shear is sub-Alfvénic, the amplified magnetic field in the shear direction is not high enough to surpass the initialised magnetic field in the shear direction for  initial || x || ì  shear (Ma1_Bx).So, Ma1_Bz ends up with the second highest in the order of magnetic field strength along the shear direction.
•  initial || x || ì  shear (Ma1_Bx) has the strongest magnetic field in shear direction, just due to the high initial magnetic strength.
• This leaves the  initial || ŷ ⊥ ì  shear ⊥ ninterface (Ma1_By): Due to the much higher overall magnetic field strength, it is harder for the resulting turbulent velocity to cause any amplification.
Hence, we get an order for the magnetic field strength along the shear direction in Sub-Alfvénic TRML simulations as: Ma1_Bx > Ma1_Bz > Ma1_By.In our simulations, we find one exception to this order, at intermediate Da, where the order of Ma1_Bz and Ma1_By is switched.

Cooling rates
According to the linear KH instability criterion, the stronger the magnetic field in the direction of the shear, the more stable the perturbation gets.This means stronger magnetic fields in the later non-linear phase may disrupt the generation or cascade of further vortices.To test this, we quantify the stability of KH perturbations, using the linear stability criterion (Chandrasekhar 1961).So, the perturbation is stable to KH instability, if  The first two panels from the left show the difference between the direction parallel to magnetic fields versus the other directions.The third panel shows the outlier case of magnetic fields normal to the interface, where both the normal ( ẑ) and shear direction show much higher fluctuations due to the presence of magnetic fields along both these directions.Hence, we choose the directions which free from these spurious fluctuations in these different cases, as denoted in Eq 13 which can be expressed as a dimensionless number We calculate  KH using profiles of all the relevant quantities along the normal to the hot/cold interface.This results in a profile of  KH , in which a value <1 denotes a tendency towards stability while a value >1 shows a tendency towards instability.Fig. 3 shows the profiles of this KH stability criterion at an advanced stage of evolution, for different magnetic field orientations and strengths.The top and middle panels show the KH stability criterion and the bottom shows the corresponding density profiles.We align the profiles so that the hot/cold interface aligns between all the cases.We do not plot the points on the profile which have a Δ < 10 −2 km/s, or if the points are out of the computational domain.We find that due to the amplification of the magnetic fields in the shear direction, the KH instability is suppressed.The order of the extent of suppression seems to follow the same order as the amplification of the magnetic fields.For both M A > 1 and M A < 1, Ma1_By and Ma10_By are the most unsuppressed as the  KH is almost entirely in the unstable regime.For the other two directions, the extent of suppression depends on the size of the portion around the mixing layer that is stable.This means, Ma1_Bx and Ma10_Bz are the most suppressed in M A < 1 and > 1, respectively.This trend in suppression is important as this can affect the cooling rates, which we check next.
We study the cooling rates using the surface brightness of TRML simulations for all the different cases, that is, different values of the Dahmköhler number (Da =  turb / cool ), Alfvénic mach number (M A = /  , where   is the Alfvén wave speed), and initial magnetic field orientation ( B0 ).We define the surface brightness as the total luminosity per unit surface area, that is,  =  total /(    ).Fig. 4 shows the calculated saturation surface brightness (along with 2 errorbars) for the different simulations on the top row, along with its temporal evolution on the bottom row.
We find that the order of amplification of magnetic fields along the shear direction, mentioned above, matches the order of suppression of surface brightness, as shown in the right panel in Fig. 4.This is due to the higher suppression of KH instability by the higher magnetic field strength along the shear direction, as also expected from the linear theory (Chandrasekhar 1961).Below, we dive deeper into this correlation.
We also find that the difference in cooling rates, due to this suppression of KH instability, is reduced for the low Da cases with M A ∼ 10, as shown in Fig. 4.This might be due to the change in the rate-limiting process.For a low Da the cooling is very slow, so the cooling rate is bottlenecked by the slow cooling rate rather than the mixing rate.This does not happen for M A ∼ 1 because the mixing rate is suppressed to such low values that mixing continues to be the rate-limiting process.

Turbulence velocity profiles
Among the different ways to mix two phases, turbulence is one of the most efficient ones.In this section, we quantify the extent of turbulence in the mixing layer in the above discussed TRML simulations and look for connections of turbulence with the rate of mixing and cooling in the system.
This dependence has been studied before, in the absence of magnetic fields.Assuming a constant pressure and cooling function for hydrodynamic TRML simulations, as shown in Tan et al. (2021), for a strong cooling regime (Da > 0.1), and for slow cooling regime (Da < 0.1), Our next step is to check these relations in the presence of magnetic fields.We use the geometric method to calculate the  ′ , similar to Tan et al. (2021).First, we calculate this bulk velocity profile as the density-weighted average of the velocity field along the other two perpendicular axes (i.e.x and ŷ).Then, we subtract the bulk velocity profile along ẑ, which is normal to the hot/cold interface, from the velocity field to obtain the turbulent component.This turbulent velocity field allows us to compute density-weighted RMS values of this field along the axes perpendicular to the hot/cold interface normal (i.e.x and ŷ), to obtain 1D profile of all three components of turbulent velocity along normal to the hot/cold interface ( ẑ).Fig. 5 shows an example of the calculated turbulent velocity profiles for a snapshot where the cooling rate has reached saturation, for different initial magnetic field orientations.There are other methods like Gaussian filtering (e.g.Brereton & Kodal 1994;Adrian et al. 2000;Abruzzo et al. 2022) to get these 1D profiles, but we find that the choice of method does not significantly influence the next steps, as shown in appendix A.
Unlike Tan et al. (2021), we cannot select a particular direction that is untouched by shear or cooling inflow, as that direction can be aligned with the magnetic field.Hence, we have to calculate the turbulent velocity using different directions for different cases.In addition to that, the turbulent velocity component along the magnetic field can have some contributions from large velocity fluctuations along the magnetic fields, as Fig 5 shows.To get around this issue, we calculate the  ′ using the other two components perpendicular to magnetic fields, except in the case of  initial || ẑ || ninterface (Ma1_Bz, Ma10_Bz).For the exceptions, where the large-scale magnetic field orientations are along the shear and normal to the interface (i.e. x and ẑ), we only consider the turbulent velocity component along the direction perpendicular to those directions, ŷ.In short, we use the following expressions to calculate the turbulent velocities in the Figure 7. Left column Density rendering at 2.6 eddy after the cold gas cloud of size 310 shatter is introduced in a turbulent box with rms velocity of Mach 0.5, 2 nd -4 th column Density projections of the same simulation, at 1.3, 2.0 and 2.6 eddy after the cold gas cloud is added.The top row panels are from the simulations with magnetic fields, and the bottom row panels are from the simulation without magnetic fields.These show the clear differences between the morphology of cold gas of gas with and without magnetic fields, while also showing the similarities in the overall evolution of the cold gas.different cases, We, furthermore, checked that the turbulent velocity components used in the equations above are the ones that have similar profiles among themselves in each case, to ensure the isotropicity of the turbulent components, as shown in Fig. 5.As in Tan et al. (2021), we consider the maximum of the obtained turbulent velocity profile as the turbulent velocity.We repeat this process for every snapshot and consider the mean of the turbulent velocity over the last 5   as the saturated turbulent velocity ( ′ ) and the standard deviation as the error.
Note that throughout we denote the Damköhler number with Da when the measured  ′ was used and Da when the theoretically expected  ′ from the Tan et al. ( 2021) was used.

Turbulence vs. cooling rates
We use the above obtained  ′ to check the relations in Eq. ( 9)-( 10), as shown in Fig. 6.The top panel of Fig. 6 shows a scatter plot of the surface brightness and turbulent velocities calculated from the simulations.We find the respective strong and weak cooling scaling relation according to Eq. ( 9)-( 10), regardless of the magnetic field strength and orientation.For a better comparison, we remove the Da dependence and show the correlation in the bottom panel of Fig. 6.We clearly show that regardless of the initial magnetic field orientation or strength, the Eq.9-10 holds true, even though the relations originally obtained for hydrodynamic systems (Tan et al. 2021).We also confirm that the method of  ′ calculation does not affect these results, as shown in Appendix A, in Fig. A2.

RESULTS: TURBULENT BOX
In the previous section, we examined the TRML setup which is considered a more idealized version of Turbulent boxes (cf.§2), and found that magnetic fields can suppress the mixing in general, regardless of their initial orientation.If we follow this conclusion, one would expect the inclusion of magnetic fields in turbulent boxes to cause significant differences in the evolution of a cold cloud.This effect can manifest either as a change in the cold gas growth rates or a change in the survival criteria.In this section, we show results from the "turbulent box" setup in which we place a cold gas clump of size  cl in a turbulent medium, with a turbulent Mach number of M s , either with (MHD) or without (HD) magnetic fields.We then look for the effect of the magnetic fields, not only on the growth rates and survival criteria but also on the morphology and overall behaviour of the cold gas.

Cold gas survival and growth
When cold gas is subject to turbulence it can either be mixed away in the hot material or the mixed gas cools sufficiently fast to ensure continuous survival of the cold gas.Gronke et al. (2022) studied this effect using hydrodynamical simulations and found a relation between the critical value of  eddy /  cool,mix (equivalently  cl / shatter ) for a given turbulent velocity.However, as we have shown in the last section, magnetic fields can suppress Kelvin-Helmholtz (KH) instability-induced mixing via , between the hot and cold phases in a TRML.Hence, one could expect a similar significant difference in a turbulent box with (MHD) and without (MHD) magnetic fields.
As mentioned in (cf.§2), the initial seed magnetic fields are such that plasma  (=  thermal / magnetic ) ≈ 100.But, due to the local dynamo effect (Schekochihin et al. 2001), the magnetic field gets amplified to reach equipartition with the turbulent kinetic energy by the end of the driving phase of the turbulence and before the cloud is introduced.We can use the fact that at equipartition, M A ∼ 1,  Gronke et al. (2022).This shows the surprising lack of difference between the survival criterion, with and without magnetic fields.The subsonic turbulent simulations agree well with the previously found survival criterion, with some deviation in trans-sonic turbulent boxes (c.f.§ 4).and the relation between the sonic (M s ) and Alfvénic (M A ) Mach numbers as follows, to get an estimate of the final plasma  for a given M s , This means the plasma  at equipartition, when the cloud is introduced, is ≈ 16, 4 and 1 for M s = 0.25, 0.5 and 0.9, respectively.
We study multiple of these turbulent box simulations at different turbulent velocities (M s ≈ 0.25, 0.5 and 0.9) and with multiple cloud radii near the critical radii found in previous hydrodynamic studies (Gronke et al. 2022).See §2 for an overview of the setup.
Fig. 7 shows one of the HD-MHD pairs of simulations with M s ≈ 0.5 and  cloud ≈ 310 shatter , where the upper row shows the simulation with magnetic fields (MHD) and the lower row shows the same simulation but without magnetic fields (HD).The leftmost column shows renderings of the number density with the view in the direction of one of the diagonals of the computational domain, at a time 2.6  eddy after the cloud is introduced.The three columns on the right show density projections of the same simulations as the first column, but at different times (1.3  eddy , 2.0  eddy and 2.6  eddy after the cloud insertion).We can see that in both the HD and the MHD simulation the crude behaviour of the cold gas is similar.The cloud survives and grows in both simulations, and the overall amount of gas in the simulation box also looks roughly similar.Fig. 7 also shows how the gas structure evolves.The cold gas seems to grow as it gets stretched, folded and transported by the turbulent motion in the hotter surrounding medium.We also see the difference in the general morphology of the cold gas in the two cases.The cold gas morphology is much more filamentary in the MHD simulation, while it is very clumpy and less dispersed in the HD simulation.We will discuss the morphology further in § 4.2.
Next, we check for the growth (or destruction) of the cold gas in all the turbulent box simulations.We define the cold gas mass as the total gas mass with temperature below 2  floor = 8 × 10 4 K and normalize the obtained value with the initial cold gas cloud  2022).This shows that there are only marginal differences between the growth and destruction rates of the cases with and without magnetic fields, compared to the differences seen in the TRML simulations.The differences further diminish as we consider cases well within the survival regime.mass.We take the obtained cold gas mass and check for survival or destruction at different Mach numbers and cloud radii.We determine the survival or destruction of the cloud using the final normalized cold gas mass values.The cases with final  cold / cold,0 > 1(< 1) are considered to show cold gas survival (destruction).We plot this survival or destruction for the different cases as a scatter plot in Fig. 8.It clearly shows the lack of difference in survival criteria between the pairs of simulations with (MHD) and without (HD) magnetic fields.This shows that the inclusion of magnetic fields does not affect the survival or destruction of the cold cloud.We also plot the survival criteria found by Gronke et al. (2022) in Fig. 8, given as We find that this survival criterion works well for our subsonic turbulent simulations, but our transonic turbulent simulations seem to deviate slightly from this survival criteria.This could be due to the difference between subsonic, transonic and supersonic turbulence due to the presence of shocks in later cases, possibly destroying the clouds which would have survived in the absence of these shocks.Regardless, this does not affect our original conclusion about the lack of significant difference in survival (or destruction) between HD and MHD simulations, hence, we leave the investigation for causes of this discrepancy to future studies.Another property which can have differences due to the inclusion of magnetic fields is the growth (destruction) rates of the cold gas.For that, we repeat the process of calculating the cold gas mass for each snapshot to obtain the temporal evolution of cold gas mass for all the different simulations with M = 0.5, and Fig. 9 shows the same.We find a lack of drastic differences in the growth rates of the simulations that are well within the survival regime.We see slight differences (within a factor of 2) in simulations close to the transition regime, but it is still less than the order of magnitude differences seen in and expected from the TRML simulations.We also plot the expected growth curve for the surviving cases, using  grow from equation ( 7) in Gronke et al. (2022) and the mass growth equation for "fragmented" growth,  cold =  cold,0  (/ grow ) , as where,  = 0.5 is a fudge factor.This lack of significant difference in the cold gas mass evolution and survival criteria between turbulent box simulations with (MHD) and without (HD) magnetic fields is surprising, in view of the results from TRML simulations in the previous (cf.§3), where we found that the mixing rates are (highly) suppressed when magnetic fields are introduced into the same simulations.This dichotomy in the results can be confusing, and we discuss a possible solution to this in the discussion section (cf.§5.1).

Cold gas distribution and morphology
In the previous subsection, we showed how the presence of magnetic fields seems to have only a marginal effect on the cold gas mass growth and survival.Still, the magnetic fields are not entirely inconsequential.The magnetic fields can affect the gas flow and vice-versa due to effects like flux-freezing.We also saw in the above section that the morphology of the cold clouds is different between the simulations with (MHD) and without (HD) magnetic fields (cf.Fig. 7).In this section, we present such differences and quantify these differences in the morphology and distribution of the cold gas.
Turbulent transport has been a long-standing field of research in fluid dynamics.In a turbulent medium, the stochastic motions can transport, break, coalesce or mix the cold gas clouds.This results in a wide variety of cold gas cloud morphology.Gronke et al. ( 2022) calculated the cloud size distribution in a hydrodynamically turbulent medium and found a power law,  (> ) ∝  −1 (which has also been found in larger scale simulations, e.g., Tan & Fielding 2023).As we saw a significant difference in the cold gas structure in Fig. 7, we check if the visual difference in cold gas morphology between HD-MHD simulations is reflected in cold gas size distributions.We calculate the cumulative number distribution of the cold clumps in a set of HD-MHD simulations with M s,hot = 0.5 and  cloud = 310 shatter .Similar to Fig. 9, we define cold gas as the gas with temperature below 2 floor = 8 × 10 4 K and use feature labelling functionality in SciPy's (Virtanen et al. 2020) ndimage as the clump finding algorithm to identify the cold clumps.We determine the volume of the obtained clumps and use it to calculate the cumulative number distribution, shown in Fig. 10.We find that the distribution is well approximated by a power-law with slope -1, i.e.  (> ) ∝  −1 , while deviating at the two extremes of the volume range due to resolution limits at lower volumes and statistical under-representation at higher volumes.This matches the results from Gronke et al. (2022).We also see that the number distribution does not show a drastic difference between the HD and MHD simulations, apart from a slight excess of small clumps in the HD simulation.This means the cold gas clumps in MHD simulation are not significantly smaller or larger in volume than its HD counterpart.
Even though the difference in number distribution is minor, Fig. 7 clearly shows significant morphological differences between the cold gas structure in HD and MHD simulations.Visually, the cold gas in MHD simulations has a much more filamentary shape, while it has a more clumpy cold gas morphology in the corresponding HD simulations.We quantify this filamentariness of the cold gas structures as the length of the longest "shortest path" possible within the clump.
To do so, we first identify the individual clouds, as done for Fig. 10, and create a "neighbourhood graph" for each clump using an adjacency matrix.In this "neighbourhood graph", each gridcell inside  ) for hot (green line) and cold (blue line) gas phases in a set of simulations with (MHD) and without (HD) magnetic fields, at M = 0.5 and  cl = 310 shatter .The dashed and solid lines show the VSF at different times,  = 1.32 eddy and 3.95 eddy after introducing the cold gas cloud.This shows the decreasing difference in VSF of the two phases with time, in both cases, which means that the two phases are kinematically well-connected.
We also find a smaller early-time difference between the hot and cold gas VSF for the MHD simulation, indicating a better kinematic connection in that case.
the cloud is a node and two nodes are connected with an edge if the two share a face.We calculate the shortest path between every node in this neighbourhood graph, and take the longest from this list of shortest paths as the required longest "shortest path".As many of the largest clumps contain ≳ 50, 000 gridcells, we have to use a faster way which can give a close enough answer instead of using the brute force method.The slowest step of the method is the shortest path calculation among each pair of nodes.Hence, to speed up this step, we only consider every 4 th node to identify the two points with the longest "shortest path", and later recalculate the path using the full graph with all nodes to get the final length.This optimised method drastically reduces the number of path calculations and makes this analysis computationally feasible.We also find negligible differences between the full brute force method and the optimised method in our tests.This is because the points with the actual longest "shortest path" is usually in close vicinity, likely within 4 grid lengths, unless the reduced graph is drastically different, which is rare.
We repeat this process for each clump in the MHD and HD simulations with M = 0.5 and  cloud = 310  shatter and plot the histogram of the obtained longest "shortest path" in Fig. 11.We find that the 90%ile of this length distribution for MHD simulations is longer by a factor of 2 compared to the corresponding HD simulation.Assuming, a constant volume of a cylindrical cold gas clump, which is reasonable as shown in Fig. 10, an increase of 2 in length corresponds to an increase of ≈ 2.8 in the length-to-width ratio of the clump.This method under-quantifies the filamentariness of the clumps, as the connected filamentary structures that are shorter than the main filament are not included.A full tree-based filament analysis will be the ideal method for this analysis, but we leave the detailed study of filamentariness to future investigations.
We find that the p-value for the two length distributions in HD and MHD is lower than 0.05, which means we can consider the filamen-tariness of HD and MHD simulations to have different underlying distributions.The KS statistic quoted in Fig. 11 quantifies the difference between the two distributions and is linked to the p-value.The higher the KS statistic, the lower the p-value.As our number of samples is limited by the number of clouds in the simulation, we will need a bigger box and longer runtime to improve the confidence level of this conclusion.
So, we conclude that, even though the cold gas clouds in the MHD simulation are similar in volume and its distribution, compared to their HD counterparts, they are significantly more filamentary in their morphology.

Cold gas entrainment
In a multiphase environment, the motion of one phase can affect the motion of another via drag forces or mixing-induced momentum transfer (see, e.g., Gronke & Oh 2020a;Tonnesen & Bryan 2021).This means the phases can be kinematically linked.In addition to that, due to flux freezing in the gas, the magnetic fields can increase the extent of this kinematic link.A good way to check for this is the first-order velocity structure function (VSF).It quantifies the average difference in velocities of gridcells separated by a given distance.The difference in the VSF of the two phases corresponds to a lack of link in the kinematics between the two, while a smaller difference corresponds to a greater kinematic link.
We calculate the velocity structure function (VSF) of the hot and cold phases of gas in simulations with (MHD) and without (HD) magnetic fields at different times.First, we calculate pairwise distances and velocity difference magnitudes between each pair of gridcells.Again, due to computational constraints, we cannot use all cells for the pairwise calculations, so we randomly choose 2 × 10 4 gridcells for this calculation.Then, the list of pairwise velocity differences is  12).The shaded regions show the corresponding 15-85%ile intervals.The figure also shows the shear on the clump boundaries is about an order of magnitude lower than the turbulence velocity in the simulations.Also, on average, clumps in the MHD simulation seem to have a marginally lower, but very similar shear, in comparison to HD simulations.
binned according to the pairwise distances and we plot the average of the velocity differences in each bin as the VSF in Fig. 12.In general, we find a higher value of VSF in the hot gas than the cold gas for both the simulations at all different times in the evolution, as seen in other idealised simulations (Gronke et al. 2022;Mohapatra et al. 2022a), and even some observations (Li et al. 2022).
For structure function calculations with steep slopes, Seta et al. (2023) had found a 2-point stencil to be unconverged and suggested the use of higher-order stencils.But, the slope of VSF for Kolmogorov turbulence (1/3) is shallow enough for a 2-point stencil to be converged.Hence, we use a 2-point stencil for all our VSF calculations.
This difference between the hot and cold VSF is much larger in HD simulation, on the left panel of Fig. 12, while in MHD simulations (right panel of Fig. 12) the difference is much more subtle with both hot and cold VSF comfortably within 16-84 %ile range of each other.This shows that the cold phase is, in general, better entrained in MHD simulations compared to HD simulations.This is likely due to the flux-freezing of the magnetic fields that can result in a more efficient kinematic connection between the hot and cold phases, as mentioned before.In presence of flux-frozen magnetic fields, any relative motion between the two phases encounters an enhanced drag force (McCourt et al. 2015).We also find that, even though the VSF of the different phases start off differently in HD-MHD simulations, they end up with very similar hot and cold VSF profiles in both cases.This means, given enough time, both simulations reach a similar extent of entrainment.Still, we do note that the entrainment is faster in MHD simulations, compared to HD, as shown by the faster decreases in the difference between hot and cold medium VSF in MHD simulations.This result indicates high entrainment of cold gas in hot gas, albeit an imperfect one.More importantly, it also points to a more efficient and faster entrainment of the cold gas in the presence of magnetic fields, with both HD and MHD simulations reaching an equivalent entrainment state, given enough time.
Further, we calculate the average shear at the cold gas clump boundaries for each snapshot in the HD and MHD simulations with M = 0.5 and  cl = 310 shatter using yt (Turk et al. 2011).Fig. 13 shows the evolution of this average and the 15-85%ile interval of the shear at clump boundaries with time.We find that the shear at the clump boundaries is, in general, about one order of magnitude lower than the turbulent velocity in the simulation.This again indicates a high entrainment of the cold gas.Also, the slightly lower values of average and 85%ile value of shear for MHD simulations suggest a more efficient and faster entrainment in the presence of magnetic fields.

Magnetic fields strength and structure
In MHD simulations, magnetic fields are kinematically very important, as the gas flows affect the magnetic fields and vice-versa.Apart from affecting the kinematics, magnetic field structure can also affect other processes like thermal conduction and cosmic ray transport in an astrophysical media (e.g., Kempski et al. 2023;Ruszkowski & Pfrommer 2023).
The turbulent motions can result in a local dynamo effect, leading to amplification of magnetic fields in MHD simulations.The extent of this amplification can vary in the different phases due to differences in the Alfvénic wave speed ( A = / √ ) and turbulent velocities.On top of the dynamo effect, the compression of gas during its cooling can also cause amplification during cold gas formation, due to fluxfreezing.
To examine these differences, we check the distribution of magnetic field strengths in the different phases.We define the cold phase as the gas with temperature  < 2 floor = 8 × 10 4 K, hot phase as  > 0.5 amb = 2 × 10 6 K, and mixed gas as the gas with temperatures in between, i.e. 8 × 10 4 K <  < 2 × 10 6 K. Fig. 14 shows the distribution of magnetic field strengths in these three gas phases, for simulations of three turbulent Mach numbers and with two cloud radii.The exact distribution has a non-trivial dependence on factors like the turbulent Mach numbers and cold gas growth rate.But in all cases, the mixed and cold gas magnetic strength distributions are centred at stronger magnetic fields, while the hot gas magnetic strength is centred around weaker magnetic fields, with the dashed line showing the equipartition magnetic field strength.This higher magnetic field strength in cold and mixed gases could be due to three possible processes: turbulent local dynamo in the dense gas as the equipartition magnetic field is higher for a denser gas mov-ing at similar velocities, flux-freezing accompanied by compression due to cooling, and magnetic draping caused by the relative motion between the dense gas and magnetic fields.We discuss more about these processes, and possible order of importance in the discussion section (c.f.§5).
Flux-freezing and the turbulent motions can result in a very tangled magnetic field structure.These tangled mangetic fields can have many consequences including reduced thermal conduction and slower cosmic ray transport.Presence of multiphase gas in a turbulent medium can add further complexity to the magnetic morphology, due the magnetic field strength distributions in different phases, as we show earlier in this section.To better understand this, we study the structure of the magnetic fields using magnetic field streamlines.We use yt (Turk et al. 2011) to calculate 10000 streamlines for 100 different streamline lengths ( stream ) between 0.01 − 1 box .We calculate the displacement between the two endpoints of the streamlines (), to get the ratio  stream /.This ratio denotes the extent of entanglement of the magnetic field.A  stream / = 1 indicates a perfectly untangled streamline, with higher values denoting a higher extent of entanglement.We repeat this process for different  stream / box , and obtain distributions of the extent of entanglement ( stream /) for each  stream / box .We calculate the mean, median and mean of the logarithmic lengths in each  stream / distribution.
The upper panel of Fig. 15 shows the trend of mean, median, mean(log) and 15-85%ile interval for each  stream / box value.We find that all the metrics of ensemble value of the ratio  stream / increase linearly with  stream / box .The upper panel of Fig. 15 shows the close approximation of the linear trend for mean and mean(log).This means the extent of the entanglement increases linearly with the length of the streamline.This property could be a sign of fractal-like behaviour of the field lines down to a certain threshold at small scales.We discuss this further in the discussion section (cf.§5).Note that the asymmetrically located 15th and 85th percentiles with respect to the median indicate a long tail towards longer  stream /.
To show this more explicitly, we further choose three different streamline lengths, shown as vertical dotted lines in the upper panel of Fig. 15, and recalculate the streamlines for 10× more (10 5 ) starting points.We repeat the above mentioned process to calculate the  stream / values and calculate the probability distribution function of the entanglement,  stream /, for these three values.We plot these probability distribution functions as solid lines in the bottom panel of Fig. 15.Even though we find some minor deviations at higher entanglement values due to insufficient counts caused by the reduced number of samples, the histogram at lower entanglement values is robust and fairly well converged with the number of streamlines.
Next, we attempt to find an analytic form for the different probability distribution functions (PDF) that we found earlier for the  stream /, in the bottom panel of Fig. 15.We calculate the variance () and mean () of log 10 ( stream /) and find a strong linear relation between the two as  = 0.242( − 0.031).For a Γ distribution, the variance and mean are given by  =  2 and  = , where  and  are the shape and scale parameters, respectively.This means, our  −  relation for log 10 ( stream /) matches the properties of a Γ distribution with  = 0.242, and  = / = /0.242.Using the linear fit for the , shown in upper panel of Fig. 15, and the equation of the Γ distribution, a fit for the PDF of  stream / for a given  stream / box is a Γ distribution for  = log 10 ( stream /) − 0.031 with  ≈ log 10 (1.5 stream / box + 1) and the  = 0.242 mentioned above.
The bottom panel of Fig. 15 shows the analytical form with the dotted lines.We find that the analytical form agrees very well with the PDF for long streamlines.But for shorter streamlines, at intermediate values of entanglement ( stream /), it overestimates the PDF at intermediate entanglement values in the tail.This can be an indication towards a different underlying analytic form for PDF, which is equivalent to the Γ distribution at longer streamlines.Or, it can be due to resolution effects as they start to become more important for short highly entangled streamlines.We leave a deeper investigation of this for future studies.
As the charged particles tend to gyrate around and follow the magnetic field lines, analytic form for the magnetic field morphology, similar to the ones we find above, can be used in developement of models for their transport through a multiphase turbulent medium.

Synthetic absorption line spectra
As shown in the results section of turbulent boxes (cf.§4), we know that the morphology and details of kinematics can differ significantly between the simulations with (MHD) and without (HD) simulations.This difference can affect observational probes like predicted quasar absorption line spectra, because the column density and Doppler shift, the two major features of lines, can be affected by these differences in morphology and kinematics.We investigate these effects and their observational consequences in this section.
Fig. 16 shows the distribution of column densities of cold ( < 10 5 K) and intermediate/mixed (10 5 <  < 10 6 K) along one of the dimensions of the box for a set of HD-MHD simulation with M hot,turb = 0.5 and  cloud = 310  shatter , same as Fig. 7 and 12.Note that because  shatter ∝  −1 McCourt et al. ( 2018) the column densities simulated can be directly compared to observations.We still find that the column density histogram of the cold gas shows a greater extent of difference between the set of HD and MHD simulations, compared to the intermediate gas.We also see that most of the difference shows up at lower values of column density.So, we expect to see some difference between the HD and MHD simulations in observational probes of cold gas, for example, LOS absorption due to MgII, at lower equivalent widths.
In order to investigate the observational consequences of the differences between HD and MHD simulations, we create mock LOS absorption spectra using Trident (Hummels et al. 2017) on the same set of HD-MHD simulations as Fig. 7, 11, 12, and 16 (M hot,turb = 0.5,  cloud = 310  shatter ).Note that both of these snapshots have similar cold gas volume filling fractions.
We sample ∼ 10000 line-of-sight (LOS) spectra along one of the axes of the computational domain on a 100 × 100 grid.Due to the isotropic nature of the system, the particular choice of the axis should not affect the statistical inferences.We use a Δ = 0.1Å (corresponding to a spectral resolution  ≡ /Δ ≈ 28, 000) to create the mock absorption features for MgII at uniform solar metallicity.We select spectra which have a maximum absorbed flux of more than 0.1, over a continuum flux of 1.0, so that we are only considering LOS that pass through significant amounts of cold gas.We also exclude the spectra which have saturated features (with a flux less than 0.1) because they correspond to unnaturally large cold gas volume filling fractions, thus, leaving the number of components ill-defined.Next, we calculate the equivalent width (EW) for the MgII line complex and the number of absorption features for each LOS, as the area under the continuum and the number of minima in a spectrum, respectively.Fig. 18 shows the 2D distribution of EW and the number of absorption features for the HD-MHD simulations.We also show the fit obtained for the same quantities from observations of MgII absorbers in quasar spectra in Churchill et al. (2020).Interestingly, the more frequent regions of both the 2D histograms in Fig. 18 are close to the relation found in Churchill et al. (2020).We also find Frequency/N cells 0.25 R cloud = 77 l shatter t/t eddy = 3.92 Histogram of magnetic field strength in gas within different temperature ranges, namely hot ( > 2 × 10 6 K), mixed (8 × 10 4 K <  < 2 × 10 6 K), and cold ( < 8 × 10 4 K) gas, for two simulations where the cloud gas cloud survives,  = 3.92 eddy after its introduction.Left M = 0.5,  cl = 310 shatter .Right M = 0.25,  cl = 77 shatter .The dashed vertical line corresponds to the equipartition magnetic field strength, achieved in the hot ambient gas at the end of driving the turbulence.This shows that the magnetic fields are significantly amplified as the gas cools down to a lower temperature.We discuss the possible causes of this amplification in § 5.2.that, for the same EW, the MgII spectra of the MHD simulation tend to have a slightly higher number of absorption features, compared to the HD simulation, but these differences are marginal.
We repeated this exercise with 10x better spectrum resolution, which is much higher than that of the observed spectra, and we also included an additional CIV 1551Å line.The corresponding distributions, analogous to Fig. 18, are shown Fig. C1 and C2 in Appendix C).In both cases, the distributions change, but still roughly follow the observed curve.

Mass transfer rates in a magnetized, turbulent medium
The two results for mixing layers and turbulent box simulations shown in Sec.2.1 and Sec. 4, respectively, present a dichotomy.On one hand, TRML simulations show a significant suppression in the mixing of two phases, and on the other hand, the turbulent box simulations show the lack of a similar difference in mixing rate, as shown by the cold gas growth rates and survival criterion.
To resolve this, we need the answer to the question, what is the primary mechanism of mixing of two phases?The mixing happens when the multiphase structures get small enough to reach the "dissipation scale" where the molecular diffusion is fast enough to mix the two phases (Obukohov-Corrsin phenomenology;Oboukhov 1949;Corrsin 1951).There can be multiple ways to reach such small scales, and one of these is via vorticity.Vorticity or vortices can stretch, fold and transport, in other words, it "stirs" and stretches the structures, eventually reaching the small scales where molecular diffusion can take over and "mix" the two phases (Villermaux 2019).In theory, this vorticity does not have to be part of turbulence, but in the high Reynold's number of fluids, as is the case in astrophysical mediums, the vortices generally become turbulent.This causes a faster stirring of the multiphase structures, and a more rapid increase in the surface area and decrease in structure width, resulting in more efficient diffusion.Hence, in an astrophysical medium, the main mechanism of mixing is expected to be turbulent mixing.In principle, the mixing should only depend on the turbulence and be independent of the source of turbulence.
In a Kelvin-Helmholtz or TRML setup, the initially structured vortices quickly give way to a turbulent mixing layer.The turbulence in this mixing layer is the key mechanism of stirring and mixing the two phases.This lead to this chain of processes: KH instability → Turbulence → Mixing.When magnetic fields are introduced in the system, depending on their orientation, they hinder the link between the KH instability and turbulence by slowing down the rate at which the turbulence in the mixing layer is driven.But as we show above, importantly, the magnetic fields do not affect the other link that connects turbulence with mixing.So, we expect to see a tight correlation between the turbulent velocity in the mixing layer and the extent of mixing/cooling that is occurring in the mixing layer, regardless of the magnetic field orientation or strength (cf.Fig. 6).Hence, we conclude that suppression of mixing in TRML simulations in magnetic fields is due to the reduced driving of turbulence in the mixing layer, which in turn leads to a reduced mixing of the two phases.
The situation is different in our turbulent box setup.There the system is the driven turbulence that cascades from the largest (boxsize) scales to the smallest (gridcell-size) scales.This implies that  ′ is fixed and since the mixing and cooling rate only depends on  ′ directly, we obtain similar growth rates in the HD and MHD cases -explaining the unaltered survival criterion and mass transfer rates found (cf.Figs. 8 and 9, respectively). stream / for different streamline lengths ( stream ).The dashed lines show the corresponding best linear fits and the shaded region shows the 15-85%ile interval.The general trend of increasing entanglement for longer and longer streamline lengths indicate a fractal-like structure of the magnetic field lines, discussed further in § 5.2.Bottom inset Points denote the mean and variance of log 10 ( stream /) and the green dashed line shows the linear fit,  = 0.24(  − 0.03).We use this relation to calculate the shown probability distribution.Bottom Solid lines show the probability distributions of different values of entanglement, log 10 ( stream /), for three values of streamline lengths.The dashed lines show the corresponding calculated Γ distributions, with the parameters mentioned in the legend.This shows the close agreement between the estimated and calculated probability distributions.There are some deviations for the probability distribution of small streamline length, which is discussed further in § 4.4.
We also find direct evidence that the turbulent, cascading  ′ is responsible for the mixing, and not the (also in the turbulent box present) hydrodynamical instabilities seeding smaller-scale turbulence.Firstly, the velocity structure functions of both the hot and the cold medium follow each other closely (cf.Fig. 12) indicating near perfect entrainment of the cold gas.Secondly, we also show explicitly the shear between cold and hot gas (c.f.§4.3) being small, i.e., the cold gas is well-entrained in the hot ambient gas.This means the shear is minimal, resulting in a lesser extent of turbulence in the This shows that the column densities for the above cases are within the observationally expected column densities for absorption spectra in a circumgalactic environment.It also shows that the lower end of column density distribution for cold temperature gas has a higher extent of difference between the HD and MHD simulations.This makes an absorption line tracing the cold gas a prime candidate for looking at observational differences between the HD and MHD simulations.
mixing layer between the two phases.If solely the shear would be responsible for the mixing and cooling, we estimate the mass doubling time to be ∼ 5 eddy for the turbulent box with M = 0.5 and  cl = 310 shatter , which is about an order of magnitude longer than actually found in the simulations (using the TRML scaling relations of Tan et al. 2021 for each surface cell on the cold gas clump).
In summary, we find the  ′ →  relation to be universal in HD and MHD (and consistent with high-resolution TRML studies; Fielding et al. 2020;Tan et al. 2021).However, magnetic fields prevent instabilities to form in the mixing layer setup leading to a lower  ′ and thus to a decreased mass transfer rate.When the extent of turbulence is fixed by larger scales -as done in the turbulent boxthe magnetic fields cannot suppress the mixing leading to comparable luminosities in the HD and MHD runs.
In realistic, astrophysical multiphase systems such as the ICM, CGM or ISM turbulence is also seeded on larger scales, then cascading downwards.In the ICM, for instance, AGN feedback is believed to play a dominant role in the stirring process leaving a characteristic imprint on the VSF (Li et al. 2022).Similarly, for the CGM where both (AGN and stellar) feedback processes as well as cosmological inflow act on ∼ 100 kpc scales 'stirring' the CGM (Chen et al. 2023).The alternative 'shearing layer' picture might occur in multiphase systems where bulk flows are dominant such as galactic winds and cold streams; however, since also there non-negligible turbulent is present which mixing channel is dominant is still unclear (Schneider et al. 2020;Tan & Fielding 2023;Rathjen et al. 2023).

Magnetic field amplification and morphology
We find that the magnetic field strengths in cold and mixed gas of our MHD turbulence simulations are higher than their equipartition values in the hot medium (c.f.§4.4).As discussed earlier, this higher value in the cold and mixed gas can be due to higher equipartition values in denser gas (as  eq ∝ √ , with the caveat of assuming similar turbulent velocities in hot and cold medium, which we discuss  16).The dotted black line shows the threshold of the minimum absorbed flux of a feature, and the red circles show the features that we consider for analysis.This figure is only for reference, as these are higher resolution spectra compared to the ones used in the analysis at Δ = 0.1Å, which is closer to observational spectral resolution.
later in the section), due to flux-freezing during compression from hot to cold medium (e.g.Sharma et al. 2010;Gronke & Oh 2020a), or due to magnetic draping around the cold gas clumps (Dursi & Pfrommer 2008;McCourt et al. 2015).
It is hard to disentangle these three processes as the extent of amplification in the simulation (≈ 6 0 ) can be achieved via all above the processes.The flux-freezing can cause an amplification up to  2/3  0 ≈ 22 0 , assuming an isotropic, isobaric collapse from  hot to  cold and conservation of magnetic flux.The local dynamo and magnetic draping can account for an amplification up to  1/2  0 ≈ 10 0 , assuming the amplification continues until equipartition is reached, i.e.M A ∼ 1, and similar  turb surrounding the cold gas means the new equilibrium magnetic field in the cold gas increases by √︁  cold / hot .The flux-freezing causes a higher magnetic field in newly formed cold or mixed gas, and the other two processes amplify the existing magnetic field in the cold or mixed gas.As the magnetic fields reach equipartition values, they start to become stiff to the gas motions and start to back-react and influence the gas motions.This means the amplification value of  1/2  0 ≈ 10 0 at equipartition gives a rough upper limit on the amplification by all the processes.And, this agrees with our results in Fig. 14.
Out of the possible processes, turbulent local dynamo and magnetic draping are less likely due to a few reasons.For the magnetic fields to be amplified to 10 0 due to turbulent local dynamo, the turbulent velocity at cold gas cloud scales has to be similar to the hot gas turbulent velocity.But, due to the small scales of the cold gas clumps, the turbulent velocities at cloud scales will be much lower at ∼  turb ( clump / box ) 2/3 .Hence, the local dynamo will not be able to cause the calculated high amplifications.
For magnetic draping to amplify the fields, there needs to be  and 17).The dash-dotted green line shows the relation found in Churchill et al. (2020).This shows that there are only marginal differences in the overall distributions of HD and MHD simulations, despite the differences in Fig. 16.We also find that they agree quite well with the observed relations from Churchill et al. (2020).
a significant relative velocity ( rel ) between the hot and cold gas, which generally is not the case, as we find a very similar VSF for hot and cold gas and low shear between the phases.This means, the  rel ≪  turb , hence the amplification of magnetic fields due to such process is probably insignificant.In addition, draping generally requires and leads to structured magnetic fields as they 'drape' around the clouds (Dursi & Pfrommer 2008) -something we do not observe in our simulations.This leaves flux-freezing and subsequent compression of magnetic fields as the only process that can cause significant amplification.Once, the amplification reaches a limit where the magnetic fields are stiff (trans/sub-Alfvénic), the gas continues to evolve along magnetic field lines, hence cold gas growth does not necessarily have to compress the magnetic fields.
Next, we consider the structure of the magnetic fields.In our study, we find that the extent of entanglement ( stream /) of the magnetic field lines increases linearly with the length of the streamline (c.f.Fig. 15).This points to a structure where the longer the streamlines are, the more relatively small-scale structures are sampled.This is possible if the magnetic streamlines have a "fractal-like" structure that goes on until a fixed small-scale, which is the grid-scale in our simulations.Hence, the longer the streamlines are, the wider the range of perturbations that are included, leading to the increasing trend of entanglement.This is analogous to the well-known problem of measuring a coastline, where the measured coastline length increases with decreasing length of the measuring stick.In this case, the roles are reversed.The measuring stick has a constant length, while we make the coastline longer.Assuming self-similarity, if we rescale this longer coastline, we effectively make the measuring stick smaller and we get back to the original coastline measuring problem.Let  ∝ 1/ stream be the effective length of the measuring stick.Hence,  =  stream / will be the rescaled coastline (streamline) length.We can use the relation found in Fig. 15 and the expression for the length of self-similar fractals, i.e.  ∝  1− (Mandelbrot 1967;Mandelbrot & Wheeler 1983) to find the fractal dimension () of the magnetic field lines.We find that the magnetic field lines have a fractal dimension,  = 2.Such naturally occurring fractal structures with a fractal dimension of 2 in 3D space are known to exist, with Brownian motion being one example (Falconer 1985).Previous studies of TRMLs have also found fractal structures, for example, Fielding et al. (2020) show that the cooling layer in a TRML has a fractal dimension of 2.5, while Tan et al. (2021) find a slightly different value but note that measured values can differ.
We also find that a Γ distribution on logarithmic entanglement (log 10  stream /) matches fairly well with the computed distribution from the simulations for longer streamlines.The Γ distribution does a poorer job for very short streamlines, which might hint towards a transition to or altogether a different distribution for the entanglement.Or, this might possibly be due to higher resolution effects on the shorter streamlines.
We hope that this analytic form for magnetic field entanglement will be helpful in development of models for transport charged particles through magnetised multiphase turbulence.

Connection to observations
Multiwavelength studies now allow the joined observational study of multiphase astrophysical media.Of the many ways to probe the properties of the multiphase gas, the absorption lines are one of the widely used methods (e.g., Steidel et al. 2010;Crighton et al. 2015;Chen 2017;Rubin et al. 2022).The different phases in the CGM of an intervening galaxy can deposit absorption features on the background quasar continuum.As the different sections of the absorbing medium can be moving with different velocities, these absorption features can be deposited at different Doppler-shifted positions near the line centre with different widths.Hence, the absorption features provide information about the kinematics and structure of the absorbing medium.
We find that there is no significant difference in mock absorption features of MgII 2796Å with and without magnetic fields.We, furthermore, show mock absorption spectra from both HD and MHD simulations agree with observed MgII absorption features from Churchill et al. (2003); Churchill et al. (2020), who established a relation between the number of 'absorbers' and the total equivalent width of the absorption.We also show in Appendix C that this agreement is approximately valid across spectral resolution and absorption lines.
As we found a universal clump mass distribution following d/d ∝  −2 (cf.Fig. 10 and  § 4.2) in both the HD and MHD cases (consistent with Gronke et al. 2022), this suggests that the Churchill et al. (2020) is a direct consequence of the clump mass distribution, and similar probes might be used to constrain it providing an interesting avenue for future work.
In addition to absorption lines, there are many studies that investigate the emission lines from multiphase media.Li et al. (2022) look at the multiphase turbulence in the ram-pressure stripped tail of ESO 137-001 using different emission lines.They find a similar velocity structure function as ours (in Fig. 12) and many other simulations (Mohapatra et al. 2021(Mohapatra et al. , 2022a)).This shows that both simulations and observations point towards a high extent of kinematic coupling between the different phases in astrophysical media.

Connection to previous studies
Due to the very high Reynolds number of astrophysical media, they are highly susceptible to turbulence.Hence, these media are expected to be turbulent in all the different scenarios in which energy is being injected into the medium, be it via supernovae, accretion or mergers.This turbulent nature of astrophysical medium has been studied before in previous studies (Schekochihin & Cowley 2007;Lancaster et al. 2021;Hu et al. 2022;Li et al. 2022;Federrath 2013;Elmegreen & Scalo 2004;Wittor & Gaspari 2020).There is also a plethora of studies that look at the different aspects of magnetohydrodynamic (MHD) turbulence, both in contexts related and unrelated to astrophysical mediums (see, e.g., review by Schekochihin 2020).
Recently, there has been a significant focus on the multiphase nature of such turbulence, with or without magnetic fields.Previous studies like Mohapatra et al. (2022c) and Gronke et al. (2022) have looked into hydrodynamic multiphase turbulence, while studies like Mohapatra et al. (2022b) and Mohapatra et al. (2022a) investigate the same with magnetic fields.And, studies like Seta & Federrath (2022) have looked at the evolution of magnetic fields in a multiphase medium.The key difference between these studies (except Gronke et al. (2022)) and ours is the thermal instability of the ambient hot medium.In our setup, we mimic a heating source and turn off the cooling for gas hotter than 0.5 amb , hence the ambient hot medium is thermally stable.Due to the absence of a thermally unstable ambient medium, mixing is the primary mechanism for creating the thermally unstable intermediate gas in our simulations.Still, results from our study will be relevant for the late evolution of simulations with thermally unstable hot medium, at which point, the further creation of cold gas is likely dominated by the cooling of mixed intermediate gas, rather than the less unstable hot medium.Importantly, the dynamics of a multiphase medium are quite different depending on which phase dominates the simulation domain.Since in most astrophysical media, the hot component is dominated by volume (see, e.g., Tumlinson et al. 2017, for the CGM), we choose to focus on the initial phase where this is also the case in our setup.Studying the full dynamic range, i.e., having a sufficiently large volume to sustain  V,c ≪ 1 for an extended period of time while resolving the small-scale structure is unfortunately computationally prohibited.
Another similar system of turbulent boxes can be the stratified turbulent boxes, as studied by Mohapatra et al. (2021), Mohapatra et al. (2021) and Wang et al. (2023).In such systems, the fundamental nature of turbulence can be different, depending on the extent of stratification.But, due to the presence of a similar hierarchy of structure and scales, we expect to see a similar growth or destruction of cold gas.In a stratified medium, there are two kinds of motions, one across the stratified layers, i.e. along the stratifying force ( strat ), and the other along the layers, i.e. perpendicular to the  strat .The growth of cold gas within the layer itself would depend on the turbulent property in the layer, roughly perpendicular to  strat , while the transport and growth of cold gas among the stratified layers would depend on the gas motion along  strat .This kind of motion can be turbulent or buoyancy-driven where the cold gas falls "down".A stronger stratification can suppress the turbulent motions across the stratified layers, while the buoyant forces and motions can get amplified.Hence, even though some of our results are relevant to a stratified system, due to the complex interplay between these different flows, further study is needed to fully understand the rich physics in play.
Apart from explicitly turbulent boxes, turbulence shows up time and again in a lot of astrophysical simulations.An example of one such system are the 'cloud-crushing' simulations modelling cold gas-wind interactions.These set of simulations, designed to study multiphase galactic outflows, have been extensively studied (e.g., Klein et al. 1994;Marinacci et al. 2010;Scannapieco & Brüggen 2015;McCourt et al. 2015;Schneider & Robertson 2017;Girichidis et al. 2021).Studies find that cold gas clouds that are bigger than a certain critical radius can not only survive against a fast-moving hot wind but even grow as they are being entrained in the wind (Gronke & Peng Oh 2018;Li et al. 2020) with the details of the critical radius still under debate (Kanjilal et al. 2021;Farber & Gronke 2021;Abruzzo et al. 2022).
Initially, when hit with the hot wind, Kelvin-Helmholtz (KH) rolls are formed near the edges facing perpendicular to the wind, where the relative velocity is the highest.These KH rolls act as one of the initial sources of turbulent motions behind the cloud in its tail and cause mixing.As the cloud gets entrained and the shear decreases, this mechanism is unable to drive any further turbulence.Still, many of the previous studies mentioned above find that the cold gas mass continues to grow even after the cloud is entrained.This points to the presence of a substitute process for driving the turbulence at later times.The nature of this substitute process is still an open question, with some suggestions being the hot gas inflow due to cooling tail (Abruzzo et al. 2022) or the pulsations of the cold clumps themselves (Gronke & Oh (2023), Gronke & Oh (2020b)).Regardless of the exact source of the late-time turbulence driving in the tails, as we show in this study, if the resulting turbulence in the tail is similar, the mixing and the cold gas evolution will be similar.Interestingly, studies with magnetic fields, like Gronke & Oh (2020a) and Hidalgo-Pineda et al. (2023) find a lack of significant difference between the growth rates of the cold gas with (MHD) and without (HD) magnetic fields.This result, in combination with what we find in our study, means that the presence of magnetic fields is not affecting the turbulence-driving mechanism.However, note that Hidalgo-Pineda et al. ( 2023) do find a significant difference in the survival criterion of clouds in a laminar flow with the inclusion of magnetic fields (∼ 2 orders of magnitude with  ∼ 1).To understand this, it is important to recall that the main difference to our turbulent setup is that for a wind tunnel setup the reduction of the drag time ( drag ∼ / cl ∼  1/2  cc ) in order to be comparable to the destruction time  cc is sufficient for survival.Hidalgo-Pineda et al. (2023) attribute this reduction to a combination of draping (Dursi & Pfrommer 2008;McCourt et al. 2015) and an altered  due to compression of magnetic fields.On the other hand, in a turbulent setup, the cold gas is never fully entrained.
Another analogous set of systems is the Ram-pressure stripped galaxies, also called jellyfish galaxies.Similar to the cloud-crushing simulations, such galaxies have a multiphase tail.And, both simulations (Roediger & Brüggen 2006;Tonnesen & Bryan 2009) and observations (Boselli et al. 2022;Li et al. 2022;Luo et al. 2023) have shown the presence of turbulence in the tails of such galaxies.Results from this study will be quite relevant to the environment in such a tail, where the extent of the turbulence in the tail will dictate the overall evolution of the multiphase gas.Even though there are some strong parallels between jellyfish galaxies and cloud-crushing simulations, there are also many differences, like the difference in overdensity, presence of self-gravity, star-formation, feedback, etc.Hence, more detailed studies are required to fully understand these systems.
One of the major sources of turbulence in the circumgalactic medium (CGM) is the galactic outflows caused by the supernova feedback in the galactic disk.In our simulations, we vary the turbulent energy injection rate in order to get a similar turbulent velocity in both HD and MHD simulations.In a more realistic system, as in isolated galaxy simulations, the energy injection is dictated by the supernova rate, and indirectly by the star formation rate (SFR).Previous studies like Hopkins et al. (2019);van de Voort et al. (2021) found that changes to SFR, stellar mass and ISM mass due to the inclusion of magnetic fields are small.This means the energy injection rate into the CGM is roughly unaltered due to the inclusion of magnetic fields.As the magnetic fields in the CGM will act as an additional energy sink, the resulting turbulent velocity in the CGM due to the outflows is expected to be lower when magnetic fields are included.This reduced turbulent velocity in the CGM can be one of the possible reasons for the lower extent of mixing of metals in CGM, resulting in the stronger angular dependence of metallicity in simulations when the magnetic fields are included (van de Voort et al. 2021).
Closer to home, multiphase MHD turbulence is also seen in the solar atmosphere.The nature of MHD turbulence in the solar atmosphere is quite different, due to the very high magnetic field intensities, leading to sub-Alfvénic turbulence.In this case, the magnetic field tension is very high, and magnetic field lines are stiff to the gas flows.Still, as the mixing of multiphase gas is fundamentally tied only to the gas flows, and in the presence of the turbulent cascade of structures, our results suggest that the evolution of the multiphase gas would primarily be affected by the overall turbulent property.One of the sources for this turbulence can be the non-linear evolution of KH instability, which has been investigated in previous studies like Hillier et al. (2023).

Caveats / future directions
Below, we mention some caveats of the study and some directions that can be explored in future studies.
• Resolution: We use a lower resolution in our TRML simulations compared to that in Tan et al. (2021).This should not affect our results because, as Tan et al. (2021) show, it is enough to properly resolve the largest eddy to get converged cooling and mixing rates, which we do.Similarly, Gronke et al. (2022) show that the growth rates and survival of cold gas clouds well within the survival regime, is converged if the cloud radius is well-resolved.As this criterion is satisfied in our simulations, we believe the results should be converged over similar resolutions.A lack of physical resistivity, viscosity or conduction means that in our simulations these are replaced by numerical resistivity, viscosity and conduction.A higher resolution will lead to a decrease in these but, as mentioned in section 5.1, the primary timescale in the problem is the turbulent eddy timescale of the largest eddy, which is unaffected by the resolution.This is similar to the analogous result in TRMLs which Tan et al. (2021) find in their study.
• Turbulent driving: In this study, we maintain a solenoidal to compressive driving ratio (  shear ) of 0.3 across all turbulent box simulations.Previous studies find that different  shear in simulations can cause differences in the turbulent power spectrum (Federrath 2013;Grete et al. 2018;Mohapatra et al. 2022c).But for our results, it is enough that the turbulent eddy timescale of the largest eddy is longer than that of smaller eddies, this remains unchanged with a different turbulent driving.The nature of turbulent driving can also affect the magnetic field amplification in MHD turbulence.The magnetic field in a turbulent box driven by an  shear > 0 is amplified much faster than a purely compressively-driven (  shear = 0) turbulent box.Still, this difference is well within an order of magnitude, and the results from our simulations should largely be applicable to the case of purely compressible turbulence.
• Subsonic vs supersonic: In this study, we restrict ourselves to the subsonic regime in both TRML and turbulent box simulations (to be more applicable to most astrophysical systems).Yang & Ji (2023) have looked at the behaviour of TRMLs with supersonic shear velocities, and find that for very high Mach numbers, the turbulent velocities in the mixing zone start to saturate with increasing shear velocities.This leads to a stagnation in the cooling rate, which is in agreement with our results from TRML simulations.Mohapatra et al. (2022b), in their simulations with supersonic turbulence, find that stronger turbulence can lead to higher compression and rarefactions in the medium.The stronger compression, along with shocks, might cause higher cold gas formation from the cooling of the ambient medium if the cooling is stronger than shock heating.This might also be valid for the supersonic multiphase turbulent boxes with non-cooling ambient medium, analogous to this study, where shocks passing through the medium might result in more efficient cooling of shocked intermediate gas regions.On the other hand, shocks in supersonic turbulence can also lead to the destruction of the cold gas, countering the additional cold gas formation.We see a hint of this more efficient destruction in our transonic (M s ≈ 0.9) turbulent simulations in Fig. 8, where the clouds larger than the subsonic critical radius get destroyed.Hence, the results in an analogous multiphase supersonic turbulent box might vary from the subsonic cases.This is further complicated by the presence of magnetic fields, where there are two kinds of shocks, and these can also lead to the amplification of the magnetic fields.
• Super-Alfvénic vs Sub-Alfvénic: Most of the large-scale astrophysical media like the ISM, CGM and ICM are usually super-Alfvénic (M A >1) in nature.Even though, most of them start with a relatively high M A , due to amplification of the magnetic fields the media reach a lower M A , but usually not equipartition due to temporal evolution, and stay Super-Alfvénic.Similarly, in our simulations, we start well within the Super-Alfvénic regime but during the turbulent driving, we reach equipartition, before we introduce the cold gas cloud.That is, M A ≳ which is tran-Alfvénic to mildly super-Alfvénic.This setup works well to understand the above mentioned astrophysical media, but there are other multiphase environments like the Solar Corona where the medium is well within the sub-Alfvénic regime and our turbulent boxes may not be analogous anymore.On the other hand, our TRML simulations include simulations with trans-Alfvénic to mildly sub-Alfvénic motions.We find that our conclusion about the relation between turbulent velocity and mixing still holds.This means, given there are turbulent motions and a turbulent cascade, the mixing will only depend on the turbulent properties and not the presence or absence of the magnetic fields.Still, we have not explicitly tested this in a turbulent box setup but can be a topic for future investigations.
• Anisotropic conduction: It is well-known that conduction is anisotropic in the presence of magnetic fields.But, as we do not have physical conduction in our simulations, the numerical conduction in the simulations is isotropic in both HD and MHD cases.While it has been shown by Tan et al. (2021); Tan & Oh (2021) that generally turbulent diffusion dominates over the laminar one (thus, explaining seemingly 'puzzling' convergence of larger scale studies such as ours), this has only recently been investigated with anisotropic conduction in an MHD setup by Zhao & Bai (2023) who corroborate our results and find similar trend in suppression of cooling (see, however, Brüggen et al. 2023;Jennings et al. 2023, who included anisotropic conduction in their 'cloud-crushing' simulations and find similar mass growth rates as the pure hydro runs).This will add an additional layer of complexity, and can also be a future direction to explore.
• Other effects neglected in this study are cosmic rays, viscosity, and geometrical variations such as stratification.Our goal here was to study mixing in MHD in a simplified setup to which we will add additional layers of complexity in future work.
Our study reconciles the seemingly contradictory results of the effect of magnetic fields in turbulent mixing layers and a fully multiphase turbulent setup.This result also implies that the presence of cold gas in multiphase media can be explained through continuous mixing and cooling -and this channel is not hindered by the presence of magnetic fields.However, the topic of multiphase MHD turbulence still remains full of many unanswered questions, like the effect of cosmic rays, thermal conduction, viscosity, etc. which we hope to tackle in future work.2021), and the right panel shows the Gaussian filtering method used in Abruzzo et al. (2022).We find only minor differences between the two methods which, at worst, stay within an order of magnitude.(ii) Construct the neighbourhood graphs from all the adjacency matrices constructed in the previous step.
(iii) Calculate the shortest path between each node in the smallest of the neighbourhood graphs created in the previous step.
(iv) Find the longest of the set of calculated shortest paths and note the nodes corresponding to that path.
(v) Recalculate the length of the longest "shortest" path between the nodes from the previous step, using the largest neighbourhood graph.
(vi) The length from the previous step gives a rough measure of the filament length in the clump.Repeat the steps for all clumps.
We tested the above method for different numbers of skipped points for calculating the shortest paths.We find a negligible difference in the calculated length of large clumps even up to the point where every 20 th point is considered.We see major deviations only when the skipped points are a big majority of the points and the resulting neighbourhood graph is not representative of the clump anymore.
In this work, we only skip every 4 th point in the clump.Fig. B1 shows an example of the calculated filament length for a clump in an MHD turbulent box simulation.

APPENDIX C: MOCK SPECTRA
As mentioned in § 4.5 and 5.3, we find only marginal differences between the statistics of the MgII mock absorption spectra of the HD and MHD simulations, despite significant differences in the column densities.This can be due to the specific property of the MgII 2796Å absorption line, like the curve of growth flattening around similar column densities, which can lead to smaller differences.Another possible reason for this lack of difference can also be the spectral resolution.To address both these points, we first increase the spectral resolution of the mock spectra tenfolds to Δ = 0.01Å and recreate the same MgII 2796Å mock absorption spectra analysis as Fig. 18.Secondly, we repeat the same analysis for CIV 1551Å at the higher We find that the increase in spectral resolution of mock MgII absorption spectra shifts the relation between the number of features and total equivalent width, but it roughly follows the same slope as the observed relation from Churchill et al. (2020).Surprisingly, the statistics of the mock CIV 1551Å absorption spectra also seem to agree with the observed MgII relation, and the HD-MHD differences are wider as expected, but to the lower number of unsaturated mock spectra, it is harder to draw concrete conclusions.This apparent robustness of the observed relation might hint towards a more fundamental origin of the relation, like the clump distribution.But, we leave it to future studies to investigate this further.

APPENDIX D: EFFECT OF STOCHASTICITY
The stochastic nature of the turbulence can cause variations in the evolution of quantities in a turbulent environment.Gronke et al. (2022) found that this stochasticity affects the cold gas mass evolution in hydrodynamic turbulent boxes with an intermediate-sized initial cold cloud.In this regime, they saw both survival and destruction of the cold cloud for different choices of random seeds for turbulent driving.This was attributed to the higher significance of the exact turbulent velocity field in the intermediate regime between cloud survival and destruction.
We repeat this test for our simulations, with and without magnetic fields.We run turbulent box simulations with 3 different random seeds at M = 0.5 and introduce clouds of different sizes to check for the effect of stochasticity of the turbulence.We use a  box / cl = 20, instead of 40, due to its lower computational costs.
Fig. D1 shows the cold gas mass evolution for the different cases.We find that the cold gas mass growth/destruction rate for cold gas clouds in intermediate and destruction regimes is sensitive to the exact choice of the random seed.We also find that this is true for both HD and MHD and with no clear order of growth rate between the HD and MHD counterparts.This high dependence on stochasticity in these regimes is due to the lack of cold gas mass.This results in a very stochastic sampling of turbulence, hence making the evolution very stochastic in nature.
On the other hand, in the survival regime, the MHD simulations seem to have a slightly lower growth rate, compared to their HD counterparts, although still a much lower difference compared to the order of magnitude difference observed in TRML simulations.We attribute this minor difference to some unavoidable systematic differences between the HD and MHD simulations.The biggest of them is the difference in dissipation rate between MHD and HD, due to the extra dissipation of magnetic energy via numerical resistivity.This higher dissipation results in a slightly hotter medium in a fully developed turbulent box, in turn resulting in a slightly deviated density distribution.These slight deviations affect the evolution via a slight difference in overdensity, mixed gas temperature, etc.
Still, as Fig. D1 shows, this difference is minor and it gets even more trivial when we take the spread due to stochasticity into account.

Figure 1 .
Figure 1.Temperature slices for different TRML simulations for  shear = 100 km s −1 (M ≈ 0.3).First column shows the hydrodynamic simulations, 2 nd to 4 th column show simulations with M A = 10, last three columns show simulations with M A = 1.Top row shows simulations with strong cooling, Da = 60, Bottom row shows simulations with weak cooling, Da = 6 × 10 −5.This shows the different ways magnetic fields evolve for different initial orientations.It also suggests that the cases with the higher magnetic field strength along the shear direction show a lesser extent of mixing.

Figure 2 .
Figure 2.  =  thermal / magnetic slices for different TRML simulations with  shear = 100km/s (corresponding to the temperature slices shown in Fig. 1). 1 st to 3 rd column show simulations with M A = 10, last three columns show simulations with M A = 1.Top row shows simulations with strong cooling,  = 60, Bottom row shows simulations with weak cooling, Da = 6 × 10 −5.This shows the extent of amplification possible in the different cases.Even with a higher initial , turbulent motions can amplify the magnetic fields to lower .For cases with lower initial  strong cooling in the mixing layer can also lead to amplification.

Figure 3 .
Figure 3. Top & middle Stability criterion of Kelvin-Helmholtz instability for different initial magnetic field orientations (cf.Eq. (8)), Bottom Density profile for the different cases shown above.This shows the difference in the stability of the mixing layers for different cases.The cases with higher magnetic field strength along the shear (initial or amplified) are more stable.The density profile shows the extent of the mixing layers.

Figure 4 .
Figure 4. Left column M A = 1, Right column M A = 10, Top row Stable values of mixing layer surface brightness for different Da, The orange and blue dashed lines on the top row panels are the expected values from hydrodynamic TRML simulations byTan et al. (2021).This shows clear suppression in cooling rates of most of the simulations with magnetic fields, in comparison to hydrodynamic TRML simulations.We discuss details about the trends in § 3. Bottom row Evolution of mixing layer surface brightness with time, for different initial magnetic field orientations at two Da values.

Figure 5 .
Figure5. ′ profiles along ẑ for different initial magnetic field orientations.The first two panels from the left show the difference between the direction parallel to magnetic fields versus the other directions.The third panel shows the outlier case of magnetic fields normal to the interface, where both the normal ( ẑ) and shear direction show much higher fluctuations due to the presence of magnetic fields along both these directions.Hence, we choose the directions which free from these spurious fluctuations in these different cases, as denoted in Eq 13

Figure 6 .
Figure 6.Top panel Scatter plot of the surface brightness () and turbulent velocity ( ′ ) calculated from the simulations.The dashed and dotted lines show the respective strong and weak cooling scaling relation according to Eq. 9-10.Bottom panel Similar to the top panel, after we remove the Da dependence.The dotted line shows the analytical expectation from Tan et al. (2021), which they find for hydrodynamic simulations.This suggests that the general relation found in hydrodynamic TRML simulations, between the turbulent velocity in the mixing layer and cooling (and hence mixing) rate, is still valid in the presence of magnetic fields.

Figure 8 .
Figure 8. Survival or destruction of the cold gas in the different turbulent boxes.The dashed line is the survival criterion fromGronke et al. (2022).This shows the surprising lack of difference between the survival criterion, with and without magnetic fields.The subsonic turbulent simulations agree well with the previously found survival criterion, with some deviation in trans-sonic turbulent boxes (c.f.§ 4).

Figure 9 .
Figure 9. Cold gas evolution with time for different simulations initiated with varying sizes of cold gas cloud in turbulence with M = 0.5.Solid lines show the simulations with magnetic fields, dashed lines show the hydrodynamic simulations and the dotted lines show the expected hydrodynamic growth rates from Gronke et al. (2022).This shows that there are only marginal differences between the growth and destruction rates of the cases with and without magnetic fields, compared to the differences seen in the TRML simulations.The differences further diminish as we consider cases well within the survival regime.

Figure 10 .
Figure 10.Cumulative number distribution for HD-MHD simulation pair with M = 0.5 and  cloud = 310 shatter .This shows the marginal difference in the overall distribution of clump sizes, and also that the distribution matches the distribution of ∝  −1 , found in previous studies.

Figure 11 .
Figure11.Histogram of longest shortest distance in the neighbourhood graph of every clump in a snapshot from the turbulent box at M turb = 0.5, with and without magnetic fields.This figure gives a lower limit on the difference in the filamentariness of the cold gas clumps in the two cases.We find at cold gas clumps can get more filamentary in the presence of magnetic fields, by about a factor of 2.

Figure 13 .
Figure13.Evolution of average shear at clump boundaries in a set of HD and MHD simulations with M = 0.5 and  cl = 310 shatter (same as Fig.12).The shaded regions show the corresponding 15-85%ile intervals.The figure also shows the shear on the clump boundaries is about an order of magnitude lower than the turbulence velocity in the simulations.Also, on average, clumps in the MHD simulation seem to have a marginally lower, but very similar shear, in comparison to HD simulations.
Figure15.Top Average, median and 10 average of logarithm of entanglement, i.e.  stream / for different streamline lengths ( stream ).The dashed lines show the corresponding best linear fits and the shaded region shows the 15-85%ile interval.The general trend of increasing entanglement for longer and longer streamline lengths indicate a fractal-like structure of the magnetic field lines, discussed further in § 5.2.Bottom inset Points denote the mean and variance of log 10 ( stream /) and the green dashed line shows the linear fit,  = 0.24(  − 0.03).We use this relation to calculate the shown probability distribution.Bottom Solid lines show the probability distributions of different values of entanglement, log 10 ( stream /), for three values of streamline lengths.The dashed lines show the corresponding calculated Γ distributions, with the parameters mentioned in the legend.This shows the close agreement between the estimated and calculated probability distributions.There are some deviations for the probability distribution of small streamline length, which is discussed further in § 4.4.

Figure 16 .
Figure 16.Column density distribution of cold ( < 10 5 K, left panel) and intermediate (10 5 K <  < 10 6 K, right panel) temperature gas in HD (in green) and MHD (in blue) simulations, with M = 0.5 and  cl = 310 shatter .This shows that the column densities for the above cases are within the observationally expected column densities for absorption spectra in a circumgalactic environment.It also shows that the lower end of column density distribution for cold temperature gas has a higher extent of difference between the HD and MHD simulations.This makes an absorption line tracing the cold gas a prime candidate for looking at observational differences between the HD and MHD simulations.

Figure 17 .
Figure 17.An example line-of-sight MgII 2796 Å absorption mock spectra with Δ = 0.01Å, from the HD (blue solid line) and MHD (green dashed line) simulations with M = 0.5 and  cl = 310 shatter (same as Fig.16).The dotted black line shows the threshold of the minimum absorbed flux of a feature, and the red circles show the features that we consider for analysis.This figure is only for reference, as these are higher resolution spectra compared to the ones used in the analysis at Δ = 0.1Å, which is closer to observational spectral resolution.

Figure 18 .
Figure 18.Contour plot of the 2D histogram of line-of-sight MgII absorption mock spectra in the number of absorption features vs. equivalent width space,for HD (solid contours) and MHD (dashed contours) simulations with M = 0.5 and  cl = 310 shatter (same as Fig.16 and 17).The dash-dotted green line shows the relation found inChurchill et al. (2020).This shows that there are only marginal differences in the overall distributions of HD and MHD simulations, despite the differences in Fig.16.We also find that they agree quite well with the observed relations fromChurchill et al. (2020).

Figure A1 .
Figure A1. ′ profiles for two different methods.The left panel shows the averaging method used in Tan et al. (2021), and the right panel shows the Gaussian filtering method used in Abruzzo et al. (2022).We find only minor differences between the two methods which, at worst, stay within an order of magnitude.

Figure A2 .
Figure A2.Same figure as Fig 6, but the  ′ is calculated using the Gaussian filtering method from Abruzzo et al. (2022).This shows that the results in Fig 6 are not sensitive to the method used to calculate the turbulent velocity.

Figure B1 .
Figure B1.Red points show the points from the skipped graph of a single clump, and the blue solid line shows the calculated filament length using every 4 th point.Axis labels correspond to the number of gridcells.

Figure C2 .
Figure C2.Same as Fig. 18 but for high resolution CIV Å absorption spectra with  = 0.01Å

Figure D1 .
Figure D1.Cold gas mass evolution for simulations with the same parameters but different random instances of turbulence.The different panels refer to different  cl / shatter in a turbulent medium with M = 0.5.The different colours denote simulations with varying random seeds for turbulent driving.The solid and dashed lines show the evolution of simulations with and without magnetic fields, respectively.