© The Authors 2026

1. Introduction

Turbulence and multiphase systems are ubiquitous in plenty of astrophysical environments, including the interstellar medium (ISM; McKee & Ostriker 1977; Audit & Hennebelle 2010), circumgalactic medium (CGM; Tumlinson et al. 2017; Faucher-Giguère & Oh 2023), galactic winds (Veilleux et al. 2020; Thompson & Heckman 2024), and intracluster medium (ICM; Markevitch & Vikhlinin 2007; Kravtsov & Borgani 2012). Therefore, a proper modeling of turbulence and multiphase gas is fundamental to understanding key processes in astrophysics.

Observations provide clear evidence of multiphase gas in the CGM, as well as in other astrophysical environments. For instance, hot gas (T ≳ 10⁶ K) in halos around galaxies has been detected in X-rays (e.g., Anderson & Bregman 2010; Bregman et al. 2018). Meanwhile, cool gas (T ∼ 10⁴ K) has been detected through Lyα emission (e.g., Steidel et al. 2010; Wisotzki et al. 2018; Arrigoni Battaia et al. 2023), as well as absorption lines in the spectra of background quasars (e.g., Tripp et al. 1998; Péroux et al. 2019).

Turbulence is inherent to astrophysics, influencing the gas dynamics across very diverse environments. It plays a crucial role in distributing the energy from large to small scales in a turbulent cascade. The energy is injected at large scales and it cascades down until it reaches the Kolmogorov scale, where the kinetic energy is dissipated into heat (Kolmogorov 1941, 1962). Turbulence is key across many different scales in astrophysics: from small scales such as stellar interiors (Brüggen & Hillebrandt 2001) and protoplanetary disks (Johansen et al. 2006) to large scales such as active galactic nucleus jets (Scannapieco & Brüggen 2008) and cold fronts (ZuHone et al. 2010).

Turbulence mixes the hot and cold gas, producing intermediate temperature gas. In nonideal gases, this intermediate gas is unstable and determines the amount of cold gas available for star formation. This mass transfer has been extensively studied in the context of turbulent radiative mixing layers (TRMLs; e.g., Begelman & Fabian 1990; Slavin et al. 1993), which can be applied to several astrophysical systems, such as galactic outflows (Brüggen & Scannapieco 2016; Schneider & Robertson 2017; Fielding & Bryan 2022) or stellar wind bubbles (El-Badry et al. 2019; Lancaster et al. 2021a,b). Based on combustion theory, Tan et al. (2021) identified a key relation between the turbulent and cooling timescales, expressed by the Damköhler number: Da = t_turb/t_cool (Damköhler 1940). This dimensionless number characterizes the cooling strength and describes a different behavior of cooling efficiency, depending on whether the system is in the weak or the strong cooling regime. Roughly speaking, when Da < 1, cooling is slow compared to turbulence and mixing dominates the gas entropy. When Da > 1, cooling is faster than turbulence and the mixed gas is cooled down rapidly, leading to a multiphase system.

The total bolometric luminosity of a system depends on how much gas can be mixed to intermediate temperatures; therefore, the suppression of turbulence and mixing might lead to differences in luminosity. This has been studied in the context of TRMLs with magnetic fields, where the magnetic field orientation is key for the suppression of the Kelvin-Helmholtz instability (KHI; Chandrasekhar 1961). Ji et al. (2019) found that regardless of the initial orientation, magnetic fields are amplified by the shear motions, leading to the suppression of mixing and, therefore, less cooling, in the slow cooling regime. However, Das & Gronke (2023) found that although magnetic fields suppress mixing, the relation between turbulent velocity and cooling rate found by Tan et al. (2021) in the fast cooling regime still holds in the magnetized case.

Transport processes such as viscosity also affect the development of instabilities and mixing. Even though they are non-magnetized, fully ionized plasmas are expected to be very viscous (Spitzer 1962; Zeldovich & Raizer 1967) and the effective viscosity is reduced in the presence of magnetic fields (Braginskii 1965; Squire et al. 2023) and plasma microinstabilities (Schekochihin & Cowley 2006; Kunz et al. 2014). Therefore, the overall effect of viscosity in astrophysical systems is still under debate (Zhuravleva et al. 2019; Heinrich et al. 2024; Marin-Gilabert et al. 2024). Although the net effective viscosity in astrophysical systems is not clear yet, it definitely influences the development of instabilities and turbulence (Roediger et al. 2013; Marin-Gilabert et al. 2022). In a similar way to magnetic fields, viscosity changes the mixing efficiency in multiphase systems. In this paper, we aim to study in detail these differences in turbulence. In particular, we consider whether this might impact the universal relation between turbulent velocity and cooling rate for inviscid TRMLs (Fielding et al. 2020; Tan et al. 2021).

Viscosity sets the Kolmogorov scale (η_L; Kolmogorov 1941), namely, the scale at which the kinetic energy is transformed into heat in the turbulent cascade. Below this scale, hydrodynamical instabilities are suppressed. On the other hand, conduction sets the Field length (δ_L; Field 1965), the maximum scale over which conductive energy transport is effective in the presence of radiative cooling (Begelman & McKee 1990). Thermal instabilities are suppressed on scales smaller than the Field length. An important dimensionless parameter in combustion theory is the Karlovitz number (Klimov 1963; Williams 1975), which compares these two scales, Ka = δ_L²/η_L². The Karlovitz number relates the flame thickness to the turbulent eddies scale, determining the shape of the flame front. If Ka < 1, the small-scale turbulence is not small enough to break up the flame and the combustion happens in a laminar flow. If Ka > 1, the small-scale turbulence grows faster than it can locally burn through its thickness, leading a turbulent flame front (Kuo & Acharya 2012). Translated to TRMLs terminology, Ka < 1 describes a scenario in which cooling takes place in a laminar regime, whereas if Ka > 1, cooling happens in a turbulent regime. In the absence of explicit conduction and viscosity, the Kolmogorov scale and the scale of numerical heat diffusion are given by a similar multiple of the resolution scale in both cases (η_L ∼ δ_L ∼ Δx) and, thus, Ka ∼ 1 (Tan et al. 2021). In this paper, we study the parameter space where Ka < 1 by explicitly including viscosity, while relying on numerical conduction, thus ensuring η_L > δ_L.

The aim of this study was to understand how viscosity affects fluid dynamics in idealized setups, for both incompressible and compressible regimes. The findings of these control cases are then applied to TRMLs, where we study the interplay between viscosity and cooling. By including an explicit viscosity, we studied the parameter space in which Ka < 1 and the implications that this might have for turbulence and mixing, as well as for cooling efficiency. In particular, we investigate whether the universal relations found for inviscid fluids (Fielding et al. 2020; Tan et al. 2021) still hold. These results have potential implications in different astrophysical systems. For instance, they might influence the star formation rate in jellyfish galaxies or the gas dynamics in galactic outflows and the CGM.

This paper is organized as follows. In Sect. 2 we introduce the relevant theory. The different numerical setups we use are explained in Sect. 3. In Sects. 4.1 and 4.2 we show the results from idealized 2D simulations, while in Sect. 4.3 we show the results for TRMLs. We expand our discussion in Sect. 5 and present our conclusions in Sect. 6.

2. Theoretical considerations

The fluid dynamics are described by the equations of hydrodynamics: conservation of mass, momentum and energy. Including viscosity, the equations are given by

$$\begin{array}{c}\hfill \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·\left(\rho \mathit{v}\right)=0,\end{array}$$(1)

$$\begin{array}{c}\hfill \frac{\partial (\rho \mathit{v})}{\partial t}+\mathrm{\nabla }·(\rho \mathit{v}\mathit{v}+P)=-\mathrm{\nabla }·\mathit{\sigma },\end{array}$$(2)

$$\begin{array}{c}\hfill \frac{\partial (\rho e)}{\partial t}+\mathrm{\nabla }·\left[(\rho e+P)\mathit{v}\right]=-\mathrm{\nabla }·(\mathit{v}·\mathit{\sigma }).\end{array}$$(3)

ρ is the gas density, v is the velocity of the gas, P is the thermal pressure given by

$$\begin{array}{c}\hfill P=(\gamma -1)\rho e_{\mathrm{int}}\end{array}$$(4)

and e is the total energy per unit mass

$$\begin{array}{c}\hfill e=e_{\mathrm{int}}+\frac{1}{2}{\mathit{v}}^2,\end{array}$$(5)

where e_int is the internal energy per unit mass. Finally, σ is the isotropic viscous stress tensor

$$\begin{array}{c}\hfill \mathit{\sigma }=\eta (\mathrm{\nabla }\mathit{v}+\mathrm{\nabla }{\mathit{v}}^{\mathrm{T}}-\frac{2}{3}\mathrm{\nabla }·\mathit{v}),\end{array}$$(6)

where we have omitted the bulk viscosity term, which is zero for monoatomic gases, since it is related to the degree of freedom of molecular rotations (see, e.g., Zeldovich & Raizer 1967; Pitaevskii & Lifshitz 1981). Then, η is the shear viscosity coefficient, which for fully ionized (non-magnetized) plasmas depends on the temperature of the gas as η ∝ T^5/2 (Spitzer 1962; Braginskii 1965). For partially ionized plasmas, the effect of viscosity is reduced and depends on the ionized fraction (Su et al. 2017; Hopkins et al. 2019).

2.1. Critical viscosity theoretical estimate

For reference, our numerical simulations will be computed as a function of the critical viscosity needed to suppress the growth of the KHI. Due to the nonlinear behavior of the KHI, there is no analytical solution to describe its evolution, so we need to estimate it based on a different set of assumptions. The theoretical estimate we used in this work comes from Roediger et al. (2013) and later used in Marin-Gilabert et al. (2022).

Viscosity smooths out the shear velocity gradient over a length scale ±d above and below the interface due to momentum diffusion. This leads to the suppression of the KHI and, thus, turbulence and mixing. The instability is suppressed if the wavelength of the perturbation (λ) is smaller than ∼10 d (Chandrasekhar 1961).

We assume that the instability will be fully suppressed when the viscous timescale, τ_ν = d²/ν, at which the velocity gradient is smoothed out a distance d = λ/10 is shorter than the Kelvin-Helmholtz timescale¹,

$$\begin{array}{c}\hfill {\tau }_{\mathrm{KH}}=\frac{\lambda }{\mathrm{\Delta }{v}_{\mathrm{shear}}}\frac{({\rho }_{\mathrm{hot}}+{\rho }_{\mathrm{cold}})}{{({\rho }_{\mathrm{hot}}{\rho }_{\mathrm{cold}})}^{1/2}}.\end{array}$$(7)

Here, ρ_cold and ρ_hot are the densities of the cold and hot medium respectively, Δv_shear is the velocity difference between the two media, and ν is the kinematic viscosity, which depends on the shear viscosity coefficient via

$$\begin{array}{c}\hfill \nu =\eta /\rho .\end{array}$$(8)

The instability becomes stable and, therefore, fully suppressed when

$$\begin{array}{c}\hfill {\tau }_{\nu }<{\tau }_{\mathrm{KH}},\end{array}$$(9)

which can be rewritten to

$$\begin{array}{c}\hfill \nu >{\nu }_{\mathrm{Crit}}=\frac{\lambda \mathrm{\Delta }{v}_{\mathrm{shear}}}{100}\frac{{({\rho }_{\mathrm{hot}}{\rho }_{\mathrm{cold}})}^{1/2}}{({\rho }_{\mathrm{hot}}+{\rho }_{\mathrm{cold}})}.\end{array}$$(10)

2.2. Radiative mixing layers

When the KHI develops, it quickly triggers turbulence, mixing the cold and hot medium into intermediate temperature gas, entering the regime where the cooling curve peaks (Begelman & Fabian 1990; Fielding et al. 2020; Tan et al. 2021). To characterize the evolution it is convenient to introduce the Damköhler number (Da; Damköhler 1940), a dimensionless number that accounts for the cooling strength:

$$\begin{array}{c}\hfill \mathrm{Da}=\frac{t_{\mathrm{turb}}}{t_{\mathrm{cool}}(T)}=\frac{\lambda }{u\prime t_{\mathrm{cool}}(T)},\end{array}$$(11)

where t_turb is the eddy turnover time, t_cool the cooling time, and u′ the turbulent velocity.

When t_turb < t_cool (Da < 1), the effect of cooling is weak compared to turbulence, and the gas entropy is dominated by mixing. This leads to a characteristic mass flux $\propto t_{\mathrm{cool}}^{-1/2}$ (Ji et al. 2019; Tan et al. 2021). However, when t_cool < t_turb (Da > 1), cooling is stronger than turbulence and, as soon as the gas is mixed, it is cooled down, leading to a multiphase medium where the gas entropy is dominated by cooling (see, e.g., Kuo & Acharya 2012; Tan et al. 2021). Specifically, in this strong cooling regime, the total luminosity of the mixing layer scales as $\propto t_{\mathrm{cool}}^{-1/4}$ (Gronke & Oh 2019; Fielding et al. 2020; Tan et al. 2021).

In summary, the emission due to cooling follows a power-law relation, which depends on the cooling strength (Damköhler number). The surface brightness Q depends on Da via

$$\begin{array}{c}\hfill Q\propto \{\begin{array}{cc}u\prime {}^{1/2}{\lambda }^{1/2}\propto u\prime {\mathrm{Da}}^{1/2},\hfill & \text{if}\mathrm{Da}<1\text{(weak cooling);}\hfill \\ u\prime {}^{3/4}{\lambda }^{1/4}\propto u\prime {\mathrm{Da}}^{1/4},\hfill & \text{if}\mathrm{Da}>1\text{(strong cooling)}.\hfill \end{array}\end{array}$$(12)

3. Numerical setup

We performed our simulations using the code ATHENA (Stone et al. 2008). We used the HLLC Riemann solver, a second-order reconstruction with slope limiters in the primitive variables, the van Leer unsplit integrator (Gardiner & Stone 2008), Cartesian geometry, and an adiabatic equation of state with γ = 5/3. In all our simulations, we kept the shear viscosity coefficient (η) constant to simplify the analysis and understand the effect of viscosity in the simplest way².

3.1. Adiabatic 2D simulations

To better understand fluid dynamics and the effects of viscosity while reducing computational costs, we used two different 2D setups without cooling. These provided a more controlled environment compared to a chaotic 3D system, while allowing us to study the KHI in detail. The first setups explore how critical viscosity depends on overdensity and Mach number by varying viscosity across different overdensity and Mach number values. We describe each 2D setup in detail below.

3.1.1. The planar slab

To understand how the critical viscosity at which the KHI is suppressed changes as a function of the overdensity, we used a modified version of the numerical setup suggested by Lecoanet et al. (2015), which provides a benchmark for the growth of the KHI in grid codes. We created a 2D domain with a resolution of 320 × 320 in the $\widehat{x}$ and $\widehat{y}$ directions, respectively, with periodic boundary conditions with dimensions L_x = 256 and L_y = 256. The density and velocity profiles satisfy

$$\begin{array}{c}\hfill \rho ={\rho }_{\mathrm{hot}}+\frac{{\rho }_{\mathrm{cold}}-{\rho }_{\mathrm{hot}}}{2}[tanh\left(\frac{y+y_{\mathrm{Int}}}{a}\right)-tanh\left(\frac{y-y_{\mathrm{Int}}}{a}\right)],\end{array}$$(13)

$$\begin{array}{c}\hfill v_x=-v_{\mathrm{shear}}\times [tanh\left(\frac{y+y_{\mathrm{Int}}}{a}\right)-tanh\left(\frac{y-y_{\mathrm{Int}}}{a}\right)-1].\end{array}$$(14)

Here, ρ_hot is a fixed value for all our simulations and we end up modifying ρ_cold as we change the value of overdensity. Then, y_Int = 64 indicates the position of the interface between the two fluids and v_shear is the shear velocity of each fluid; namely, the relative motion can be expressed as ℳ ≡ 2v_shear/c_c, hot.

In this setup, we kept ℳ = 0.34 constant, as well as T_cold = 10⁴ K, while T_hot was given by the overdensity to ensure pressure equilibrium in the whole domain. Overall, the value of a = λ/10 ensures smooth initial conditions, fundamental for the growth of the desired instability mode without spurious secondary instabilities in grid codes; λ = 128 indicates the wavelength of the perturbation used to trigger the instability. We used the same initial perturbation as in Marin-Gilabert et al. (2022), a modified version of the one introduced by Read et al. (2010), expressed as

$$\begin{array}{cc}\hfill v_y& =-\delta v_y[sin\left(\frac{2\pi (x+\lambda /2)}{\lambda }\right)exp(-{\left(\frac{y-y_{\mathrm{Int}}}{\sigma }\right)}^2)\hfill \\ \hfill & +sin\left(\frac{2\pi x}{\lambda }\right)exp(-{\left(\frac{y+y_{\mathrm{Int}}}{\sigma }\right)}^2)],\hfill \end{array}$$(15)

where δv_y = v_shear/10 is the initial amplitude of the perturbation and σ = 0.2λ is a scaling parameter to control the width of the perturbation layer.

To study how the critical viscosity changes with the overdensity of the system, we used different values of χ ≡ ρ_cold/ρ_hot from 1 to 100 and tested different values of viscosity for each overdensity. For reference, the Spitzer viscosity of a fluid with T = 10⁶ K is η_Sp ≃ 0.11g(cm s)⁻¹, which corresponds to η_Sp ≃ 1809.28 in our unit system.

3.1.2. The planar sheet

When increasing the Mach number, body modes caused by pressure waves moving across the domain between the stream boundaries appear in the planar slab setup (Mandelker et al. 2019). To avoid this, we modified our initial conditions to study the dependence of the critical viscosity as a function of the Mach number. We used a similar setup as before, but with just one interface between the two fluids at y_Int = 0 and open boundaries in the $\widehat{y}$ direction. The density and velocity profiles are expressed as

$$\begin{array}{c}\hfill \rho =\frac{1}{2}\times [{\rho }_{\mathrm{cold}}(1-tanh\left(\frac{y}{a}\right))+{\rho }_{\mathrm{hot}}(1+tanh\left(\frac{y}{a}\right))],\end{array}$$(16)

$$\begin{array}{c}\hfill v_x=v_{\mathrm{shear}}\times tanh\left(\frac{y}{a}\right),\end{array}$$(17)

and the initial perturbation is

$$\begin{array}{c}\hfill v_y=-\delta v_{y}sin\left(\frac{2\pi (x+\lambda /2)}{\lambda }\right)exp(-{\left(\frac{y}{\sigma }\right)}^2).\end{array}$$(18)

With this setup, we want to study the dependence of the critical viscosity on the Mach number. Therefore, we kept χ = 10 constant and we modified v_shear to obtain Mach values ranging from 0.01 to 1.8. As with the previous setup, we tested different values of viscosity for each Mach number to identify the critical viscosity in each case.

3.2. 3D simulations: Turbulent radiative mixing layers

Once we obtained the value for the critical viscosity in a system with an overdensity of 100, we made use of this value to study how different values of viscosity affect a TRML system. For a proper development of turbulence, a 3D setup is needed; therefore, we used a modified version of the setup suggested in Ji et al. (2019), where the density profile is the same as in Eq. (16), but the velocity profile is

$$\begin{array}{c}\hfill v_x=\frac{v_{\mathrm{shear}}}{2}\times [1-tanh(-\frac{y}{a})],\end{array}$$(19)

while the initial perturbation is

$$\begin{array}{c}\hfill v_y=\delta v_{y}sin\left(\frac{2\pi x}{\lambda }\right)sin\left(\frac{2\pi z}{\lambda }\right)exp(-{\left(\frac{y}{a}\right)}^2).\end{array}$$(20)

In this case, v_shear corresponds to the velocity difference between the media, where the lower (colder) medium is stationary and the upper one moves at v_shear. We used a Mach number of ℳ = 0.5 and δv_y = v_shear/100. In this case our domain is ten times larger in the $\widehat{y}$ direction, whereas in the $\widehat{x}$ and $\widehat{z}$ direction is equal to the wavelength of the perturbation L_y = 10L_x = 10L_z = 10λ. We kept the hot gas density and the cold gas temperature constant, ρ_hot and T_cold = 10⁴ K, respectively. In this setup, we use a fixed overdensity χ = 100 (i.e., T_hot = 10⁶ K). The resolution is the same as in Das & Gronke (2023), consisting of 64 × 640 × 64 in the $\widehat{x}$, $\widehat{y}$, and $\widehat{z}$ directions respectively, with outflow boundaries in the $\widehat{y}$ direction and periodic boundaries in $\widehat{x}$ and $\widehat{z}$.

We used the Townsend radiative cooling algorithm (Townsend 2009, as, e.g., in Gronke et al. 2021), with a solar metallicity cooling curve (Sutherland & Dopita 1993) and a T_floor = 10⁴ K. In contrast to Tan et al. (2021), we did not modify the cooling strength. Instead, we modified the wavelength to have different values of Damköhler number. Since Da = t_turb/t_cool, and t_turb ∝ λ (see Eq. (11)), by varying the wavelength (and also the size of the box), we got different values of Da. We ran a total 12 simulations with different Da for each amount of viscosity, from Da_mix ≃ 2.3 × 10⁻³ to Da_mix ≃ 1.1 × 10³, where the temperature of the intermediate mixed gas chosen to compute the cooling time is T_mix = 2 × 10⁵ K (to be consistent with Tan et al. 2021). The values of Da were calculated assuming that t_turb ∼ τ_KH. We know that in the presence of cooling, the box is filled with cold gas, causing the mixing layer leaves the domain, especially in the strong cooling regime. To prevent this, we can calculate the velocity at which the box is being filled with cold gas, adding a velocity in the opposite direction to the whole box, following Das & Gronke (2023).

4. Results

4.1. Critical viscosity versus overdensity

We provide an analytical estimate for the critical viscosity in a KHI setup in Sect. 2.1 (see Eq. (10), based on the growth time of the KHI and the viscous time. Numerically, this critical viscosity is computed by measuring the growth or decay of the instability within 1τ_KH (see Appendix A for details).

The analytical estimate presented in Sect. 2.1 sets a criterion for the kinematic viscosity (ν). However, in our simulations, we used a constant dynamic viscosity (η) for both fluids. In a multi-fluid system with different densities but the same η, each fluid will exhibit its own ν. Therefore, to be able to compute the critical η for our setup, it is necessary to define an appropriate average of these ν to find a general $\overline{\nu }$ of the whole system. Thus, for a system characterized by a single dynamic viscosity (η), two densities (ρ_h and ρ_c), and consequently two different kinematic viscosities (ν_h and ν_c), we define an averaged kinematic viscosity using the harmonic mean, namely,

$$\begin{array}{c}\hfill {\nu }_{\mathrm{h}}=\frac{\eta }{{\rho }_{\mathrm{h}}},{\nu }_{\mathrm{c}}=\frac{\eta }{{\rho }_{\mathrm{c}}},\end{array}$$(21)

$$\begin{array}{c}\hfill \overline{\nu }=\frac{{\sum }_i{\eta }_i}{{\sum }_i{\rho }_i}=\frac{2\eta }{{\rho }_{\mathrm{h}}+{\rho }_{\mathrm{c}}}=\frac{2\eta }{\frac{\eta }{{\nu }_{\mathrm{h}}}+\frac{\eta }{{\nu }_{\mathrm{c}}}}=\frac{2}{\frac{1}{{\nu }_{\mathrm{h}}}+\frac{1}{{\nu }_{\mathrm{c}}}},\end{array}$$(22)

which is equivalent to a density weighted mean,

$$\begin{array}{c}\hfill \overline{\nu }=\frac{{\nu }_{\mathrm{h}}{\rho }_{\mathrm{h}}+{\nu }_{\mathrm{c}}{\rho }_{\mathrm{c}}}{{\rho }_{\mathrm{h}}+{\rho }_{\mathrm{c}}}=\frac{2\eta }{{\rho }_{\mathrm{h}}+{\rho }_{\mathrm{c}}}\to \eta =\frac{\overline{\nu }({\rho }_{\mathrm{h}}+{\rho }_{\mathrm{c}})}{2}.\end{array}$$(23)

Taking the expression of critical kinematic viscosity (see Eq. (10)), and assuming that $\overline{\nu }$ corresponds to the averaged value of the two systems, we can express the value of the critical dynamic viscosity of our system as

$$\begin{array}{c}\hfill {\eta }_{\mathrm{Crit}}=\frac{\lambda {\rho }_{\mathrm{h}}{\mathcal{M}}_{\mathrm{h}}{c}_{\mathrm{s},\mathrm{cold}}\chi }{200},\end{array}$$(24)

where the wavelength (λ), the density and Mach number of the hot gas (ρ_h,  ℳ_h) and the sound speed of the cold gas (c_s, cold) are kept constant.

Previous studies have tested the critical viscosity estimate for different overdensities, but the results deviated from the predictions at high overdensities since they did not consider the difference in momentum (e.g., Esch 1957; Roediger et al. 2013). For high overdensities each fluid carries a very different momentum (inertia) which viscosity needs to compensate for. By using an harmonic (weighted) mean, we are taking this effect into account.

To test the dependence of a critical viscosity on the overdensity χ, we ran different simulations changing the value of χ and measure the critical viscosity in each case, using the 2D setup described in Sect. 3.1.1. Since we included a unique dynamic viscosity in our setup, we quantified the critical dynamic viscosity (Eq. (24)). Figure 1 shows the numerical results (red dots) normalized to the Spitzer value of viscosity at T = 10⁶ K, together with different methods to average the kinematic viscosity besides the harmonic mean (black line): considering only the low density fluid $\overline{\nu }=\eta /{\rho }_{\mathrm{h}}$ (blue line); considering only the high density fluid $\overline{\nu }=\eta /{\rho }_{\mathrm{c}}$ (red line); and finding an arithmetic mean $\overline{\nu }=({\nu }_{\mathrm{h}}+{\nu }_{\mathrm{c}})/2$ (green line).

Fig. 1. Dependence of ηCrit on the overdensity, χ. Red dots show the numerical results. The blue line shows the theoretical expectation assuming that only the low density is relevant for kinematic viscosity. The red line shows that only the high density is relevant. The green line assumes an arithmetic mean of kinematic viscosity between the fluids. The black line shows a harmonic (weighted) mean. Note that we normalized the values of viscosity to the Spitzer value at 106 K.

The cold medium becomes denser when increasing the density contrast. This has two effects: (i) the cold medium has a much lower kinematic viscosity, ν_c = η/ρ_c. (ii) The cold medium carries much more momentum than the hot medium. For the high overdensity flows χ ≫ 1 we consider here, it is a good approximation to assume that the viscosity of the cold gas is negligible (as we show explicitly in Sect. 5.1, the hot gas viscosity is what matters), but that it contains almost all the momentum. Thus, the critical viscosity increases with overdensity, η_crit ∝ χ (see Eq. (24) and Fig. 1) because the total momentum in the flow is higher and stronger friction is required to have an effect.

From the values of critical viscosity obtained and using the harmonic mean, we can compute their corresponding Reynolds number as

$$\begin{array}{c}\hfill {\mathrm{Re}}_{\mathrm{Crit}}=\frac{\lambda v_{\mathrm{shear}}}{{\overline{\nu }}_{\mathrm{Crit}}}=\frac{\lambda v_{\mathrm{shear}}({\rho }_{\mathrm{h}}+{\rho }_{\mathrm{c}})}{2{\eta }_{\mathrm{Crit}}}=100\frac{1+\chi }{{\chi }^{1/2}},\end{array}$$(25)

namely, Re_Crit ≈ 100χ^1/2 for large χ.

4.2. Critical viscosity versus Mach number

According to our analytical estimate, the critical viscosity depends linearly on the Mach number for a given overdensity (see Eq. (24)). To test this relation, we ran different simulations varying ℳ_h, while keeping χ = 10 fixed and we measured the critical viscosity in each case, using the setup described in Sect. 3.1.2. However, the theoretical analysis was made assuming incompressible fluids, which is a good approximation for low Mach number, but breaks down at high Mach numbers, where compressibility effects become significant.

When compressibility is taken into account, the classical dispersion relation of the KHI derived from linear analysis is no longer valid (Landau 1944; Landau & Lifshitz 1987). Instead, the dispersion relation for compressive fluids in a planar sheet is

$$\begin{array}{c}\hfill [\chi {(\varpi -1)}^2-{\varpi }^2]·[\chi {(\varpi -1)}^2({\mathcal{M}}_{\mathrm{h}}^2{\varpi }^2-1)-{\varpi }^2]=0,\end{array}$$(26)

where ϖ ≡ ω/kΔv_shear (see Mandelker et al. 2016 for details).

The two roots of the quadratic part of the equation are always real; therefore, they do not describe the growth of the KHI. Two of the four roots in the quartic factor are also real, while two are complex conjugates. We took the imaginary part of both complex conjugates, which describe an exponential growth (positive sign) or decay (negative sign) of the instability. The positive imaginary numerical solutions of Eq. (26) for χ = 5, 10, 100 are shown in Fig. 2.

Fig. 2. Numerical solution for the values of Im(ϖ) of Eq. (26) as a function of ℳh for χ = 5, 10, 100. The dashed lines show the ℳh above which the KHI is stable.

In the case of χ = 10, the instability grows (positive Im(ϖ)) for ℳ_h ≲ 1.77. Above this critical value of Mach number (ℳ_Crit), the solution becomes real and the KHI is stable (vertical dashed lines). The general solution for ℳ_Crit (Mandelker et al. 2016) is given by

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{Crit}}={(1+{\chi }^{-1/3})}^{3/2}.\end{array}$$(27)

The value of ℳ_Crit for χ = 10 is computed numerically in Fig. 3 for nonviscous fluids using the method described in Appendix A, finding a ℳ_Crit = 1.65 ± 0.5.

Fig. 3. Growth of the KHI with time for χ = 10, color-coded by ℳh. The solid lines indicate the simulations in which the KHI is still unstable and the dash-dotted lines where the instability is suppressed (values above ℳCrit), finding a ℳCrit = 1.65 ± 0.5.

For ℳ_h > ℳ_Crit, the instability is suppressed and the growth rate decays. An intuitive explanation for this is that when the instability starts growing, areas of high and low-pressure are formed around the perturbation at a rate of approximately the sound-crossing time. These areas feed the instability, leading to its growth. At large mach number, the shear motion of the fluids is faster than the sound-speed, preventing these areas of high and low pressure from forming and leading to the flattening and stabilization of the wave.

For ℳ_h > ℳ_Crit, the instability is fully suppressed by the compressive nature of the fluids and the viscosity needed to keep the system stable is zero. Therefore, there is a transition from the incompressible limit where the suppression is dominated by viscosity to the compressible limit set by ℳ_Crit. An analytical solution for η_Crit as a function of ℳ_h can be estimated in a similar way, as in Sect. 2.1.

In the incompressible limit (ℳ_h → 0), the positive imaginary solution of Eq. (26) yields

$$\begin{array}{c}\hfill \mathrm{Im}(\varpi )=\frac{\sqrt{\chi }}{(\chi +1)},\mathrm{Im}(\omega )=k\mathrm{\Delta }v\frac{\sqrt{\chi }}{(\chi +1)},\end{array}$$(28)

recovering the KH timescale (Eq. (7)), which, in its general form, is defined as

$$\begin{array}{c}\hfill {\tau }_{\mathrm{KH}}=\frac{1}{k\mathrm{\Delta }{v}_{\mathrm{shear}}\mathrm{Im}(\varpi )}.\end{array}$$(29)

By using the numerical solution for Im(ϖ) (shown in Fig. 2) and evaluating τ_KH > τ_ν (i.e., the same analysis as in Sect. 2.1 for the incompressible case), we obtained an analytical solution for a critical viscosity in the compressible case,

$$\begin{array}{c}\hfill {\eta }_{\mathrm{Crit}}=\frac{\lambda {\rho }_{\mathrm{h}}{\mathcal{M}}_{\mathrm{h}}{c}_{\mathrm{s},\mathrm{cold}}\sqrt{\chi }(\chi +1)\mathrm{Im}(\varpi )}{200}.\end{array}$$(30)

In Fig. 4 we show the analytical estimate of η_Crit, normalized to the Spitzer value at T = 10⁶ K, as a function of ℳ_h for incompressible fluids (dashed black line) and for compressible fluids (dashed blue line). The red data points are the numerical values of η_Crit obtained from the simulations. Although the fit is not perfect due to the large assumptions made in the estimate, the initial growth at low ℳ_h, and the decay when ℳ_h ∼ ℳ_Crit match the theoretical estimate. However, there is a drop at ℳ_h ∼ 1, indicating that the KHI becomes stable at transonic speeds. Although the origin of this drop is not clear, one possibility is that higher order modes of the KHI have different values of η_Crit. We will explore this purely adiabatic suppression in a future work.

Fig. 4. Dependence of ηCrit on the Mach number, ℳh. The numerical results (red dots) deviate from the analytical incompressible ηCrit(ℳh) (Eq. (24); dashed black line). However, this is corrected by taking into account the suppression due to high ℳh (Eq. (30), dashed blue line). Note that we normalized the viscosities to the Spitzer value at 106 K.

4.3. Turbulent radiative mixing layers

To study how viscosity affects the TRMLs, we made use of the value of critical viscosity and Reynolds number obtained in Sect. 4.1 for χ = 100: η_Crit = 6425 ± 25 and Re_Crit = 248.6 ± 2.4 respectively. In terms of Spitzer viscosity, the critical viscosity measured corresponds to η_Crit ≃ 3.55η_Sp at a temperature of T = 10⁶ K. We express the amount of viscosity in our setup as a function of the critical viscosity needed to suppress the instability. Using the 3D setup described in Sect. 3.2, we ran a simulation with no viscosity for each Damköhler number Da (labeled “No visc”), with 5% of the critical viscosity (“0.05η_Crit”), 10% of critical viscosity (“0.1η_Crit”) and 50% of critical viscosity (“0.5η_Crit”). Although the definition of critical viscosity considers only the growth of the KHI within 1τ_KH, in this step, we ran the simulations for longer times for the proper development of turbulence. The set of simulations with different viscosities, as well as the value of the Reynolds number in each case are shown in Table 1.

It is important to note that in the “No visc” run, although we did not include any explicit physical viscosity, there is still numerical viscosity given by resolution. Since our resolution consists of 64 × 64 in the $\widehat{x}\widehat{z}$ plane and the length of the box is L_x = L_z = λ, the cell length is Δx = λ/64, which sets the Kolmogorov scale of our “No visc” runs. The runs with explicit viscosity set the Kolmogorov scale at larger values, given by percentage of critical viscosity taken in each case. Since we have not included explicit thermal conduction, the Field length of thermal conduction is given by δ = Δx = λ/64. This allows us to explore the Karlovitz number (see Sect. 1) parameter space from Ka = 1 (“No visc” run) toward smaller values of Ka.

4.3.1. Front morphology

In Fig. 5 we show how viscosity affects mixing in the different runs for a low Da (weak cooling, top panels) and high Da (strong cooling, bottom panels). Each panel shows a temperature slice through the box, centered at the interface between the two fluids at t = 2.5τ_KH. In the weak cooling regime (low Da, top panels), turbulence has developed in the run without explicit viscosity whereas it is more suppressed the more viscous the media are (left to right).

Fig. 5. Temperature slices for the different TRML simulations. From left to right: Non-viscous case, 5% of critical viscosity, 10% of critical viscosity and 50% of critical viscosity. Top: Very weak cooling with Damix ≃ 2.3 × 10−3. Bottom: Very strong cooling with Damix ≃ 1.1 × 103.

It is important to note that the computation of critical viscosity is done during t < τ_KH, but naturally viscosity is affecting the fluid movement for the entire simulation. This implies that, although during t < τ_KH the instability can grow (η < η_Crit), it might eventually lead to a quasi-laminar flow. This can be seen in the top-right panel of Fig. 5, where a laminar regime has already been reached.

In the strong cooling regime (high Da, bottom panels), the results look similar regardless of the amount of viscosity. In all cases the hot and cold medium are clearly separated and the front consists of (approximately) one cell only. Surprisingly, in contrast to what we see in the weak cooling regime, the simulations show some movement in the interface independently of the amount of viscosity. This is analyzed in Sect. 4.3.5.

4.3.2. Effect of viscosity on the cooling efficiency

The luminosity of TRMLs is predominantly due to enthalpy flux (i.e., mass conversion from hot to cold medium) and not due to dissipation of kinetic energy (see Ji et al. 2019; Tan et al. 2021). We checked that this is also true in the case with viscosity (see Appendix B), specifically also that viscous heating plays a negligible role.³

To study the difference in results with and without viscosity in a quantitative way, we need to measure the surface brightness Q in each case for each Da_mix. This surface brightness is calculated by measuring the total luminosity of the box L_rad and dividing it by the area of the section of the box (Q = L_rad/λ²).

In Fig. 6 we show the effect of viscosity for the simulations with weakest cooling (left panel) and the simulations with the strongest cooling (right panel) considered. Note that the dark blue lines correspond to the “No visc” run, considering the estimated numerical viscosity. Although visually viscosity suppresses the development of turbulence in the weakest regime (see the top panel of Fig. 5), the surface brightness obtained after several τ_KH is very similar in all cases.

Fig. 6. Evolution of the surface brightness due to cooling in the mixing layers. Left: Evolution of Q color-coded by viscosity in a very weak cooling regime with Damix ≃ 2.3 × 10−3. Right: Evolution of Q for different viscosities in a very strong cooling regime with Damix ≃ 1.1 × 103.

In the case of strong cooling (right panel of Fig. 6)⁴, the differences become greater, with an increase of Q with viscosity. This might be counterintuitive if we are thinking of viscosity as a force suppressing the turbulent movement and, therefore, the cooling process. This is analyzed in more detail in Sect. 4.3.5.

In Fig. 7 we show the relation between surface brightness and Da (Eq. (12)). Our simulations without viscosity reproduce the expected broken power-law behavior (see Eq. (12)) with the change of regime from weak to strong cooling occurring at Da_mix ≃ 1⁵. Although in the weak cooling regime, the differences are small, in the strong cooling regime viscosity leads to larger values of Q. This implies that the transition from weak to strong cooling regime is shifted toward larger values of Da_mix when increasing the amount of viscosity.

Fig. 7. Surface brightness, Q, as a function of Damix color-coded by viscosity. The darker dashed line shows the best fit in each case in the weak cooling regime (Da < 1) with a slope of α = 1/2. The lighter dashed line shows the best fit in each case in the strong cooling regime (Da > 1) with a slope of α = 1/4.

In summary, both Figs. 6 and 7 reveal two counterintuitive behaviors: while we expect viscosity to suppress the mixing process and, thus, the amount of intermediate temperature gas available for cooling, we find (i) for the weak cooling regime a mass transfer rate between the phases independent of viscosity, and (ii) for the strong cooling regime even an increased ṁ with larger values of viscosity. This implies that either the naive expectation of how viscosity affects turbulence is wrong, or the u′–Q relation is changed due to viscosity. We went on to explore this by analyzing the turbulence inside the mixing layer in the next sections.

4.3.3. Effect of viscosity on the turbulence

From visual inspection of Fig. 5, viscosity seems to suppress turbulence in the weak cooling regime, but not in the strong cooling regime. To test this in a quantitative way, we computed the turbulent velocity in our domain following Tan et al. (2021); namely, to avoid the effect of bulk motions, we measured the velocity dispersion in the $\widehat{z}$ direction, by taking the maximum of v_z profile for each snapshot and then averaging over several τ_KH⁶.

Figure 8 shows relation 12, between Q and the turbulent velocity normalized to the cold gas sound speed, color-coded by Da_mix. To get rid of the Da_mix dependence, we normalize the value of Q to the value of ${\stackrel{\sim }{\mathit{Da}}}^{\alpha }$, as done in Das & Gronke (2023). $\stackrel{\sim }{\mathit{Da}}$ is calculated from Eq. (11), using the turbulent velocity measured from the simulations and α corresponds to the index from Eq. (12). By doing this normalization, $Q/{\stackrel{\sim }{\mathit{Da}}}^{\alpha }\propto u\prime$. This relation is satisfied in the strong cooling regime (Da_mix > 1), regardless of the amount of viscosity. However, in the weak cooling regime (Da_mix < 1), the relation is satisfied for the “No visc” (circles) and 0.05η_Crit (triangles) cases. The turbulent velocity of the runs with 0.1η_Crit (squares) are one order of magnitude lower than expected and the runs with 0.5η_Crit (crosses), more than two orders of magnitude lower.

Fig. 8. Relation between Q, normalized by the numerical measured Damköhler number ($\stackrel{\sim }{\mathit{Da}}$) in each regime, and the turbulent velocity for different viscosities and color-coded by Damix. The dashed line shows the relation $Q/{\stackrel{\sim }{\mathit{Da}}}^{\alpha }\propto u\prime$ expected.

This result indicates that, as previously noted by visual inspection in Fig. 5, viscosity suppresses turbulence, but it does not reduce surface brightness. However, the suppression occurs in the weak cooling regime, whereas in the strong cooling regime, the relation between Q and u′ is not modified. This is surprising, due to the weak dependence of u′ on Da (u′∝Da^1/10; see Tan et al. 2021).

We also tested this weak dependence of the cooling on turbulence explicitly and give the results in Fig. 9. The simulations without viscosity fall on the u′∝Da^1/10 relation found in Tan et al. (2021). In the weak cooling regime (Da_mix < 1) the relation is not satisfied in the viscous runs. The turbulent velocities drop until they reach the u′ in adiabatic simulations, saturating at this value (dashed blue lines). The adiabatic system sets the floor for turbulent velocities depending on viscosity.

Fig. 9. Dependence of the turbulent velocity on Damix for different amounts of viscosity and color-coded by the surface brightness. The dashed black line indicates the dependency of $u\prime \propto {{\mathrm{Da}}_{\mathrm{mix}}}^{1/10}$. The dashed blue lines show the turbulent velocity values of the adiabatic system with different viscosities.

In this weak cooling regime, the simulations show the expected behavior: the more viscous the system is, the less turbulent it becomes. However, for larger values of Da_mix > 1, all the runs tend to follow the u′∝Da^1/10; namely, the runs mostly display a level of turbulence (a) that is much higher than expected and (b) fairly independent of viscosity.

Furthermore, Fig. 9 also shows that the characteristic Damköhler number at which the simulations follow the expected u′∝Da^1/10 relation increases with increasing viscosity (i.e., the runs with lower viscosity start following the relation at a lower Da_mix), while the more viscous runs occur at a higher Da_mix. In the strong cooling regime, all the simulations follow the u′∝Da^1/10 relation, suggesting that the effect of viscosity is reduced due to cooling.

Viscosity acts as a friction, reducing the velocity of the gas and, therefore, it is expected to produce lower values of turbulent velocity. However, two questions arise, which we explore in the following subsections.

1. Viscosity suppresses turbulence in the weak cooling regime (Da < 1) – as expected. However, we need to investigate why the amount of cooling is only weakly affected (see Fig. 7), thereby effectively breaking the universal u′−Q relation (see Fig. 8).

2. On the other hand, viscosity has a weak effect on both the level of turbulence as well as the amount of cooling in the strong cooling regime (see Figs. 8 and 9). This is particularly puzzling, since turbulence only depends very weakly on cooling (see Fig. 9). Thus, even with very strong cooling turbulence should not directly “overpower” the effect of viscosity.

4.3.4. Weak cooling regime

To address the first question, namely, why for Da < 1 the cooling rate is not affected by viscosity (but u′ is), we need to understand the behavior of a viscous gas in the weak cooling regime. Due to the long cooling time, the gas behaves similarly to an adiabatic gas, which has been previously studied in detail (e.g., Roediger et al. 2013; Marin-Gilabert et al. 2022). Viscosity smooths out the shear velocity gradient and suppresses the growth of instabilities, leading to a reduced turbulent velocity (see the top panel of Fig. D.1). However, if turbulence is not the driver for mixing, we must explore why Q is on the same order of magnitude regardless of the amount of viscosity.

When viscosity is strong enough, the instabilities are suppressed, leading to a laminar flow, rather than a turbulent flow. The theoretical relation in Fig. 8 assumes turbulent mixing (e.g., Damköhler 1940; Shchelkin 1943), therefore it is not surprising that for high viscosities the simulations do not follow the expected relation. Viscous diffusion dominates over convection and cooling, increasing the thickness of the interface, producing that the theoretical analysis fails (e.g., Klimov 1963; Libby & Williams 1982). Although the viscous runs present a laminar flow, the amount of gas at intermediate temperatures where the cooling curve peaks is similar among the viscous cases (see the top row of Fig. C.1). The nonviscous case leads to a higher amount of intermediate temperature gas due to turbulent mixing. However, the long cooling time in the weak cooling regime produces a bottleneck effect, which explains the constant Q obtained, despite the different diffusion mechanisms. Note that our criterion for critical viscosity is defined within 1τ_KH, but our simulations in the weak cooling regime are run until 50τ_KH. This means that, although within 1τ_KH the instability might grow, the system might become laminar after several τ_KH.

4.3.5. Strong cooling regime

To answer the second question, namely, why neither u′, nor the cooling are strongly affected by viscosity for Da > 1, we need to study the strong cooling regime in detail. Viscosity is expected to reduce turbulence, although in the strong cooling regime, the Q − u′ relation follows the analytical expectations from turbulent theory. Thus, if viscosity is indeed expected to suppress turbulence, we must consider why u′ is similar in all cases regardless of the amount of viscosity.

In the strong cooling regime, cooling acts much faster than turbulent mixing and ends up quickly cooling all the intermediate gas and producing a front where cooling occurs that is ∼1 cell thick. As a consequence, cooling tends to keep the temperature profile sharp and, if the system is in pressure equilibrium, also the momentum and shear velocity profiles. This is the opposite effect to viscosity, which tends to smooth out the shear velocity gradient. This produces that, in the strong cooling regime, cooling dominates over viscosity, keeping the shear velocity profile sharp and, as a consequence, being able to induce turbulence in a similar way as the nonviscous case and reducing the effect of viscosity (see Fig. D.1).

Keeping the shear profile sharp leads to continuous KHIs forming, and hence to turbulence, as we can see in the top panel of Fig. 10, where we show the dependence of turbulence on the scale d at which the shear velocity profile is smoothed after t = 1.25τ_KH (see Appendix E for details). For “No visc” and 0.05η_Crit, d is always small, producing larger values of u′ (i.e., turbulence). For the cases of 0.1η_Crit and 0.5η_Crit, the weak cooling regime results show the effect of viscosity in smoothing out the gradient (large values of d), leading to a low u′. However, in the strong cooling regime the shear velocity gradient becomes sharp (low values of d), inducing turbulence. In summary, while cooling does not (strongly) affect the turbulence directly, it sharpens the density gradient required to seed KHI, and hence combats the smoothing effect of viscosity.

Fig. 10. Top: Turbulent velocity as a function of the velocity gradient width, color-coded by the surface brightness. Bottom: Surface brightness as a function of the velocity gradient width, color-coded by Damix.

In other words, a sharp velocity gradient induces turbulence and, therefore, leads to larger surface brightness. This can be seen in the bottom panel of Fig. 10, where larger values of Da_mix lead to sharper gradients (lower d) and, therefore, larger Q.

To test our assumption that the system follows the expected u′∝Da^1/10 when cooling dominates over viscosity, we estimate analytically at which value of Da_mix this should happen. Cooling dominates over viscosity when t_cool < τ_ν, where t_cool is measured at T_mix = 2 × 10⁵ K and τ_ν corresponds to the time that a given viscosity takes to smooth out the velocity gradient a distance d = λ/10 (see Sect. 2.1).

In our case, we express the different amounts of viscosity as a function of the critical viscosity as η = ζ η_Crit, where ζ is the fraction of the critical viscosity used. Using the equation for critical viscosity (Eq. (10)), KH time (Eq. (7)), and ν = η/ρ, we can express η as

$$\begin{array}{c}\hfill \eta =\zeta {\eta }_{\mathrm{Crit}}=\zeta \frac{{\lambda }^2\rho }{100{\tau }_{\mathrm{KH}}}.\end{array}$$(31)

We can also express t_cool as

$$\begin{array}{c}\hfill t_{\mathrm{cool}}=\frac{{\tau }_{\mathrm{KH}}}{\mathrm{Da}},\end{array}$$(32)

where t_turb = λ/u′∼τ_KH and the viscous timescale as

$$\begin{array}{c}\hfill {\tau }_{\nu }=\frac{d^2}{\nu }=\frac{d^2\rho }{\eta }=\frac{{\lambda }^2\rho }{100\eta }=\frac{{\tau }_{\mathrm{KH}}}{\zeta }.\end{array}$$(33)

Therefore, cooling will dominate over viscosity when

$$\begin{array}{c}\hfill \mathrm{Da}>\zeta .\end{array}$$(34)

We tested this estimate with our simulations by checking the critical value of Da_mix (which we denote as Da^*), when the simulations start to follow the u′∝Da^1/10 relation for the different amounts of viscosity (see Fig. 9 in which can be read off of Da^* directly). We included two more sets of simulations for completeness and we plot the computed values of Da^* as a function of ζ in Fig. 11.

Fig. 11. Measured critical Da (Da∗) vs. the fraction of critical viscosity (ζ) for different sets of simulations. We included two more sets for completeness that are not shown in the previous plots. The dashed line indicates the expected one-to-one relation.

Although the theoretical estimate was based on broad assumptions, the computed values of Da^* from our simulations fit the theoretical estimate. This agrees with our previous assumption that cooling effectively reduces the effect of viscosity by keeping the temperature gradient sharp.

To understand why cooling is increased (by a factor of ∼2) for the more viscous runs at Da > Da^∗, we start by recalling that while cooling reduces the effect of viscosity, the viscous runs still lead to a broader front (see the bottom row of Fig. 5). As a consequence, the amount of gas with an intermediate temperature is larger with viscosity (see the bottom row of Fig. C.1). This produces that the overall amount of cooling is larger with viscosity than without viscosity, leading to a slightly higher Q in the strong cooling regime in the cases with higher viscosities⁷.

To study how u′ depends on viscosity in more detail, Fig. 12 shows how turbulent velocities change with viscosity. In the nonviscous case, the data points cluster together, corresponding to the weak dependence of u′ on Da. For higher viscosities, the two regimes defined by Da^∗(η) can be identified in each case: the data points group separately in the viscous-dominated regime (Da_mix < ζ) and in the cooling-dominated regime (Da_mix > ζ).

Fig. 12. Turbulent velocity as a function of the amount of viscosity of the medium (Karlovitz number), color-coded by Damix.

Assuming that u′ depends on viscosity by a power law u′∝η^α, we can fit the slopes for individual Da_mix in Fig. 12 and check the two different regimes. Figure 13 shows the different values of α depending on Da_mix, reflecting the change of regime from viscous-dominated to cooling-dominated at Da_mix ∼ 0.5. In the viscous-dominated regime (low Da_mix) the turbulent velocity decays with viscosity as u′∝η⁻³, while in the cooling-dominated regime (high Da_mix) the exponent is close to zero, indicating that the system behaves similarly regardless of the amount of viscosity⁸.

Fig. 13. Exponent of the dependence of the turbulent velocity on the amount of viscosity (Karlovitz number) calculated from Fig. 12 as a function of Damix.

5. Discussion

5.1. Temperature-dependent viscosity

Considering a more realistic Spitzer viscosity dependent on the temperature of the gas, each fluid will have a viscosity η. In the case of the cold fluid: ${\eta }_{\mathrm{Cold}}\propto T_{\mathrm{Cold}}^{5/2}$; while in the hot medium: ${\eta }_{\mathrm{Hot}}\propto T_{\mathrm{Hot}}^{5/2}$. Therefore, by comparing both viscosities, we have

$$\begin{array}{c}\hfill {\eta }_{\mathrm{Hot}}={\eta }_{\mathrm{Cold}}{\left(\frac{T_{\mathrm{Hot}}}{T_{\mathrm{Cold}}}\right)}^{5/2}={\eta }_{\mathrm{Cold}}{\chi }^{5/2}.\end{array}$$(35)

In the case of an overdensity of χ = 100: η_Hot = 10⁵ η_Cold, rendering the viscosity of the cold medium negligible.

Additionally, the diffusivity of a system depends on the diffusion coefficient, which in the viscous case is set by the kinematic viscosity (ν). For a system with a constant dynamic viscosity (η, as in our setup) but two different densities, each fluid will have a different diffusivity (see Eq. (8)). Since the kinematic viscosity depends on the inverse of the density, for a given dynamic viscosity, the cold medium will have a lower diffusivity compared to the hot medium. As a consequence, for systems with a high density contrast, the diffusivity of the system is set by the hot gas.

Taking into account that: with a more realistic temperature-dependent Spitzer viscosity, the dynamic viscosity of the cold medium is negligible and that the viscous diffusion is set by the hot gas, a viscous cold medium should have a small impact compared to a viscous hot medium. We tested this using 1D simulations, with the following setup:

$$\begin{array}{c}\hfill (\rho ,T,v,\eta )=\{\begin{array}{cc}({\rho }_{\mathrm{Cold}},T_{\mathrm{Cold}},-v_x,{\eta }_{\mathrm{Cold}})\hfill & \text{for}y<0,\hfill \\ ({\rho }_{\mathrm{Hot}},T_{\mathrm{Hot}},v_x,{\eta }_{\mathrm{Hot}})\hfill & \text{for}y>0,\hfill \end{array}\end{array}$$(36)

where Δv_x/c_s = 0.24, we keep T_Hot and ρ_Hot constant and change density (and temperature) of the cold medium as T_Hot/T_Cold = ρ_Cold/ρ_Hot = 2, 5, 10, 100, 1000. To analyze the diffusivity depending on the medium, we set a constant viscosity of η = 25 (in internal units) first in the cold medium and then in the hot medium. Figure 14 shows the v_x profile after Δv_xt/L = 2.50 for the two cases: the viscous cold medium on the top panel and the viscous hot medium on the bottom panel.

Fig. 14. Velocity profile of 1D simulations. Top: Only the cold gas (y < 0) is viscous. Bottom: Only the hot gas (y > 0) is viscous.

When the cold medium is viscous (top panel in Fig. 14), diffusivity is decreased as the cold medium becomes denser. Thus, if the initial (i.e., lowest overdensity) value of viscosity is too low to have a discernible effect, it becomes increasingly negligible at higher overdensity. When the hot medium is viscous (bottom panel in Fig. 14), diffusivity remains constant, since the density of the hot medium is kept constant. However, the velocity profile becomes smoother as the cold medium becomes denser, despite the cold medium not being viscous. This happens because, although the cold medium is not viscous, it carries a lot more momentum than the hot (and viscous) medium. Therefore, when the cold medium transfers momentum to the hot medium, viscosity changes the momentum of the hot medium very effectively, smoothing the velocity gradient. The denser the cold medium is, the more momentum is transferred to the hot medium, leading to a smoother velocity profile. These results corroborate the assertion reported in Sect. 4.1, namely, that the hot phase dominates viscosity and the cold phase dominates momentum. We note that this leads to the interesting effect that higher overdensities imply a larger degree of smoothing of the velocity profile, however, as η_crit also rises (see Eq. (24) and Fig. 1), a higher viscosity is required to suppress the KHI.

Considering a constant η for the whole system leads to a diminishing momentum diffusion in the cold medium, similarly to a realistic temperature-dependent Spitzer viscosity. Therefore, a constant η would lead to more realistic results than considering a constant ν, where both media would be equally diffusive.

5.2. Super-critical viscosity

In this paper, we focus on the sub-critical viscosity regime, where all the cases studied have a viscosity smaller than the critical viscosity. The instability is able to grow within one Kelvin-Helmholtz time, being enhanced by strong cooling and allowing to keep the turbulence (see Sect. 4.3.5).

With a super-critical viscosity (η > η_Crit) the instability is fully suppressed and there is no initial growth. This means that cooling cannot enhance the growth of the instability and dominate over viscosity. To test this, we ran a simulation with η = 2η_Crit at Da_mix ≃ 1.1 × 10 ², which is an order of magnitude larger than the estimated Damköhler number required to overpower viscosity Da^∗. Since the viscous length is larger than the wavelength for super-critical viscosity, we increase the size of the box to L_x = L_z = 3λ.

In Fig. 15 we show the shear velocity profile after t = 1.25τ_KH for different levels of viscosity. For η < η_Crit the effect of cooling keeps the velocity profile sharp, while for η = 2η_Crit viscosity is able to smooth out the gradient, although Da_mix ≫ Da^∗.

Fig. 15. vshear profiles at Damix ≃ 1.1 × 102 and t = 1.25τKH including the run with 2ηCrit. The distance, d, at which the shear profile has been smoothed in the 2ηCrit case is large compared to the other cases.

As a consequence, the turbulent velocity is suppressed and viscosity still dominates over cooling. The turbulent velocity of the run with 2η_Crit is $u\prime /{\mathrm{c}}_{\mathrm{s},\mathrm{cold}}=9.4\times {10}^{-7}$, behaving as viscous-dominated (see Fig. 9 for reference).

When η < η_Crit, there is an initial growth of the KHI that is later enhanced by cooling, leading to a cooling-dominated regime where the system behaves as nonviscous. However, if η > η_Crit there is not initial growth that can be enhanced by cooling, therefore the system will always be viscous-dominated.

In the temperature domain we consider in this paper, where the cooling function peaks (10⁴ K < T < 10⁶ K), the Spitzer value for the 2η_Crit is 2η_Crit ≃ 7.10η_Sp (see Sect. 4.3 and Table 1), which is not realistic in plasma physics (Kunz et al. 2022). Therefore, an astrophysical scenario with such a high level of viscosity is very unlikely.

5.3. Applications

Viscosity becomes relevant in very hot and ionized media, with temperatures typically above 10⁶ K. On the other hand, cooling efficiency peaks between 10⁴ − 10⁶ K. This means that both cooling and viscosity will act together in astrophysical scenarios where very hot and cold media co-exist and mix, which is the case for the ICM, ISM, and CGM – but also in other multiphase systems such as galactic winds or coronal rain. This interaction might happen in different astrophysical processes, where the Reynolds number is comparable to the ones used in this paper (see Table 1). Similarly, mixing and radiative cooling are physical processes that occur frequently in these multiphase systems. In fact, viscosity, mixing and cooling are commonly thought of as being crucial in determining their long-term evolution. Examples of such multiphase systems include:

• Galaxies moving through the ICM experience a ram pressure that strips the cold ISM (e.g., Gisler 1976; Nulsen 1982). The growth of instabilities in the tails mixes the hot ICM with the cold ISM (e.g., Iapichino et al. 2008; Ghosh et al. 2024). As a consequence, the hot ICM cools down, providing new gas and enhancing star formation (Müller et al. 2021). Viscosity not only affects the morphology and mixing of the tail with the medium (e.g., Kraft et al. 2017; Marin-Gilabert et al. 2024), but also the cooling efficiency and, therefore, it might change the star formation rate.

• The existence of multiphase gas within galactic winds is commonly often attributed to the effects of continuous cooling (see the reviews by Veilleux et al. 2020 and Thompson & Heckman 2024). Similarly, cold clouds infalling through the CGM is another example where the hot gas of the CGM mixes with the cold gas of the cloud being disrupted (Gronke & Oh 2018; Gronke et al. 2021). When the CGM is hot enough (T > 10⁶ K), transport processes become relevant and shape the morphology of the tails of the cloud (Brüggen et al. 2023). Viscosity affects the cloud dynamics, as well as the evolution and evaporation of the tail (Jennings & Li 2021). Therefore, it also affects the amount of cold gas and the mixing efficiency of the cloud and the CGM. The mixing between the cold stripped gas and the medium can be observed and studied using different ion lines (e.g., Kwak & Shelton 2010).

• Cold streams in hot galactic halos are likely crucial to our understanding of galaxy fueling and growth (Dekel et al. 2009). As other classical example of multiphase gas interactions, mixing, cooling, and the effect of viscosity are discussed here in the context of cold streams (Brüggen & Ruszkowski 2005; Mandelker et al. 2020).

In each of these systems, turbulence, mixing, and cooling are crucial to understanding the dynamics and ultimate fate of the gas. As a result, key findings from studies of TRMLs (e.g., Ji et al. 2019; Fielding et al. 2020; Tan et al. 2021) are commonly invoked to interpret a wide range of observational diagnostics (e.g., Qu et al. 2022; Lin et al. 2025). However, most previous works have neglected the role of viscosity, meaning that such interpretations could be significantly biased; in particular, if viscosity were to strongly suppress turbulent mixing and, consequently, cooling, as might be naively expected.

In this study, however, we show that (contrary to expectations) viscosity has only a modest impact on the dynamics of turbulent mixing layers. In particular, for systems in the fast cooling regime (Da > 1), which includes most of the astrophysical environments described above, not only does the u′–Q relation remain unchanged, but also the amount of turbulence expected from KHI. This allows us to continue using established theoretical predictions to relate the line widths to expected emissivities and other properties.

Another important application of our work relates to the growing number of multiphase subgrid models, which aim to address the long-standing resolution problem in large-scale (e.g., cosmological) simulations. Such simulations are unable to fully resolve multiphase structures, such as those discussed above, and therefore rely on sub-grid prescriptions to approximate their behavior (Huang et al. 2020; Smith et al. 2023; Butsky et al. 2024; Das et al. 2025). These models commonly adopt a numerical multi-fluid framework (e.g., Weinberger & Hernquist 2022), with source and sink terms for each fluid component informed by small-scale simulations of turbulent mixing layers (e.g., Tan et al. 2021).

However, to date, none of these models have accounted for the combined effects of radiative cooling and viscosity. Instead, they implicitly assume that the omitted effect (e.g., viscosity) has a negligible influence on the relevant dynamics. In this work, we provide a robust physical foundation for current and future implementations of multiphase sub-grid models by explicitly testing this assumption. In particular, we show how mass transfer rates between phases are (largely) unaffected by viscosity, thus justifying the use of existing cooling-based prescriptions in multi-fluid models.

Taken together, our results show that viscosity, while potentially important in shaping large-scale morphology, does not significantly alter the fundamental scaling relations governing turbulent mixing and radiative cooling. Although we might naively expect viscosity to suppress mixing and thereby invalidate commonly used theoretical predictions, we find (perhaps surprisingly) that this is not the case. This lends confidence to both the interpretation of observational diagnostics and the construction of multiphase subgrid models, which often rely on idealized, inviscid small-scale simulations. By clarifying the limited role of viscosity, our work helps bridge the gap between detailed local physics and the broader astrophysical systems in which these processes operate.

5.4. Caveats and future work

This paper focuses on a very idealized setup for simplicity, to isolate the effect of viscosity and cooling. However, future projects can expand this work further, for instance, by focusing on:

• Temperature-dependent viscosity: In this study, we looked at the physical effects at play and, thus, we used a constant viscosity in our simulations for simplicity. However, Spitzer viscosity is expected to depend on the temperature of the medium, leading to a more complex system. Although viscosity is expected to be low in the cold medium, the momentum transfer from the cold and denser medium to the viscous hot medium can lead to non-negligible effects (see Sect. 5.1). Additionally, the mass inflow due to cooling might also be affected, since gas with different temperatures will have different viscosities. Temperature-dependent viscosity could also potentially affect the steady-state amount of gas at different temperatures (i.e., the temperature probability density function), similarly to how the temperature dependence of thermal conduction strongly affects the temperature probability density function (PDF) (Tan & Oh 2021). The temperature PDF is crucial to predicting absorption and emission line ratios. Interestingly, the absolute value of viscosity appears to affect the temperature PDF (Fig. C.1), whereas the absolute value of conduction does not appear to affect the temperature PDF, which only depends on how conduction scales with temperature (Tan & Oh 2021). Further investigations of such effects are beyond the scope of this paper.

• Anisotropic viscosity: In magnetized media such as the CGM, where the particles’ gyro-radius is small compared to the mean free path, the effect of viscosity is expected to be anisotropic. This leads to a lower viscosity perpendicular to the magnetic field lines, decreasing the total effect of viscosity. A more realistic setup should include magnetic fields, which can suppress the growth of instabilities (Das & Gronke 2023), and an anisotropic viscosity, which will be dependent on the magnetic field direction (ZuHone et al. 2015; Berlok et al. 2019; Marin-Gilabert et al. 2025). The combination of both might modify the morphology of the interface, depending on the direction of the magnetic field, potentially changing the cooling efficiency of the system.

• Larger-scale simulations: Here, we focused on an individual mixing layer. However, more complex multiphase setups, such as the “cloud crushing” setup, exhibit turbulent mixing in several locations. Viscosity affects the growth of instabilities, which is essential for determining the survival of a cold cloud (e.g., Klein et al. 1994) subject to a uniform wind. However, it also affects cooling, which is essential in the survival and the mass growth of the cloud in a radiative turbulent medium (e.g., Gronke & Oh 2018; Kanjilal et al. 2020). The change in morphology and cooling efficiency might lead to different criteria to determine the survival of a cloud or the mass growth. We note that Li et al. (2019) conducted such “cloud crushing” experiments including viscosity. They found a prolonged lifetime of the cloud; however, they did not study the (potential change of the) mass transfer rate between the phases in detail.

• Effect of thermal conduction: While thermal conduction is generally subdominant over turbulent diffusion and, thus, it is not thought to affect the overall level of cooling found in TRMLs (Tan et al. 2021), the inclusion of conduction can shift the Karlovitz number (see Sect. 1) to Ka > 1. This might potentially affect the microphysics of the cooling front. Hence, we want to include both viscosity and conduction in the future to put the mass transfer scalings on solid ground.

6. Conclusions

The role of viscosity in TRMLs is key for understanding the suppression of hydrodynamical instabilities and its effect on multiphase gas evolution in astrophysical environments. In this work, we investigated the impact of viscosity on the KHI and its implications for mixing and cooling processes. Using numerical simulations, we systematically varied viscosity, overdensity, and Mach numbers to study the transition from turbulent to laminar regimes and characterized the interplay between viscosity and radiative cooling. Our key findings can be summarized as follows:

• We confirm previous results (e.g., Roediger et al. 2013; Marin-Gilabert et al. 2022) stating that in adiabatic setups, the KHI is suppressed when the dynamical viscosity exceeds a critical threshold. We derived an analytical solution based on the comparison between the viscous timescale and the instability growth time. We numerically confirmed that the critical viscosity depends linearly on overdensity.

• At high Mach numbers in compressive-adiabatic setups, the KHI transitions from a viscosity-dominated suppression to a compressibility-induced stabilization. We verified numerically that the instability is fully suppressed above a critical Mach number, in agreement with a previous analytical work (Mandelker et al. 2016).

• In a multiphase system, viscosity alters the turbulent cascade, modifying the mixing layer properties. In the weak cooling regime, viscosity suppresses turbulence, transitioning the mixing layer into a laminar state dominated by diffusion rather than convection (leading to an altered u′−Q relation). Nevertheless, because cooling is the bottleneck in such a case, the overall luminosity remains unchanged.

• We find (perhaps surprisingly) that in the strong cooling regime, all the TRML simulations have a similar (same order of magnitude) surface brightness and net cooling rates as inviscid simulations, regardless of the amount of viscosity. This counterintuitive behavior arises from cooling sharpening the temperature profile, thereby directly compensating for the smoothing effect of viscosity, allowing for KHI to develop, and keeping the scaling relations in TRMLs intact. We derived an analytical estimate for when viscosity does not affect cooling rates based on the relative timescales of turbulence, viscosity, and cooling. Our numerical results match this theoretical prediction, confirming that cooling can compensate for the viscous suppression in sufficiently strong cooling conditions.

• The suppression of turbulence only mildly affects the energy dissipation in TRMLs. In the weak cooling regime, viscosity reduces turbulent velocities, but it has a weak effect on the surface brightness. In the strong cooling regime, the surface brightness increases up to a factor of ∼2 with viscosity, suggesting that the smoother temperature gradient leads to a larger amount of intermediate temperature gas, thereby enhancing cooling efficiency.

These results quantify the impact of viscosity in astrophysical TRMLs, with important implications for the dynamics of ram-pressure stripping galaxies and galaxy halos. In addition, it also provides key insights into the dynamics of the CGM and ICM.

Acknowledgments

TM would like to thank Hitesh Das, Eugene Churazov, Elke Roediger, Rajsekhar Mohapatra and Klaus Dolag for the intense discussions that motivated this work. The authors also want to thank the referee for their very useful comments. TM acknowledges support by the COMPLEX project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679. MG thanks the Max Planck Society for support through the Max Planck Research Group, and the European Union for support through ERC-2024-STG 101165038 (ReMMU). Most computations were performed on the HPC system Freya at the Max Planck Computing and Data Facility (MPCDF). This research was supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP).

References

Appendix A: Growth of the KHI

The growth or decay of the KHI is calculated using a discrete convolution of the sinusoidal perturbation (McNally et al. 2012) to measure the rate of growth of the KHI. A positive slope within 1τ_KH is considered a growth, while a negative slope a decay. The numerical value of critical viscosity will be the value found between the run with the shallowest positive slope and the one with the shallowest negative slope. Figure A.1 illustrates the measurement of the critical viscosity for an overdensity of 10, where we obtain a critical viscosity of η_Crit = 515 ± 5 in code units. The initial decay is due to the loss of kinetic energy of the fluid moving perpendicular to the shear flow through the flow moving in the opposite direction (Junk et al. 2010).

Fig. A.1. Growth of the instability as a function of time for an overdensity of χ = 10, color-coded for different values of viscosity. The value of the critical viscosity is equal to ηCrit = 515 ± 5, where the error is given by the average between the viscous run with a positive and a negative slope. The solid lines indicate the simulations in which the KHI is still unstable and the dash-dotted line where the instability is suppressed (η > ηCrit).

Appendix B: Mass conversion and luminosity

In radiative mixing layers, the net radiative losses due to cooling are balanced by the enthalpy flux, producing a mass flow from the hot to the cold medium. This means that the total luminosity is proportional to the mass transfer from the hot to the cold medium: Ṁ ∝ L_rad (Ji et al. 2019; Gronke et al. 2021). We checked that the total luminosity measured in our system was due to enthalpy flux and not to dissipation of kinetic energy by measuring the mass transfer from hot to cold gas and the total luminosity (see Fig. B.1). The linear dependence (dashed black line) shows that the luminosity in our system is due to enthalpy flux (mass transfer from hot to cold gas), and not due to dissipation of kinetic energy or viscous heating, regardless of the amount of viscosity and the Da.

Fig. B.1. Mass conversion from hot to cold gas as a function of the total luminosity, color-coded by Da. The mass transfer is proportional to the total luminosity, indicating that the luminosity is due to enthalpy flux.

Appendix C: Intermediate temperature histogram

The cooling curve reaches the highest values at intermediate temperatures, making cooling more effective in this range. This means that the amount of Q for a given Da depends on the amount of gas with an intermediate temperature. Quantifying the amount of intermediate gas gives us a hint on how Q responds in different simulations. Figure C.1 shows the histogram of intermediate temperature gas for different Da_mix (top to bottom: Da_mix ≃ 2.3 × 10⁻³, Da_mix ≃ 2.3 and Da_mix ≃ 1.1 × 10³) and different times (left to right: t = 2.5τ_KH, t = 5τ_KH and t = 6.25τ_KH).

Fig. C.1. Intermediate temperature histograms color-coded by viscosity. From top to bottom: Weak cooling (Damix ≃ 2.3 × 10−3), intermediate cooling (Damix ≃ 2.3), and strong cooling (Damix ≃ 1.1 × 103). From left to right: Histogram measured at t = 2.5τKH, t = 5τKH and t = 6.25τKH.

In the weak cooling regime (top row), the "No visc" run mixes the gas more effectively, leading to a larger amount of intermediate temperature gas. However, since cooling is the bottleneck in this case, this does not make a big impact on Q (see Fig. 8). In the viscous runs, the amount of intermediate temperature gas is very similar. This is due to the smooth momentum profile as a consequence of viscosity. Due to pressure equilibrium, this smoother momentum is translated into a smoother temperature profile, leading to more intermediate temperature gas. Note the flat histogram in the case of 0.5η_Crit due to its laminar flow, where the temperatures are distributed approximately equally along the y-axis. For Da_mix ≃ 2.3 (middle row), the histograms look very similar, since at this value of Da_mix all the simulations are already cooling-dominated, although at later times the "No visc" run seems to produce more intermediate temperature gas. In the strong cooling regime (bottom row), viscosity leads to a slightly broader temperature profile, and therefore more intermediate temperature gas. Although the effect is only a factor of ∼2, the strong cooling produces that this effect can be visible when Q is calculated (see the right panel of Fig. 6).

Appendix D: Shear velocity profiles

In a planar slab (or sheet) setup, the effect of viscosity in the velocity gradient is obtained by solving the Rayleigh problem (also known as Stokes’ first problem; Stokes 1851; Lord Rayleigh 1911):

$$\begin{array}{c}\hfill v_{\mathrm{Shear}}(y)=|v_{{\mathrm{Shear}}_0}|\mathrm{erf}\left(\frac{y}{2\sqrt{\nu t}}\right).\end{array}$$(D.1)

Here, v_Shear0 is the initial shear velocity, y is the direction perpendicular to the interface, t is the time and ν is the kinematic viscosity. The more viscous the medium is, the faster the velocity gradient is smoothed out. This can be seen in the top panel of Fig. D.1, where in the weak cooling regime a higher viscosity leads to a smoother velocity profile within the same amount of time. The effect of turbulence can be seen in the "No visc" case, where the profile is not completely smooth. However, in the strong cooling regime all the profiles remain sharp after the same amount of time due to the effect of cooling.

Fig. D.1. vshear profiles after t = 1.25τKH color-coded by viscosity. Top: Results for Damix ≃ 2.3 × 10−3, where the velocity gradient is smoothed depending on the amount of viscosity. Bottom: Results for Damix ≃ 1.1 × 103, where cooling dominates over viscosity and keeps the profiles sharp.

Appendix E: Measurement of the smoothing distance, d

To quantify the distance, d, at which the velocity gradient has been smoothed due to the effect of viscosity, we focus on the hot medium (y > 0; see Fig. E.1). The reason for this is that, in a system with a constant dynamic viscosity (η), the high density of the cold medium leads to a lower diffusion coefficient (ν), therefore the overall diffusion of the system is dominated by the hot gas (see Sect. 5.1). The effect of viscosity is given by Eq. D.1, which tells us how much the velocity gradient has been smoothed after a time t for a given viscosity. By fitting this equation to our data (thick dashed lines), we have a smooth function for our system, avoiding spurious artifacts that can affect the results. Once we have fitted the function, we calculate the distance d from the center to the point where the velocity gradient has been smoothed 5% with respect to the initial value (vertical dashed lines).

Fig. E.1. vshear profiles at Damix ≃ 2.3 × 10−3 and t = 1.25τKH showing the fit of Eq. D.1 for the different viscosities (thick dashed lines). The distance d at which the shear profile has been smoothed 5% with respect to the original value is given by the vertical dashed lines.

Appendix F: Convergence test

To make sure that our results are robust against numerical resolution, we have performed a sample of simulations for the weak and strong cooling regimes Da_mix ≃ 2.3 × 10⁻³ and Da_mix ≃ 1.1 × 10³, respectively, varying the resolution: from (16 × 160 × 16) to (128 × 1280 × 128).

Figure F.1 shows the evolution of Q over time for the weak cooling regime (top panel) and for the strong cooling regime (bottom panel), for the nonviscous (blue) and 0.5η_Crit (red) cases at different resolutions. In all cases, the results shown in the paper hold: no significant difference between viscous and inviscid cases in the weak cooling regime, and a small increase (factor of ∼2) in the strong cooling regime. It is important to note two details:

• The runs with the highest resolution in the weak cooling regime lead to a slightly larger surface brightness compared to the runs with lower resolution. However, far from making it worse, an increase in Q means that the data points in Fig. 7 of the paper become closer to the expected value at that given Da.

• The viscous run with the highest resolution in the strong cooling regime leads to a slight increase in Q. We attribute this to numerical noise associated with the frame boost (see Sect. 3.2) that we found when increasing resolution, and not a lack of convergence (note that the rest of the runs are perfectly converged).

Another important note here is that, since we have changed the resolution, as a consequence, we have modified the numerical conduction. This means that, in this analysis, the Karlovitz number is expected to be different than the one given throughout the whole paper. Changing both the viscosity and conduction in a controlled manner, thus, assessing the impact of mixing and cooling in the full Ka range, is a logical next step but beyond the scope of the current study.

Fig. F.1. Evolution of the surface brightness due to cooling varying resolution from (16 × 160 × 16) to (128 × 1280 × 128) for the nonviscous (blue) and 0.5ηCrit (red) cases. Top: Weak cooling regime, Damix ≃ 2.3 × 10−3. Bottom: Strong cooling regime, Damix ≃ 1.1 × 103.

All Tables

All Figures

Fig. 1. Dependence of ηCrit on the overdensity, χ. Red dots show the numerical results. The blue line shows the theoretical expectation assuming that only the low density is relevant for kinematic viscosity. The red line shows that only the high density is relevant. The green line assumes an arithmetic mean of kinematic viscosity between the fluids. The black line shows a harmonic (weighted) mean. Note that we normalized the values of viscosity to the Spitzer value at 106 K.

In the text

	Fig. 2. Numerical solution for the values of Im(ϖ) of Eq. (26) as a function of ℳh for χ = 5, 10, 100. The dashed lines show the ℳh above which the KHI is stable.
In the text

	Fig. 3. Growth of the KHI with time for χ = 10, color-coded by ℳh. The solid lines indicate the simulations in which the KHI is still unstable and the dash-dotted lines where the instability is suppressed (values above ℳCrit), finding a ℳCrit = 1.65 ± 0.5.
In the text

	Fig. 4. Dependence of ηCrit on the Mach number, ℳh. The numerical results (red dots) deviate from the analytical incompressible ηCrit(ℳh) (Eq. (24); dashed black line). However, this is corrected by taking into account the suppression due to high ℳh (Eq. (30), dashed blue line). Note that we normalized the viscosities to the Spitzer value at 106 K.
In the text

	Fig. 5. Temperature slices for the different TRML simulations. From left to right: Non-viscous case, 5% of critical viscosity, 10% of critical viscosity and 50% of critical viscosity. Top: Very weak cooling with Damix ≃ 2.3 × 10−3. Bottom: Very strong cooling with Damix ≃ 1.1 × 103.
In the text

	Fig. 6. Evolution of the surface brightness due to cooling in the mixing layers. Left: Evolution of Q color-coded by viscosity in a very weak cooling regime with Damix ≃ 2.3 × 10−3. Right: Evolution of Q for different viscosities in a very strong cooling regime with Damix ≃ 1.1 × 103.
In the text

	Fig. 7. Surface brightness, Q, as a function of Damix color-coded by viscosity. The darker dashed line shows the best fit in each case in the weak cooling regime (Da < 1) with a slope of α = 1/2. The lighter dashed line shows the best fit in each case in the strong cooling regime (Da > 1) with a slope of α = 1/4.
In the text

	Fig. 8. Relation between Q, normalized by the numerical measured Damköhler number ($\stackrel{\sim }{\mathit{Da}}$) in each regime, and the turbulent velocity for different viscosities and color-coded by Damix. The dashed line shows the relation $Q/{\stackrel{\sim }{\mathit{Da}}}^{\alpha }\propto u\prime$ expected.
In the text

	Fig. 9. Dependence of the turbulent velocity on Damix for different amounts of viscosity and color-coded by the surface brightness. The dashed black line indicates the dependency of $u\prime \propto {{\mathrm{Da}}_{\mathrm{mix}}}^{1/10}$. The dashed blue lines show the turbulent velocity values of the adiabatic system with different viscosities.
In the text

	Fig. 10. Top: Turbulent velocity as a function of the velocity gradient width, color-coded by the surface brightness. Bottom: Surface brightness as a function of the velocity gradient width, color-coded by Damix.
In the text

	Fig. 11. Measured critical Da (Da∗) vs. the fraction of critical viscosity (ζ) for different sets of simulations. We included two more sets for completeness that are not shown in the previous plots. The dashed line indicates the expected one-to-one relation.
In the text

	Fig. 12. Turbulent velocity as a function of the amount of viscosity of the medium (Karlovitz number), color-coded by Damix.
In the text

	Fig. 13. Exponent of the dependence of the turbulent velocity on the amount of viscosity (Karlovitz number) calculated from Fig. 12 as a function of Damix.
In the text

	Fig. 14. Velocity profile of 1D simulations. Top: Only the cold gas (y < 0) is viscous. Bottom: Only the hot gas (y > 0) is viscous.
In the text

	Fig. 15. vshear profiles at Damix ≃ 1.1 × 102 and t = 1.25τKH including the run with 2ηCrit. The distance, d, at which the shear profile has been smoothed in the 2ηCrit case is large compared to the other cases.
In the text

Fig. A.1. Growth of the instability as a function of time for an overdensity of χ = 10, color-coded for different values of viscosity. The value of the critical viscosity is equal to ηCrit = 515 ± 5, where the error is given by the average between the viscous run with a positive and a negative slope. The solid lines indicate the simulations in which the KHI is still unstable and the dash-dotted line where the instability is suppressed (η > ηCrit).

In the text

	Fig. B.1. Mass conversion from hot to cold gas as a function of the total luminosity, color-coded by Da. The mass transfer is proportional to the total luminosity, indicating that the luminosity is due to enthalpy flux.
In the text

	Fig. C.1. Intermediate temperature histograms color-coded by viscosity. From top to bottom: Weak cooling (Damix ≃ 2.3 × 10−3), intermediate cooling (Damix ≃ 2.3), and strong cooling (Damix ≃ 1.1 × 103). From left to right: Histogram measured at t = 2.5τKH, t = 5τKH and t = 6.25τKH.
In the text

	Fig. D.1. vshear profiles after t = 1.25τKH color-coded by viscosity. Top: Results for Damix ≃ 2.3 × 10−3, where the velocity gradient is smoothed depending on the amount of viscosity. Bottom: Results for Damix ≃ 1.1 × 103, where cooling dominates over viscosity and keeps the profiles sharp.
In the text

	Fig. E.1. vshear profiles at Damix ≃ 2.3 × 10−3 and t = 1.25τKH showing the fit of Eq. D.1 for the different viscosities (thick dashed lines). The distance d at which the shear profile has been smoothed 5% with respect to the original value is given by the vertical dashed lines.
In the text

	Fig. F.1. Evolution of the surface brightness due to cooling varying resolution from (16 × 160 × 16) to (128 × 1280 × 128) for the nonviscous (blue) and 0.5ηCrit (red) cases. Top: Weak cooling regime, Damix ≃ 2.3 × 10−3. Bottom: Strong cooling regime, Damix ≃ 1.1 × 103.
In the text