The Anatomy of a Turbulent Radiative Mixing Layer: Insights from an Analytic Model with Turbulent Conduction and Viscosity

Turbulent Radiative Mixing Layers (TRMLs) form at the interface of cold, dense gas and hot, diffuse gas in motion with each other. TRMLs are ubiquitous in and around galaxies on a variety of scales, including galactic winds and the circumgalactic medium. They host the intermediate temperature gases that are efficient in radiative cooling, thus play a crucial role in controlling the cold gas supply, phase structure, and spectral features of galaxies. In this work, we develop an intuitive analytic 1.5 dimensional model for TRMLs that includes a simple parameterization of the effective turbulent conductivity and viscosity and a piece-wise power-law cooling curve. Our analytic model reproduces the mass flux, total cooling, and phase structure of 3D simulations of TRMLs at a fraction of the computational cost. It also reveals essential insights into the physics of TRMLs, particularly the importance of the viscous dissipation of relative kinetic energy in balancing radiative cooling as the shear Mach number approaches unity. This dissipation takes place both in the intermediate temperature phase, which reduces the enthalpy flux from the hot phase, and in the cold phase, which enhances radiative cooling. Additionally, our model provides a fast and easy way of computing the column density and surface brightness of TRMLs, which can be directly linked to observations.


INTRODUCTION
The interaction and mixing of a hot wind and cold gas clouds is a ubiquitous situation in many astrophysical systems. In particular, both simulations (e.g. Kim & Ostriker 2018;Schneider et al. 2020) and observations (e.g. Strickland & Heckman 2009;Rupke 2018) suggest that astrophysical gases in and around galaxies are highly multiphase. That is, cold and hot gases coexist on a wide range of scales. These cold and hot phases are often in pressure equilibrium and move relative to each other. At the interface of the two phases Kelvin-Helmholtz instability (KHI) drives turbulent mixing which populates the intermediate temperatures that radiate efficiently, thus forming turbulent radiative mixing layers (TRMLs).
TRMLs are essential in controlling the phase structure and evolution of galaxies on a variety of scales, including the interstellar medium Audit, E. & Hennebelle, P. 2010) predict that cold clouds should be shredded apart prior to being accelerated to high velocities that they are observed to have in galactic winds and the circumgalactic medium (e.g., Zhang et al. 2017). Intermediate temperature gases in TRMLs consist of partially ionized atoms that are most conducive to radiative cooling, which means they convert hot gas to cold gas rapidly, allowing the original cold cloud to survive or even grow when it gets entrained in the hot wind. Later works confirmed the conclusion of cloud growth in the presence of radiative cooling when t cool,mix / t cc 1, where t cool,mix is the characteristic cooling time of the mixed gas 1 , and t cc is the cloud crushing time (e.g., Kanjilal et al. 2021;Abruzzo et al. 2021;Gronke et al. 2021, see Li et al. 2020Sparre et al. 2020 for a somewhat different interpretation predicated on the central importance of the hot phase cooling time).
Thus, TRMLs dictate the fate of cold clouds and play a crucial role in regulating the amount of cold gas that is available for star formation and black hole growth in galaxies. Turbulent Radiative Mixing Layers are also essential to observations because they host the intermediate temperature gases that contribute a significant amount of the spectral features (see e.g. Gronke & Oh 2020 on how a significant amount of emission from a cloud comes from TRMLs) commonly seen in both emission (e.g., Hayes et al. 2016;Reichardt Chu et al. 2022) and absorption (e.g., Werk et al. 2014;Qu et al. 2022).
There has been significant work focusing on the interface to study TRMLs specifically. Begelman & Fabian (1990) were amongst the earliest, and constructed a model where the advected enthalpy from the hot phase is balanced by cooling and discussed the rich emission and absorption features that are produced by the mixing layers. More recently, Kwak and collaborators (Kwak & Shelton 2010;Kwak et al. 2011;Henley et al. 2012;Kwak et al. 2015) performed 2D simulations of TRMLs. They compared their results with observables, including line ratios and column densities. Later on, Ji et al. (2019) ran 3D simulations of TRMLs to compare with previous analytic theories, finding surprisingly low inflow and turbulent velocities, as well as different scalings for the mixing layer width and surface brightness with cooling time. Recent TRML models based on high-resolution, three-dimensional simulations (Fielding et al. 2020;) have further elucidated the physical processes responsible for regulating the rate at which mass is exchanged through the fractal surface separating the hot and cold phases. Besides three-dimensional simulations on a variety of scales that have yielded essential insights into the physics of TRMLs, another valuable approach to study TRMLs is through analytic, one-dimensional models. McKee & Cowie (1977) were among the first to develop analytic solutions for the gas transfer between a spherical cloud and the surrounding hot gas. They modeled thermal conductivity as a powerlaw of temperature (taken from Massey et al. (1975) and Spitzer (1962)) and identified a critical radius that determines whether the cold cloud evaporates or condenses.  and Inoue et al. (2006) both analyzed the structure and stability of TRMLs (which they called the transition layer that connects the warm and cold neutral medium) by finding one-dimensional, steady-state solutions to the fluid equations in the reference frame where the transition layer is stationary. Binette et al. (2009) also developed a one-dimensional analytic model for TRML to understand its impact on observational signatures like line profile shapes. Notably, Binette et al. (2009) adopted a mixing length approach similar to what we will discuss in § 2.3 and uses a constant turbulent viscosity. More recently,  and  argued that understanding the competition between thermal conduction and cooling is a crucial step towards capturing the detailed phase structure of TRMLs, including the temperature distribution and column density. They demonstrated that despite its complex fractal structure, TRMLs can actually be described by a simple 1D conductive-cooling front model that quantitatively reproduces 3D simulation results including column densities and line ratios. These works suggest that studying TRMLs through 1D models is a promising approach that has the potential of reducing computational cost and helping us develop physical intuition.
In this work, we develop an analytic, 1.5-dimensional 2 model for Turbulent Radiative Mixing Layers (TRMLs). The crucial difference between our model and the models mentioned above is that we now explicitly include the impact of turbulence driven by the shear flow (which is likely to dominate over the physical conduction and viscosity in most regimes). We do this by introducing a simple parameterization of the effective turbulent conductivity and viscosity that is proportional to the shear velocity gradient. We also adopt a piece-wise power-law cooling curve that takes into account of both collisional ionization and photo-ionization equilibrium. Our model allows us to better quantify the role viscosity plays in TRMLs. It also reproduces the mass, energy, momentum transfer and shape of the phase distributions from 3D simulations (Fielding et al. 2020). In § 2, we recast the steady state fluid equations into our desired form and introduce the cooling curve and effective turbulent conductivity and viscosity in our model, as well as the simulations we compare against. In § 3, we present results of our model, including: the relationship between the mass fluxṁ and total cooling Q cool derived from energy conservation, the relative importance of energy sources and sinks in TRMLs as a function of parameters of our model, a demonstration that we match key features of our analytic model with 3D simulations, and a determination of the turbulent Prandtl number Pr . In § 4, we compare with previous works and discuss observational implications.

METHODS
We consider a 1.5 dimensional model for a turbulent radiative mixing layer. The coordinate system we will adopt has the mass fluxṁ between the phases in thê z direction and the shear motion between the phases in thex direction. We will adopt a convenient notation in which the pressure P = ρT , where we have absorbed the constant factor k B / (µm p ) into our definition of T (with this convention T has units of velocity squared, and can be thought of as the isothermal sound speed squared). Since this is a 1.5 dimensional model the derivatives perpendicular toẑ are all zero.
The equations of mass, momentum, and energy conservation are ∂ρ ∂t + ∇ · (ρv) = 0 (1) where the total internal energy is given by E = 1/2ρv 2 +P/(γ−1), Π = µ (∇v+∇v −2/3I∇ · v) is the viscous stress tensor, Q = −κ∇T is the conductive heat flux, andĖ cool is the cooling rate. In reality, the conductivity κ and viscosity µ of a fluid are mediated by the discrete collisions of the particles that comprise the fluid. In practice, in addition to the physical κ and µ, turbulence can effectively act as an additional form of "eddy" conductivity and viscosity, as has long been studied (e.g. Frisch 1995). In this paper, we assume that the physical κ and µ values are negligible and that the properties of TRMLs are dominated by the turbulent κ and µ. We will discuss how we model µ and κ in more detail in § 2.3. Throughout this paper, we use 3D simulations of TRMLs as a point of comparison to our analytic model. The 3D simulations are described in detail in § 2.4. For now, we utilize two pieces of information from the 3D simulations to facilitate our derivation of the fluid equations. First, we show an annotated 2D projection from the 3D simulations in Fig. 1 that clearly illustrates the basic geometry of the system we are dealing with. Second, we make use of a fundamental property of rapidly cooling TRMLs that has been uncovered in 3D simulations. Namely, over time the horizontally averaged shear velocity and momentum profiles continually expand due to the effective turbulent viscosity. The horizontally averaged temperature and inflow velocity profiles, however, do not. Instead these profiles quickly reach an apparent equilibrium configuration and maintain it indefinitely. An example of this behavior is shown in Fig. 2 for the same 3D simulation shown in Fig. 1. Physically, this can be understood by comparing the timescales over which these properties change. The shear profile spreads out on a timescale proportional to the shear time t sh = L/v rel , where L is the characteristic lengthscale in the direction of the shear flow. The temperature and inflow velocity, on the other hand, evolve on a timescale proportional to the cooling time t cool . Therefore, in the rapid cooling limit in which t cool t sh , the temperature and inflow velocity profiles evolve much more rapidly than the shear velocity and, thus, can, to good approximation, be considered to have ample time to reach a quasi-steady state equilibrium configuration. Note that in the slowly cooling limit (t cool t sh ), which we do not consider here, this assumption breaks down and all of the quantities will evolve on similar timescales.
In our model we make use of this quasi-steady state nature of all horizontally averaged quantities except the shear velocity to explicitly assume that the time derivative of ρ, v z , and T is zero, and only the time derivative of v x is non-zero. Using this we can now expand the 4 equations that fully describe our system. Figure 1. Example of the structure of a turbulent radiative mixing layer (TRML) three dimensional simulation similar to what was published by Fielding et al. (2020) but using the cooling curve we have adopted here. Going from left to right the panels show the mass weighted average temperature T , the root mean square pressure Prms, the volume weighted density ρ, the average cooling time t cool (defined to be the ratio of integral of the thermal energy divided by the integral of the radiative cooling along theŷ direction), the mass weighted averagex-velocity vx, the root mean squareŷ-velocity vy,rms, and the mass weighted averagê z-velocity vz. Overplotted in grey is the horizontally averaged profile of each of these quantities. This demonstrates the general configuration in which the hot phase is flowing in toward the cold phase and undergoes a rapid transition from hot to cold, which is where the vast majority of the cooling occurs. This is also the locus of the turbulence, which is traced out by vy,rms. The shear velocity vx has a thicker transition region.
Which then gives uṡ where κ is the conductivity, µ is the viscosity given by µ = Pr κ, andĖ cool is the radiative cooling function. At this point we are at a crossroads for how to solve the above system of equations. In principle, one could solve the system numerically as it is, which involves dealing with both spatial and temporal derivatives. Instead, we rely on the fact that all the time derivatives are applied on v x , which is slowly varying. We exploit this feature to solve for a snapshot of the system in time by prescribing a profile for v x by hand. The details of how we do that will be described in § 2.5. For now, we rearrange the above system to solve for expressions for the second derivatives of v z and T . The general form for these expressions for some arbitrary choice of κ and a constant Pr is Since with our assumptions T and v z are now only functions of z, we change all partial derivatives to full derivatives to emphasize that we are solving for a time-independent slice of the solution. Note that for the dv x / dz term that appears in Eq. 13, we plug in our prescribed v x profile, which is informed by results from 3D simulations. This precludes the need of simultaneously solving a third v x equation that involves time derivatives.
With the system of Eq. 12 and Eq. 13, we can easily plug in different choices for κ and Pr , and then solve the resulting equations numerically.
It is worth clarifying that the time-independence of Eq. 12, Eq. 13, and our prescribed v x profile does not mean that the evolution of TRMLs is time-independent. Solving Eq. 12 and Eq. 13 only gives us a "snapshot" of the TRML in time. To solve for the time evolution, one 0.0 0.5 Figure 2. Time evolution of the horizontally averaged profiles of temperature (top) and shear velocity (bottom) weighted by mass in the same 3D turbulent radiative mixing layer (TRML) shown in Fig. 1. At each time the profiles are aligned to the z location where T = T hot /2, which for this simulation corresponds to boosting into a reference frame shifted by a small z velocity equal to −v rel /200. While the shear velocity profiles continues to broaden in time the temperature profile quickly reaches an apparent equilibrium configuration and maintains it throughout the long duration of this simulation. would need to find ∂v x / ∂t using which comes from solving for the second derivative of v x in Eq. 9, Eq. 10, and Eq. 11. Using Eq. 14 and an appropriate time step, one can time-evolve the original prescribed v x profile and plug the new v x profile into the system of Eq. 12 and Eq. 13 to solve for a new snapshot of the TRML at a later time. This process is akin to the approach that has been adopted for many other systems with multi-scale processes, such as stellar structure (e.g., Kippenhahn & Weigert 1990). Turning our attention back to the system we are solving (Eq. 12 and Eq. 13), it is instructive to rearrange Eq. 13 into a more intuitive form: where we have labeled each term by its physical meaning using the same colors that will be used in subsequent figures. The work term comes from the following relation If we integrate this from hot to cold then we geṫ where Q cool = Ė cool dz (18) Recall that we have absorbed the constants into our definition of T , so γ∆T = c 2 s,hot − c 2 s,cold = ∆c 2 s . It is useful to think of Eq. 17 as a statement of energy sources and sinks. It tells us that the advected enthalpy from the hot phase (LHS of Eq. 17) receives three contributions (RHS of Eq. 17): radiated away by cooling, gained through viscous dissipation of the relative kinetic energy between the phases, or work done by the pressure gradient.

Conventions and Unit Conversion
In the remainder of this paper, many physical quantities will be expressed as dimensionless ratios. This keeps our approach as general as possible, but it can be useful to provide physical conversions to guide intuition and comparison to observations. In this section, we explicitly state the conventions we are using and discuss how to convert between code units and physical units.

Dimensionless Numbers
The structure of a TRML is strongly affected by three dimensionless numbers that characterize the physics of the system, which are where v rel is the relative velocity between the hot and cold phases, c s,hot is the hot phase sound speed, γ is the adiabatic index that we take to be 5/3 throughout, t cool,min is the minimum cooling time that occurs at intermediate temperatures, L 0 , t 0 , and v 0 are the characteristic length, time, and velocity of the system, µ is the mean molecular weight and is equal to 0.62 in an ionized medium with solar metallicity, and m p is the mass of a proton. We set these characteristic scales by setting (i) the hot phase pressure P hot , (ii) the hot phase temperature T hot , and (iii) the minimum cooling time t cool,min . With these the velocity scale is v 0 = P hot / ρ hot = k B T hot / µm p , which differs from c s,hot by √ γ.
As we will see later in § 2.2, the dimensionless number τ serves as a normalization to the cooling rate and the cooling time. In 3D simulations, this normalization is set relative to the shear time t sh = L/ v rel , where L is the shear layer length, which is well-defined in 3D TRML simulations. Thus by setting a physical scale for L, and the ratio t cool /t sh , the physical cooling rate is determined. However, here in our 1.5D model we do not have a transverse lengthscale, which precludes us from implementing the same setup. We, therefore, take the opposite approach and set the physical timescale which in turn sets the physical lengthscale. We set the physical timescale using the following relation: The numerical value in Eq. 24 comes from using an approximate fit to the cooling rates from Wiersma et al. (2009) assuming solar metallicity. In reality there should be a metallicity dependence on the minimum cooling time, however, because hydrogen is the primary coolant at the temperatures where the cooling time is minimized (∼ 2 × 10 4 K ) the metallicity dependence is very weak, so we leave it out for sake of simplicity. Thus, at a fixed pressure, adjusting τ implies adjusting the characteristic time t 0 , which is related to the characteristic length of the system by where v 0 is the characteristic velocity, and t 0 is the characteristic timescale that comes from Eq. 24. In other words, by adjusting τ we are effectively adjusting the lengthscale of the system. This characteristic length L 0 allows us to convert other lengthscales in the problem into physical units. For example, the thickness of the mixing layer, which we will later denote as h, is given by h × L 0 in physical units.
Whenever possible, we also express physical quantities in a dimensionless form. For example, we speak in terms of the dimensionless T / T hot instead of T . We can trivially make the conversion to dimension-full quantities by multiplying through by the hot phase temperature T hot , which can be chosen freely based on the physical scenario at hand. The exact same logic holds for pressure (which we normalize by the hot phase pressure P hot ) and density (which can be calculated from temperature and pressure using the ideal gas law).

The Phase Distributions
To understand the detailed phase structure of TRMLs, it is instructive to examine the distribution of mass and the amount of cooling as a function of temperature. We call these distributions the phase distributions. Here we discuss the definition and calculation of the cooling flux and column density distribution of TRMLs.
The cooling flux distribution, or dF cool / dlogT, tells us how the cooling flux is distributed in log temperature space. The cooling flux distribution can be calculated by whereĖ cool (z) is rate of energy loss from radiative cooling. Similarly, the temperature distribution of mass per unit area, or more conveniently column density N , is given by where m H is the mass of a hydrogen atom. The cooling flux and column density distributions are also normalized to be dimensionless as follows:

Cooling Curve
Since we are interested in general properties and behavior of TRML, we adopt a simplified cooling rate that depends only on temperature; the overall shape is intended to be similar to more realistic curves but is both simpler and customizable. In particular, we use a cooling curve that is a piece-wise power-law. The powerlaw slopes are β cold = −2 if T ≤ T peak and β hot = 1 if T > T peak . Additionally, we include a heating term such that the cooling curve has two equilibria at both the cold and the hot phase temperatures. In realistic systems P/k B = 10 3 K cm −3 Z = Z β cold = −2 β hot = 1 Figure 3. Comparison of the realistic cooling time profile with our piece-wise power-law cooling time prescription. The realistic cooling time profile (in black) assumes P / k b = 10 3 Kcm −3 and solar metallicity. The piece-wise powerlaw profile (in blue) has β cold = −2 and β hot = 1, where β cold and β hot are the power-law slopes at T ≤ T peak and T > T peak , respectively. We discuss our piece-wise power-law prescription in more details in § 2.2. As shown in the figure, this prescription closely resembles the realistic profile and has the advantage of being simpler and easily customizable.
the hot phase does not have a stable equilibrium and is instead very slowly cooling. We adopt this artificial hot phase equilibrium to improve numerical performance but in practice it has no appreciable impact on our results. We discuss the cooling curve we use in more quantitative details in appendix. In Fig. 3, we compare the cooling time profile obtained from our piece-wise power-law cooling curve with the realistic cooling time profile. For most of the remainder of this paper, we are going to use our piece-wise powerlaw cooling curve as a proxy to the realistic cooling curve. However, in § 4.3 and § 4.4 where we calculate the ion column density and emission line surface brightness of TRMLs from results of our analytic solutions, we use the realistic cooling curve to obtain more accurate results. Smagorinsky (1963) is one of the first to develop a framework for modeling viscosity from turbulence in eddies. In his pioneering subgrid-scale model for large eddy simulations, Smagorinsky (1963) defined what he called the eddy viscosity as

Conductivity & Viscosity
where is the lengthscale of the large eddies,S ij is the strain rate tensor, and S ijSij gives the velocity scale corresponding to the lengthscale . Prandtl (1925) also took the large-scale eddy viscosity (comparable to what we will later refer to as the effective turbulent viscosity in our analytic model) into account in his mixing length theory. He recognized this eddy viscosity as characterizing the transport and dissipation of momentum and energy due to turbulence on small scales. Ultimately, the physics of turbulence and viscosity determines the distance a gas cloud can travel with its original temperature and density before blending into the surrounding. Prandtl (1925) called this distance "the mixing length".
Mixing length theory (Prandtl 1925) has been used in previous works on TRMLs. For instance,  developed a 1D mixing length model of TRMLs (with a constant or temperature-dependent conductivity and no viscosity) and computed the resulting mean density and temperature profiles that matches well with 3D simulation results.
Motivated by Smagorinsky (1963) and Prandtl (1925), we introduce a model for effective turbulent conductivity and viscosity. The conductivity is given by where ρ and h are the density and thickness of the mixing layer, v x is the shear velocity, and κ 0 is a small additive constant used to prevent singularities when dvx dz → 0. The κ 0 term is somewhat analogous to the numerical dissipation in the 3D simulations and we find that as long as it is small compared to the first term our results are insensitive to the exact choice (we adopt κ 0 = 10 −6 throughout). We include the absolute value on dv x / dz because κ must always be greater than 0. Note that in the definition of κ we introduce the constant, dimensionless prefactor f ν , which can be used to adjust the magnitude of κ. Another useful way to understand the utility of f ν is to examine the combination f ν h 2 , which has dimensions of length squared. By defining f ν we are effectively introducing a new lengthscale given by f ν h 2 1/2 = f 1/2 ν h that is equivalent to in the Smagorinsky/mixing length formalism. Since in our 1.5D model we do not have the freedom to set a transverse lengthscale as in the 3D simulation, this new lengthscale defined by f ν is introduced to compensate for the lost of that degree of freedom.
We assume that the effective turbulent viscosity has the same functional form as the conductivity, and is thus defined as where Pr is the effective turbulent Prandtl number. We note that in previous works on analytic, 1D models of TRMLs (e.g., , the conductivity κ is frequently modelled as either a constant or a function of temperature. This temperature dependence is either taken to follow the Spitzer conductivity scaling (κ ∝ T 5/2 ) (Spitzer 1962), or is given an arbitrary power law scaling that is calibrated to match the results of 3D simulations. In these models the viscosity µ is not taken into account. In this way, our setup is fundamentally different from these previous works as here we compute the effective turbulent conductivity and viscosity based on the shear velocity profile (see Section 2.5), which enables us to connect more closely to the underlying physical processes.

Description of Benchmark 3D Simulations
In order to better motivate the choices inherent in our model and to calibrate the parameters we have introduced, f ν and Pr , we compare to a series of 3D simulations.
These simulations are similar to what were presented in Fielding et al. (2020) with the primary difference being that they use the cooling curve as we have adopted here. A full description of these simulations as well as a detailed analysis will be presented in an upcoming work. We now briefly summarize these simulations and their analysis.
The simulations were performed with the athena++ code framework (Stone et al. 2019), which solves the standard fluid equations (see equations 1-3). The simulations were performed with no explicit conduction or viscosity. As mentioned, the simulations include the exact same cooling and heating as given in Eq. A1.
We now briefly describe the simulation setup and the suite of parameter variations, a more comprehensive explanation will be included in a forthcoming companion paper. The simulations are initialized in pressure equilibrium between the hot and cold phase with a smooth temperature transition spread out over 4 cells. The hot and cold phases are initialized so that they are moving relative to each other in thex direction with a velocity given by v rel . As with the temperature, the shear velocity transition is resolved by 4 cells. All three components of the velocity at the interface between the hot and cold phase is perturbed in the initial conditions. The root-mean-square of the velocity perturbation at the interface is set to v rel /20 and falls off exponentially with distance from the interface with a scale length of 16 cells. The perturbations are generated in Fourier space and have equal power on wave numbers 1 ≤ kL/2π ≤ 16.
As the simulation progresses in time the shear velocity/momentum transition layer spreads.
The temperature, density, and inflow velocity respond quickly to the gradual spreading of the shear velocity and maintain nearly the same horizontally averaged profile. This is shown explicitly in Fig. 2. This motivates our assumption that these quantities can be taken to be in steady state while the shear velocity cannot. We must, therefore, assume some shear velocity profile in order to solve our 1.5D model. We discuss our approach to doing so in the next subsection. We need to first, however, decide when to measure the simulation properties since they evolve with time. In order to standardize the measurements across the range of parameters all 3D simulation properties are measured when the thickness of the mixing layer is equal to the transverse width of the box, i.e. h = L. 3 We measure the thickness h by calculating the distance between the places where thex andŷ mass-weighted average shear velocity v x is within 5 percent of the initial values (i.e. between 0.05 and 0.95 v rel ). We tested the sensitivity to the exact definition and found only minor quantitative differences. Because this definition is somewhat imprecise in practice we take the time average of when h is within 5 percent of L. Again the exact choice has no substantive impact on our results, but helps avoid fluke temporal fluctuations. When comparing to simulation properties in later sections we use error bars to denote these temporal fluctuations about the mean.

Taking a Prescribed Shear Velocity (v x ) Profile
from the 3D Simulation Instead of solving Eq. 14, which contains a time derivative of v x , together with Eq. 12 and Eq. 13, we resort to the 3D simulations to obtain a prescribed v x profile and solve only for Eq. 12 and Eq. 13.
We explore two different ways of extracting the v x profile from 3D simulations.
First, the most straightforward way is to directly take the mass weighted horizontally averaged v x profile from the 3D simulation. The physics of the system at hand dictates that the v x profile we use in our analytic model has to be continuous and perfectly flat at both the hot and cold phases ( dv x / dz = 0 at both ends of the solution.) Additionally, it is important to remember that the ultimate purpose of extracting the v x profile is to feed values of dv x / dz into our definition of the effective turbulent conductivity and viscosity in § 2.3, which means whatever v x profile we come up with must also be differentiable. Due to these considerations, we choose to fit the 3D simulation v x profile with a cosine 4 , as shown in the top panel of Fig. 4. We find the cosine fit convenient because it satisfies all the requirements Figure 4. Comparison of the vx and ρvx profiles from 3D simulation and the cosine fits that we use in our analytic model. The solid curve on the top panel is the mass weighted x-velocity from the 3D simulation, and the solid curve on the bottom panel is the volume weighted density times the mass weighted x-velocity from the 3D simulation, which gives us a x-momentum profile. The dotted curves on the top and bottom panels are the cosine fits to the vx and ρvx profiles, respectively. We discussed these cosine fits in more details in § 2.6 and § 2.7. mentioned above, has a simple functional form, and provides a reasonable fit all at the same time.
Alternatively, we also tried extracting the xmomentum profile ρv x from 3D simulations and fit it with a cosine. Then, the v x profile we need can be obtained by dividing the x-momentum profile by ρ, oṙ m/ v z . Although this approach is more algebraically complex, one immediate benefit is that the cosine fit is much more suitable for the x-momentum profile, as shown in the bottom panel of Fig. 4. As we will see later, the x-momentum profile setup also exhibits some characteristics we see from the 3D simulations but not in the v x profile setup, which makes these two approaches complementary.
Next, we rigorously define the v x profile and ρv x profile we are using and discuss the algebraic steps that need to be taken to simplify Eq. 12 and Eq. 13 in each of the two setups. In § 3, we begin our investigation with the v x profile prescription and later compare results from the 3D simulations with results from both the v x and the ρv x profiles.

v x Profile as a Cosine
Here where h is the thickness of the mixing layer, and v rel = |v Using Eq. 31 and ρ =ṁ/ v z , we can express dκ/ dz as: Plugging Eq. 31 and Eq. 36 into Eq. 12 and Eq. 13 gives us a system we can numerically integrate.

ρv x Profile as a Cosine
Next, let's look at prescribing the x-momentum ρv x as a cosine profile. Before defining the profile, we first note a crucial constraint due to the expression of the shear velocity (v x ) gradient. The shear velocity gradient is given by where p x = ρv x is the x-momentum profile (which we will later define as a cosine). The constraint here is that the two terms in the expression for dv x / dz must have the same signs. Otherwise, dv x / dz would evaluate to 0 somewhere in the solution, causing a singularity and thus difficulty for the integration to proceed. Since we know for sure that p x > 0, dv z / dz < 0, and v z > 0 throughout the solution, this implies that dp x / dz must be less than 0. With this in mind, we define p x (z) as where ρ hot = P hot / T hot is a known boundary condition at the hot end (see Appendix B for more details). Under this setup, we have Note that the last term of Eq. 40 has a second derivative of v z in it, which means after we plug Eq. 39 and Eq. 40 into Eq. 12 and Eq. 13, we still need to rearrange and solve for the second derivatives in order to obtain a system that is numerically integrable. After some algebra, we get

Carrying Out the Numerical Integration
After prescribing the v x or the x-momentum profile, we are left to numerically integrate two coupled differential equations to solve for the z-velocity (v z ) and temperature profiles. We set appropriate initial conditions and use scipy's solve-ivp integrator to perform the numerical integration.
An additional constraint of this problem is that the temperature profile needs to be flat ( dT / dz = 0) at both the hot and cold phases by definition. For any given choice of parameters, there's only one value of the mass fluxṁ that allows for that (as shown in Fig. 5). We call that the eigenvalue ofṁ. To find this eigenvalue, we bisect the parameter space ofṁ and use scipy's root finder optimize.root-scalar to pinpoint the eigenvalue.
In Appendix B, we provide more details on the numerical integration procedure.

RESULTS
In this section, we present results from our analytic model of TRMLs. In § 3.1, we examine a fiducial solution of TRML produced by our analytic model. We discuss how changing the relative shear velocity changes the phase structure of TRMLs to help the . There is only one value of the mass fluxṁ that allows the temperature profile to reach the cold phase temperature smoothly ( dT / dz = 0 at the cold phase). We call it the eigenvalue ofṁ (ṁeigen). Here we plot 9 temperature profiles, varying the value of the mass fluxṁ by ±10% aroundṁeigen. The T profile corresponding toṁeigen is shown in black. Whenṁ >ṁeigen, the T profile shoots down to the cold phase with a steep slope, with dT / dz < 0 at the end of the solution (the red group of solutions in the figure); whenṁ <ṁeigen, the T profile barely reaches the cold phase before bouncing back up, with dT / dz > 0 at the end of the solution (the blue group of solutions in the figure). The different signs of dT / dz at the end of the solution in these two cases allow us to perform a bisection iṅ m to pinpointṁeigen. Our bisection method resolvesṁeigen to an accuracy of 1 part in 10 15 .
readers build intuition on the physical processes at play. In § 3.2 and § 3.3, we discuss the energy budgeting in TRMLs, identify the energy conservation criterion H visc =ṁv 2 rel 2, and demonstrate how this criterion constrains the choices of model specific parameters Pr and f ν . In § 3.4 and § 3.5, we examine how the relative importance of the energy sources and sinks in TRMLs change as a function of the relative shear velocity for the two variations of our analytic model (cosine v x profile and cosine ρv x profile). Up until this point, all the physical insights are obtained purely through our analytic model. In § 3.6, § 3.7, and § 3.8, we discuss how to match the results from our analytic model with the 3D simulations, including matchingṁ, Q cool , and the shapes of the phase distributions.

Structure of TRMLs (with the Cosine v x Profile)
We first turn our attention to analyzing the detailed structure of solutions we obtained with the cosine v x profile prescription. We argue that the relative shear velocity (M rel ) between the hot and cold phase affects the solutions the most in our model. Varying the choice of M rel has two noteworthy effects. First, the x-component of viscous heating is proportional to ( dv x / dz) 2 , which means larger M rel leads to more viscous heating. Second, as defined in Eq. 31, the conductivity κ is proportional to dv x / dz, so larger M rel also implies higher conductivity.
In what follows, we first analyze a solution with a fiducial choice of M rel to identify some key characteristics in the structure of the mixing layer. Next, we compare three solutions with low, medium and high values of M rel to understand how the previously identified characteristics change across different M rel .
3.1.1. A TRML with a Fiducial Choice of the Relative Shear Velocity M rel In this section, we analyze the structure a mixing layer with M rel = 0.75.
The first step towards understanding the mixing layer is to look at the behavior of the key physical quantities that are involved. In our case, we are most interested in the temperature (T ), velocity (v z ), shear velocity (v x ), and pressure (P ) within the mixing layer. The top left panel of Fig. 6 shows these 4 quantities plotted against the vertical distance into the mixing layer. To begin with, we note that the shear velocity (v x ) profile is prescribed to be a cosine function, as discussed in § 2.5. The T and v x profiles change separately. The T profile drops to the cold end at as early as z ∼ 0.2 and hovers just above the cold phase temperature when the v x profile just starts to change at about the same point in space. In other words, this solution cools first (in position space), then mixes, meaning that there is significant change in v x within the cold phase. The pressure profile provides important insights into the physics that is driving the mixing layer. Naively, one would expect that the rapid cooling in the mixing layer leads to a region of low pressure at intermediate temperatures, which creates a pressure gradient that accelerates the inflow of hot gas. However, the top left panel of Fig. 6 shows that this is clearly not the case. This suggests that it is turbulent mixing (which is in the form of effective turbulent viscosity and conductivity in our model), not pressure gradient, that controls the mass inflow rate and thus the radiative cooling rate. We further note that the hot phase pressure is slightly higher than that of the cold phase. This difference is accounted for by the additional ram pressure at the hot phase.
To better understand the physics that is controlling the behavior of the mixing layer, it is instructive to return to the fluid equations and deconstruct the advected enthalpy from the hot phase, or the temperature gradient, dT / dz, into its constituent terms. This can be done by solving for dT / dz in Eq. 13. dT / dz can be decomposed into terms that represent radiative cooling, conduction, viscous heating, and work done by the pressure gradient. By examining the relative importance of these terms, we can better understand the underlying physics that controls the mixing layer.
The two bottom panels of Fig. 6 shows the "terms decomposition" in position (z) and temperature space. There are three parts to the solution. Starting from the hot end, the mixing layer is initially dominated by conductive cooling. At intermediate temperatures, radiative cooling becomes significant and is balanced by conductive heating. This balance holds true until we reach the cold end, where the x-component of viscous heating takes over from conductive heating to balance the radiative cooling. This final part of the solution takes up a tiny fraction in temperature space but a significant portion of the position space. This is consistent with our previous observation that the T and v x profiles change separately.
The top right panel of Fig. 6 shows the column density and cooling flux distributions of the mixing layer (see § 2.1.2 for more details). Corresponding to the strong viscous heating zone near the cold end, we see a spike in the cooling flux distribution at low temperatures. This is a feature that is commonly seen in 3D simulations, and our analytic model is able to give us physical intuition towards the origin of this spike: radiating away the viscously heated material in the cold phase.

Comparing Three TRMLs with different Choices of the Relative Shear Velocity M rel
Now that we are familiar with the key characteristics of a fiducial mixing layer, we expand our analysis to include three mixing layers with different choices of the relative shear velocity M rel (including a supersonic M rel ) and understand how the key characteristics change across these mixing layers.
In Fig. 7, we present the temperature, pressure, velocity profiles (top panels) and the terms decomposition of dT / dz in position space (bottom panels) of three mixing layers with M rel = 0.25, 0.75, and 1.25. The most notable difference is that the magnitude of the x-component of viscous heating increases together with M rel . This is expected because the x-component of viscous heating is proportional to ( dv x / dz) 2 , and with a fixed thickness of the mixing layer and higher M rel , dv x / dz is bound to increase in magnitude. On the other hand, most of the qualitative features of the mixing layer remain the same regardless of the choice of M rel . All three mixing layers have a mixing zone where conduction and radiative cooling drives most of the temperature drop, and a viscous heating zone at the cold phase where the majority of the viscous heating takes place.
The left column of Fig. 8 shows the cooling flux and column density distributions of the three mixing layers we just examined. As M rel increases, the spike at the cold phase of the cooling flux distribution increases in height, corresponding to the increased magnitude of the x-component of viscous heating. Furthermore, these spikes also gradually widen and shift towards higher temperature, indicating that the edges of the cooling and mixing zones in the mixing layers are starting to blend into each other as M rel increases. For the sake of completeness, we also show how the the cooling flux   and column density distributions change as we vary τ and χ in the middle and right column of Fig. 8. In general, lower values of τ and χ correspond to a higher and narrower spike in the cooling flux distribution.
As the mixing zone at the cold phase widens and undergoes more viscous heating with increasing M rel , its mass content also increases, leading to the cold phase spikes in the column density distribution. Another interesting point to note is that the column density distributions clearly do not have the same shape as the t cool curve as is predicted in pure cooling flows (e.g. Fabian 1994). Above the temperature where the cooling curve peaks (T peak ), the power law slope of the column density distribution remains negative and gradually levels off, but certainly does not turn positive, as does the t cool curve.
3.2. Bernoulli Flux and Requiring H visc =ṁv 2 rel 2 In § 3.1 we looked at the detailed structure of the mixing layers as a function of position and temperature   by examining the relative importance of the terms in the fluid equations. Alternatively, it is also possible to treat the entire mixing layer as a whole. After deriving the equations we numerically solved at the beginning of § 2, we integrate from hot to cold to obtain a relationship between mass fluxṁ and total cooling Q cool given by Eq. 17. This is a statement of the energy sources and sinks in the entire mixing layer, and it helps us in understanding the energy budget of the system and allows us to explore how the importance of these energy sources and sinks change as parameters of our model are varied.
Besides the relationship between energy sources and sinks given by Eq. 17, it is also crucial to enforce energy conservation on the mixing layer. The enthalpy flux (due to temperature difference between the hot and cold phase) and relative kinetic energy flux (due to shear velocity difference between the hot and cold phase) that is advected into the mixing layer must be radiated away by cooling to ensure energy conservation and prevent energy build-up within any part of the mixing layer. (We note that the kinetic energy influx must be first dissipated into heat, and then radiated away by cooling.) In other words, Q cool as defined in Eq. 18 must be balanced by what we define to be the "Bernoulli flux", given byṁ ∆c 2 s (γ − 1) + v 2 rel 2 (Fielding & Bryan 2022). Another way to think about this is that the two phases are experiencing a perfectly inelastic collision through the mixing layer, and the relative kinetic energy between the phases must be dissipated into heat and radiated away in the mixing layer. Combining the Bernoulli flux equation with Eq. 17, we demand that H visc +W =ṁv 2 rel 2. In practice, W is several orders of magnitude smaller than H visc , which means our energy conservation condition simplifies to H visc =ṁv 2 rel 2. Now that we have established why energy conservation dictates that H visc =ṁv 2 rel 2 must be true for TRMLs, it is useful to examine what parameter choices for our analytic model satisfy this constraint. Given choices of M rel , χ, and τ , which are the "physical" parameters that are set by the physics of the mixing layer of interest, it is crucial to realize that not all combinations of Pr and f ν produce solutions to TRML that satisfy H visc =ṁv 2 rel 2. In other words, Pr and f ν are "model-specific" parameters that are constrained by the energy conservation criterion, and the problem at hand is: how to pick the values for Pr and f ν ? In § 3.3, we will show that enforcing energy conservation reduces the two degrees of freedom in Pr and f ν by one. In § 3.6, we will compare our analytic model with 3D simulations produced by Fielding et al. (2020) to clear up the remaining degree of freedom.

Energy Sources & Sinks vs. Pr and f ν
In this section, we keep all other parameters fixed at fiducial values and focus on how Pr and f ν affect the energy conservation condition by plotting the energy sources and sinks in Eq. 17 as we vary Pr and f ν . vary with respect to Pr and f ν . Each vertical group of dots represents a solution to a mixing layer calculated using our model. In each panel, the parameter choice that satisfies the energy conservation condition H visc = mv 2 rel 2 can be found by identifying the intersection of H visc (the pink curve) andṁv 2 rel 2 (the black dashed curve). For example, in the left panel we fix f ν = 0.01 and vary Pr . H visc andṁv 2 rel 2 intersects at around Pr = 0.07. This is the only value of Pr that yields a physically sensible solution for the given choice of f ν . A similar situation holds true in the right panels. Once we fix either Pr or f ν , there is only one value of the other parameter that is consistent with energy conservation. In conclusion, demanding energy conservation (H visc =ṁv 2 rel 2) reduces one degree of freedom in the parameter space of our analytic model, providing an important constraint on the choices of Pr and f ν . We stress that this constraint does not come from simulations and is a natural result of enforcing energy conservation. In fact, in § 3.6, we will discuss how to use simulation results to further constrain the parameter space such that we can obtain best fit values of Pr and f ν , the two parameters that are introduced by our analytic model. Our controlled experiments so far in which we varying Pr and f ν are mainly for providing constraints on the viable choices of these parameters. Alternatively, it is also instructive to vary the other parameters of our analytic model that are set by the physics of the mixing layer of interest. Physically, this corresponds to comparing the energy budget across an ensemble of mixing layers and can potentially help us gain insights into what is controlling the phase structure. In this section, we explore the consequences of vary M rel on the energy sources and sinks in a mixing layer.
As shown in Fig. 10, with τ = 10 −1 , χ = 10 2 , Pr = 0.07, and f ν = 10 −2 , energy conservation is satisfied (i.e.,ṁ ∆c 2 s (γ − 1) + v 2 rel 2 = Q cool + H visc + W) for a range of M rel . The crucial point to note is that the choice of M rel affects the energy conservation criterion very weakly. Graphically, this is shown by the overlap of the curves representing H visc andṁv 2 rel 2 in the entire range of M rel plotted in Fig. 10. This characteristic ensures that the ensemble of mixing layers at different M rel (with all other parameters held fixed) all satisfy energy conservation, which makes comparing between them physically meaningful. Figure 9. The energy sources and sinks in Eq. 17 as a function of Pr (left) and fν (right). As in Fig. 6, blue is related to radiative cooling, pink is related to viscous heating, and green is related to (work done by) pressure. Additionally, we show the Bernoulli flux (as defined in § 3.2) in dark brown, the enthalpy flux in light brown, and the relative kinetic energy flux in dotted black. Each vertical group of dots represents a solution to a mixing layer calculated using our analytic model with the cosine vx profile. All mixing layers have M rel = 1.5, τ = 10 −1 , χ = 10 2 , Pr = 0.07, fν = 10 −2 unless either Pr or fν is explicitly varied. In each panel, there is only one parameter choice, denoted by the black arrow, that satisfies the energy conservation condition Hvisc =ṁv 2 rel 2 as discussed in § 3.2. Figure 10. The setup is similar to Fig. 9, but here we vary M rel for a series of mixing layers with the cosine vx profile, τ = 10 −1 , χ = 10 2 , Pr = 0.07, and fν = 10 −2 . Increasing M rel implies increasing Hvisc, as shown by the upward sloping Hvisc curve. As we saw in Fig. 6, most of the contribution to Hvisc comes from the region near the cold phase where viscous heating balances radiative cooling. Thus, as Hvisc (pink) continues to increase with M rel and becomes comparable to the enthalpy flux, we see an upturn in Q cool (blue) as a response.
As discussed in § 3.1, higher M rel implies more viscous heating (H visc ). This can be seen in Fig. 10 through the positively sloped H visc line (in pink). While analyzing the detailed structure of these mixing layers with the cosine v x profile prescription in § 3.1.1, we observed that most of the viscous heating is concentrated in the third ("viscous") zone within the cold phase, where viscous heating is balanced by radiative cooling. Thus, as H visc increases with M rel , Q cool increases together with it. This effect is most pronounced at the high M rel end of Fig. 10, where the curves for the enthalpy flux and the Bernoulli flux (see § 3.2) diverge and Q cool follows the Bernoulli flux. The importance of including the relative kinetic energy flux in our energy conservation argument is most apparent in the high M rel regime in which H visc can exceed the enthalpy flux (ṁ∆c 2 s / (γ − 1)). The essential point is that at low M the enthalpy flux far exceed the viscous heating, but the enthalpy flux scales with M 3/4 rel , as is found in 3D simulations (Fielding et al. 2020;, whereas the viscous heating scales with M 2 rel , so at sufficiently large M rel the viscous heating will always dominate. When the viscous heating begins to dominate the Q cool , which had been primarily following the enthalpy flux scaling, now must steepen to follow the viscous heating scaling. This upturn in Q cool is an important feature of the 3D simulations of the mixing layers (this point will be further explored in § 3.7). We have seen that this feature, as well as the cold end spike in the cooling flux distribution discussed in § 3.1.1, are both produced by the extended viscous zone at the cold phase of the mixing layer. We believe the same basic physics applies to the 3D simulations because they exhibit very similar features.
Although Fig. 10 does provide many useful insights, it does not reproduce all of the key features that are seen in the 3D simulations. Namely, in the 3D simulations we also see a downturn in the enthalpy flux with increasing M rel (we will later show this explicitly in Fig. 17), which is not a feature of our analytic model with the cosine v x profile setup. We therefore turn to the alternate version of our analytic model in which we prescribe the x-momentum ρv x as a cosine and infer the value of v x from that. The detailed setup of this model is introduced in § 2.7. Now we explore key differences of the mixing layer structure we obtained from the v x and ρv x profiles to see what the ρv x profile setup can tell us about the underlying physics of the 3D simulations.

Energy Sources & Sinks vs. M rel (with the Cosine ρv x Profile)
To better understand how the energy sources and sinks vary with M rel under the cosine ρv x profile prescription, it is instructive to first go back and examine the structure of a mixing layer under this new setup. Fig. 11 shows a mixing layer with the cosine xmomentum profile and fiducial choice of parameters. Compared to a mixing layer with the cosine v x profile (Fig. 6), the crucial difference here is that the T and v x profiles change almost simultaneously in position space. This is because we are tying dρ/ dz and dv x / dz together through the cosine ρv x profile, and since dρ/ dz and dT / dz are inherently tied through the ideal gas law, dv x / dz and dT / dz are also related, leading to the simultaneous change of T and v x . A direct consequence is that we lose the extended viscous zone at the cold phase and the corresponding spike in the cooling flux distribution. Instead, both the x-and z-component of viscous heating are now strongest at intermediate temperatures in the mixing layer.
Recall that our energy conservation argument introduced in § 3.2 states that energy radiated away through cooling is balanced by the enthalpy flux plus the viscous dissipation. As we turn up viscous dissipation by increasing M rel , as shown in Fig. 12, the amount of radiative cooling is unaffected because at intermediate temperatures, where viscous heating is the strongest, radiative cooling is significantly larger in magnitude. Thus, in order to maintain energy conservation, the enthalpy flux decreases with increasing M rel to compensate for the increase in viscous dissipation. Graphically, this is manifested by the downturn in the brown curve for enthalpy flux in Fig. 12.
We note that this downturn in enthalpy flux is another important feature of the 3D simulations. This investigation with the x-momentum profile allows us to link this feature with significant viscous heating at intermediate temperatures.

Finding Pr and f ν Values that Best Matchṁ and Q cool from 3D Simulations
We have argued that our analytic model with both the cosine v x and the cosine ρv x profiles are able to reproduce key features of the 3D simulations (upturn in Q cool discussed in § 3.4 and downturn in enthalpy flux discussed in § 3.5). To provide concrete evidence for this argument, we seek to directly compare our analytic model with the 3D simulations. To do so, we need to return to the question of finding appropriate values for Pr and f ν , the two additional parameters introduced by our analytic model. Remember in § 3.3 we have already discussed how we can reduce one degree of freedom in the parameter space of Pr and f ν by enforcing energy conservation. Here we constrain the remaining degree of freedom by requiring that our analytic solution with some M rel , χ, and tau reproduces key results of the 3D simulation with the same parameter choices. We are most interested in the mass, momentum, and energy transfer between the two phases and the shape of the phase distributions. Thus, it is the most useful to compare the values of the mass fluxṁ and the total amount of cooling in the mixing layer Q cool . In particular, we seek an analytic model that producesṁ and Q cool that are identical to its 3D simulation counterpart with the same set of physical parameters (M rel , χ, and τ ).
Suppose we were given a 3D simulation with some relative shear velocity M rel , characteristic cooling time τ , density contrast χ, and we select a time in the simulation such that we have a well-defined value for h, the thickness of the mixing layer. Then we can plug the values of these parameters into our analytic model to obtain the corresponding 1.5D solutions. Note that since there is still one degree of freedom in the Pr vs. f ν parameters, there is a suite of 1.5D solutions (with various combinations of Pr and f ν ) that satisfy energy conservation and can potentially match the 3D simulation we started with. Each of these 1.5D solutions, or each "energy-conserving" pair of Pr and f ν , yields a mixing layer with some values ofṁ and Q cool . Together this suite of mixing layers form a linear subset in theṁ-Q cool space. To determine which pair of Pr and f ν best matches the 3D simulation and how good the match is, we simply take the values ofṁ and Q cool from the 3D simulation as a point in theṁ-Q cool space and plot it together with linear subset in theṁ-Q cool space formed by the suite of mixing layers from our analytic model. The point on the linear subset that is the closest to the 3D simulation point corresponds to the best fit Pr and f ν values. Fig. 13 is a graphical representation of this idea. We have chosen the cosine v x profile setup and fiducial choices of parameters (M rel = 0.75, τ = 10 −1.5 , χ = 10 2 ) for this demonstration. For any given value of f ν , the corresponding value of Pr that satisfies the energy  Figure 11. A fiducial solution of our analytic model using the cosine ρvx profile. The parameter choices for this solution are: M rel = 0.75, χ = 10 2 , τ = 10 −1.5 , fν = 10 −3 , Pr = 0.3. The information in each panel and the color choices is the same as Fig. 6. However, a crucial difference is that mixing layers with the cosine ρvx profile has the temperature and velocity profiles change simultaneously, which means there is no significant viscous heating near the cold phase and thus no spike in the cooling flux distribution. Now viscous heating is mostly concentrated at intermediate temperatures and never dominant amongst the various contributions to the advected enthalpy.
conservation condition H visc =ṁv 2 rel 2 can be found by performing a bisection on Pr . We are able to perform this bisection because, as shown in the left panel of Fig. 9, when Pr is greater than the energy-conserving value, H visc >ṁv 2 rel 2, and when Pr is smaller than the energy-conserving value, H visc <ṁv 2 rel 2. (Notice on the right panel of Fig. 9 that the same holds true in the f ν parameter space, which means this bisection can also be perform on f ν for a given value of Pr .) Note that this bisection on Pr is different from the bisection on the mass fluxṁ as discussed in § 2.8. For each "guess" of Pr , we first bisect onṁ to find a temperature profile that reaches the cold phase smoothly. Then we compare H visc andṁv 2 rel 2 for the mixing layer with the eigenvalue ofṁ to inform our next guess of Pr . We repeat this bisection on Pr for many values of f ν to obtain a suite of energy-conserving pairs of Pr and f ν . Together these pairs form a curve on the bottom panel of Fig. 13.
Each energy-conserving pair of Pr and f ν corresponds to a mixing layer with someṁ and Q cool . Thus, for each point on the Pr -f ν curve on the bottom panel of Fig. 13, we can find corresponding values ofṁ and Q cool , allowing us to make plots ofṁ versus Q cool , f ν , and Pr , as shown in the top panels of Fig. 13.
We compare to a 3D simulation with the exact same parameter choices and include the result as a point in theṁ-Q cool space, with error bars denoting the 2σ temporal variations around when h in the 3D simulation is within 5 % of h = 1. The 3D simulation point appears to lie almost exactly on theṁ-Q cool curve from our analytic model (top left panel of Fig. 13), which means there exists some energy-conserving pair of Pr and f ν that, when plugged into our analytic model, very closely reproduces bothṁ and Q cool of the 3D simulation. The best-fit values of Pr and f ν (with error bars) can be found by projected theṁ value from the 3D simulation onto theṁ versus f ν and Pr spaces, as shown in the top middle and top right panels of Fig. 13. The inferred values of Pr and f ν are also shown on the Pr versus f ν curve in the bottom panel.  Figure 12. Similar to Fig. 10 but with the cosine ρvx profile, τ = 10 −1.5 , Pr = 0.3, χ = 10 2 , and fν = 10 −3 . As shown in Fig. 11, with the cosine ρvx profile we lose the region where viscous heating balances radiative cooling, which means increasing Hvisc with M rel does not affect Q cool . To maintain energy conservation (energy radiated away=enthalpy flux from hot phase+viscous heating), there is a downturn in enthalpy flux to account for the enhanced viscous heating.
In Fig. 14 and Fig. 15, we repeat the same exercise for several choices of M rel , χ, and τ , with both the v x and the x-momentum profile. We see that in all cases the 3D simulation points lie very close to the curves obtained from our analytic model, which means both the v x and ρv x profiles are able to reproduceṁ and Q cool values from the 3D simulations. The analytic model also has the additional benefit of helping us determine the turbulent Prandtl number, which is information not immediately available from the 3D simulations. The best fit Pr and f ν values for each choice of parameters can be determined using the approach illustrated in Fig. 13 and is shown in Fig. 16 as a function of M rel , χ, and τ . For the cosine v x profile setup, the best fit f ν values stay constant across different M rel , χ, and τ , and the best fit Pr values can be closely approximated by power-laws with respect to M rel and χ, and a logistic function in log space with respect to τ . As for the cosine ρv x profile, Pr stays almost constant (there is a very weak dependence on τ ) and f ν shows power-law variations with respect to all three parameters.
To summarize results from Fig. 16, we provide equations that allow us to calculate the best fit values of Pr and f ν as a function of M rel , χ, and τ . For the cosine v x profile, the best fit values of Pr and f ν are given by  Figure 13. Schematic of finding the best fit Pr and fν for a 3D simulation. Consider a 3D simulation with M rel = 0.75, τ = 10 −1.5 , and χ = 10 2 . We seek to find a solution using our analytic model with the cosine vx profile that best matches the 3D simulation. The complication is that our analytic model takes in two additional parameters, Pr and fν , as inputs. Imposing the energy conservation condition Hvisc =ṁv 2 rel 2 reduces one degree of freedom in the parameter space. The energy conserving pairs of Pr and fν are plotted in the bottom panel. Each pair of Pr and fν produces an analytic solution with its own mass flux (ṁ) and Q cool . The correspondingṁ vs. Q cool , fν , and Pr plots are shown in the top row. In the top left panel,ṁ is normalized byṁcrit = ρ hot c s,hot , and Q cool is normalized by Q0 = ρ hot c 3 s,hot L0. The grey point (with error bars) on the top left panel shows the 3D simulation values ofṁ and Q cool . Note that the 3D simulation point lies exactly on the curve generated from our analytic model, meaning that there is a combination of Pr and fν that, when plugged into our analytic model, is able to accurately reproduce bothṁ and Q cool from the 3D simulation. The green shaded region denotes the 3D simulation value ofṁ with 2σ uncertainty, which can be extended to theṁ vs. fν and Pr plots to find the corresponding best fit values of fν and Pr with 2σ uncertainty, which are the blue and orange shaded regions, respectively. The blue and orange shaded regions are also plotted together in the bottom panel.
For the cosine ρv x profile, the best fit values of Pr and f ν are given by  Figure 14. Theṁ vs. Q cool curves for solutions with the cosine vx profile at different choices of M rel , τ , and χ. As described in Fig. 13, theṁ vs. Q cool curves are obtained from the Pr vs. fν curve that satisfies energy conservation.ṁ is normalized bẏ mcrit = ρ hot c s,hot , and Q cool is normalized by Q0 = ρ hot c 3 s,hot L0. In each panel, theṁ vs. Q cool curve with the fiducial choice of parameters, M rel = 0.75, τ = 10 −1.5 , and χ = 10 2 , are colored in black. We vary χ (left), τ (middle), and M rel (right) and compare the resultingṁ vs. Q cool curves. Theṁ and Q cool combinations of the corresponding 3D simulations are also plotted. The 3D simulation points all lie on the curves from our analytic model, which means there is a combination of Pr and fν that, when plugged into our analytic model, is able to accurately reproduce bothṁ and Q cool from the 3D simulation. This also means that the 3D simulations satisfies the energy conservation condition Q cool = ∆c 2 s (γ − 1) + v 2 rel 2 ṁ, which is imposed to be true along the curves from our analytic model. we vary M rel , summarizes our analysis in § 3.4 and § 3.5. The cosine v x profile is characterized by the upturn in radiative cooling due to viscous heating near the cold phase, and the cosine ρv x profile is characterized by the downturn in enthalpy flux due to viscous heating at intermediate temperature, and the 3D simulations lie somewhere in between, exhibiting both features. More specifically, at low M rel both the enthalpy flux and radiative cooling scales as M 3 4 rel , while at high M rel radiative cooling turns up to scale as M 2 rel because the relative kinetic energy term is dominant. We see that, while results from both the v x and the ρv x profiles are close to the 3D simulations, neither is a perfect match. We experimented with a few other choices (e.g., a cosine shear kinetic energy profile) but did not find any to be significantly better; for the rest of this paper, we will explore both profiles with the expectation that they bound reality.
Although we do not yet have a perfect analytic model that is able to encapsulate both features and match the 3D simulation perfectly, this investigation still gave us many insights and allowed us to build physical intuition. Most notably, we conclude that there are two means by which viscous dissipation can impact the mixing layer: (i) by heating gas at intermediate temperatures, which reduces the enthalpy flux from the hot phase as the shear Mach number increases while leaving the scaling between total cooling and the shear Mach number unchanged (see § 3.5 and Fig. 12), and (ii) by heating gas in the cold phase, which causes radiative cooling to follow a steeper scaling relationship with the shear Mach number as the shear Mach number approaches unity while leaving the scaling of enthalpy flux with the shear Mach number unchanged (see § 3.4 and Fig. 10). These two modes of viscous heating play a crucial role in setting the phase structure of the mixing layer and the mass, momentum, and energy transfer between the two phases.
Besides understanding how the enthalpy flux and radiative cooling in the mixing layers scale as M rel , we also performed the same investigation while varying χ and τ . The results are summarized in the middle and right columns of Fig. 17. The key point here is that results from the two variations of the analytic model bracket the 3D simulation data point, which support our earlier claim that the analytic models bound reality. More specifically, both the enthalpy flux and radiative cooling show very weak dependence with respect to χ  Figure 17. Here we juxtapose the enthalpy flux and radiative cooling in mixing layers as a function of M rel (left column), χ (middle column), and τ (right column) from our analytic model with the cosine vx profile (dash-dotted lines), the cosine ρvx profile (solid lines), and from the 3D simulations (scattered data points). To find the 3D simulation data points and error bars, we extract the time-dependent 3D simulation results between t = 2tmix and t = 4tmix (where tmix = χ 0.5 L v turb , v turb is measured in 3D simulations), which is when the thickness of the 3D TRML matches the thickness defined in our analytic model. The data points are the median data between these times, and the error bars denote the 2σ temporal variations. For each point on the analytic model curves, we use values of Pr and fν that best fits the corresponding 3D simulation using the scaling relationships shown in Fig. 16. The fiducial choices of parameters are M rel = 0.75, χ = 10 2 , and τ = 10 −1.5 . Across all six panels, the curves calculated from the two variations of the analytic model lie near the 3D simulation results and appears to bound reality. In the left column, the cosine vx profile is characterized by the upturn in radiative cooling, and the cosine ρvx profile is characterized by the downturn in enthalpy flux. In the middle column, both the enthalpy flux and radiative cooling show weak dependence on χ. In the right column, both the enthalpy flux and radiative cooling scales as τ −1/4 , which is accurately reproduced by the cosine ρvx profile and the vx profile with constant Pr and fν (dash-dotted, semi-transparent lines). and scales as τ −1/4 . The τ −1/4 scaling is accurately reproduced by the cosine ρv x profile prescription but not by the cosine v x profile prescription. We recognize this as a limitation of our cosine v x profile. However, in the two panels in the right column, we have also included two semi-transparent curves that are calculated using the cosine v x profile at constant (not best fit) values of Pr and f ν (Pr = 0.07, f ν = 0.01). The semi-transparent curves follow the τ −1/4 scaling perfectly.
3.8. Matching the Shapes of the Phase Distributions with 3D Simulations Being able to match the values ofṁ and Q cool from 3D simulations tells us that we can successfully reproduce the mass transfer between the phases and the normalization of the cooling flux distribution. However, to be able to predict the detailed spectral features from the mixing layer, we usually need to rely on higher level features that depend on the shape of the cooling flux and column density distributions. In this section we explore the effectiveness of our analytic model in matching the shapes of the cooling flux and column density distributions from 3D simulations, which is a prerequisite for accurately predicting observational metrics.
In § 3.6, we have already described the procedure for choosing parameters in our analytic model to obtain the best fit to a 3D simulation. Now we simply compare the resulting cooling flux and column density distributions of our analytic model with that of the 3D simulation. Fig. 18 shows the cooling flux and column density distributions from a fiducial 3D simulation over-plotted  Figure 18. Comparison of the cooling flux and column density distributions from the 3D simulation and the corresponding best fit analytic solution with the cosine vx and ρvx profiles. The cooling flux distributions are normalized by F0 = ρ hot c 3 s,hot , and the column density distributions are normalized by N0 = ρ hot L0/ mH . The 3D simulation has M rel = 0.75, χ = 10 2 , τ = 10 −1.5 , and its cooling flux and column density distributions are plotted in black. The analytic models have the same M rel , χ, and τ , and the best fit Pr and fν are obtained using the procedure described by Fig. 13. For the cosine vx profile, the best fit Pr = 0.79 and best fit fν = 0.013. For the cosine ρvx profile, the best fit Pr = 0.40 and best fit fν = 0.00073. The phase distributions from the cosine vx profile (in blue) provides a better fit.
with the best fit analytic models using both the v x and the ρv x profile. The v x profile produces a close match except for a small mismatch in the slope of the cooling flux distribution at temperatures lower than the peak of the cooling curve and higher than the cold phase spike. The ρv x profile, on the other hand, produces a less accurate match. This is expected because the ρv x profile setup does not encapsulate the cold phase spike in the cooling flux distribution, which is an important feature that contributes significantly to the total amount of cooling Q cool . To compensate for missing the cold phase spike, the cooling flux distribution produced by the ρv x profile setup has a shallower slope at temperatures above the peak of the cooling curve, and as we saw in § 3.6, the resulting Q cool is remarkably accurate. Both profiles produce reasonable matches to the column density distribution, although the overall shape of the v x profile is a better match, for the reasons just mentioned. Much recent research has provided essential insights to relevant problems using both 3D hydrodynamical simulations and analytic models. We now put our theory in the context of these existing works. Gronke & Oh (2018) explained the importance of radiative cooling in TRMLs through cloud-crushing simulations (this had been seen phenomenologically before by Marinacci et al. 2010;Armillotta et al. 2016, and many times after, e.g., Li et al. 2020;Sparre et al. 2020;Kanjilal et al. 2021;Abruzzo et al. 2021). They presented a physical model to explain why sufficiently large clouds are able to grow in mass through radiative cooling of hot gas in the TRMLs. Similar results are also found in simulations of cold filaments embedded in a hot medium by Mandelker et al. (2020). Fielding et al. (2020) demonstrated that the cooling is isobaric in fully resolved 3D simulations, leading them to argue that the hot gas inflow is driven by turbulence mixing rather than pressure gradients. In solutions produced by our analytic models, there is generally little or no pressure gradient, consistent with the argument of Fielding et al. (2020).
Additionally, Fielding et al. (2020) found that the cold gas growth and entrainment is stronger under: (i) higher relative shear velocities (M rel ), (ii) higher density contrast (χ), and (iii) more rapid cooling (smaller τ ). In the language of our analytic model, stronger cold gas growth is equivalent to a larger value of the mass fluxṁ. We see thatṁ is positively correlated with M rel and χ, and negatively correlated with τ . These are consistent with the findings of Fielding et al. (2020).  and Tan & Oh (2021) explored a 1D model for TRML in parallel with their 3D simulation results. Their pioneering attempt of encapsulating the physics of complex TRMLs within a 1D model revealed a tremendous amount of insights. In particular, they experimented with a constant conductivity (with the constant value taken from Spitzer conductivity) and a temperature-dependent conductivity. Building from their work, the setup of our analytic model involves two main differences. First, we have included viscosity. Second, we use an effective turbulent conductivity and viscosity (as introduced in § 2.3). These different choices reflects different understandings of the fundamental physical process at play. While  and  believe that heat diffusion is dominated by thermal conduction, we believe that turbulence plays a more important role. Yang & Ji (2022) also identified the importance of kinetic energy in the energy balancing argument. Their result is similar to the Bernoulli flux argument we introduced in § 3.2. However, a crucial difference is that Yang & Ji (2022) focused on the slowly cooling regime while we are primarily interested in the rapidly cooling regime. Additionally, Yang & Ji (2022) found a separation of turbulent and mixing zones in their TRMLs, which they argue to be a feature that develops in the high Mach number limit. In this work, we generally focus on Mach numbers that are much lower than what Yang & Ji (2022) explored, but we do observe a qualitatively similar separation of the turbulent and mixing zone in our cosine v x profile setup. However, the underlying physics is probably fundamentally different in the fast cooling regime, and we emphasize that the turbulent zone (or what we called "the viscous heating zone" in § 3.1) that we observe is in the cold phase, while the turbulent zone of Yang & Ji (2022) is in the hot phase.
We emphasize that our 1.5D analytic model has proven to be successful in matching the values ofṁ and Q cool and the shape of the phase distributions from 3D simulations. This gives us confidence that the physical insights we gained from our analytic model also applies to the 3D simulations.

A User's Guide to Our analytic Model
On the most practical level, the utility of our analytic model is that it allows users to obtain the key features of any TRML with any combination of parameters in a matter of seconds, without having to run the much more computationally expensive 3D simulations. For the convenience of the reader, we have put together several annotated notebooks that allows the reader to easily implement our analytic model. Following this URL, the reader can find a series of four notebooks: two for calculating an analytic solution of a TRML for any given set of parameters (using both the cosine v x and ρv x profile), and two for finding energy conserving pairs of Pr and f ν using the method outlined in Fig. 13 (again, using both the cosine v x and ρv x profile).
Finally, we emphasize here that besides being able to reproduce key features of the 3D simulations at a fraction of the computational cost, our analytic model is also valuable as a flexible, stand-alone model. Under our analytic model, it is easy to test how changing different parameters of our model affects the physics of the resulting mixing layer. For example, in Appendix C, we explore how changing β cold and β hot of the cooling curve (see § 2.2 for more details) affect the phase distributions of the mixing layer. Similar exercises on changing other parameters of the model can be easily performed as well.

Links to Observations -Calculating Ion
Fractions and Column Densities Using Our Analytic Model Due to its diffuse nature, the circumgalactic medium (CGM) is often mapped by absorption-line spectroscopy, where diffuse gas is detected through its absorption of light from bright background objects (for a review see Tumlinson et al. 2017). Absorption line observations of the CGM have, for example, revealed that the CGM is a dominant reservoir of baryons on galactic scales and can help explain the deficiency of baryons within the virial radius as inferred from the stellar mass to dark matter mass ratio (Werk et al. 2014). Simultaneous measurements of multiple ions provide a valuable window into the turbulent and multiphase nature of the CGM across a broad range of redshifts (e.g. Werk et al. 2016;Rudie et al. 2019;Qu et al. 2022). Furthermore, measuring the column density of these tracer ions can reveal key characteristics of the galaxy formation process. For example, Tumlinson et al. (2011) found that there is a correlation between the column density of O VI in the CGM and the specific star formation rate of galaxies, indicating that the CGM plays an important role in the evolution of galaxies.
Absorption-line spectroscopy has been similarly powerful for understanding galactic winds (for a review see Veilleux et al. 2020). These observations have made it clear that galactic winds are a ubiquitous feature of star-forming galaxies (Weiner et al. 2009;Rubin et al. 2014), both in the local universe (e.g. Heckman et al. 1990;Martin 1999), and at high redshifts of z ∼ 2 (e.g., Rudie et al. 2012). These observations have been incredibly useful in determining how much mass is carried out of the galaxy by these winds and how this mass flux depends on galactic properties (e.g., McQuinn et al. 2019). Galactic winds are known to be highly multiphase (Strickland & Heckman 2009), which makes it essential to have a robust model for the phase structure in order to extract as much information as possible from these observations.
In this Section, we discuss how our analytic model provides a quick and accurate way of calculating the ion fraction and column densities in a given TRML appropriate for either the CGM or a galactic wind. These calculations allow us to make predictions on the observable absorption features and better understand the composition and characteristics of the host galaxy.
Our analytic model provides us with the temperature, pressure, and density profiles of the mixing layer, which we can then use to calculate quantities that are directly tied to observations, including the ion fractions and ion column densities in the mixing layer.
To facilitate these calculations, we make two changes to our model. First, instead of the piece-wise power-law cooling curve introduced in § 2.2, we adopt a realistic cooling curve. This reduces the flexibility of our model but allows for more accurate calculations. Second, we need to convert our results from a dimensionless form into physical units to ensure that our calculation is physically sensible. We have discussed how to do this conversion in § 2.1.
For any given analytic solution to a TRML, we can proceed to calculate the ion fraction as a function of position in the mixing layer. We use the ion abundances table from TRIDENT (Hummels et al. 2017), which is tabulated using CLOUDY (Ferland et al. 1998), under the assumption of ionization equilibrium assuming a background radiation field from Haardt & Madau (2012). We note that this background radiation field is more applicable to the (outer) CGM and might need to be modified if one considers scenarios that are close to the disk or near quasars. The ion abundances are functions of redshift, metallicity, temperature, and density. Using the temperature and density profiles of a TRML obtained from our 1.5D model, one can calculate the ion fraction at any position in the mixing layer after assuming some redshift and metallicity values based on the physics of the scenario at hand. For the purposes of demonstration, we assume a redshift of 0 and solar metallicity in the remainder of the section and compute the ion fractions and column densities in a fiducial analytic solution of a TRML. We emphasize that this serves as a demonstration of the model capabilities rather than a definitive prediction. In fact, the solar metallicity assumption is more accurate for galactic winds and might need to be modified for the CGM.
Under these assumptions, we can calculate the ion fraction at any position in the mixing layer. Furthermore, the number density of an ion can be calculated by where n is the total number density, A element is the abundance of the atomic species of interest (assuming solar metallicity), and f ion is the ion fraction. The column density of that ion can then be calculated by where the integral runs through the entire mixing layer. In Fig. 19, we plot the column density distribution generated by our analytic model, the ion fractions, and the column densities of a selected list of ions. We note that the column density of the ions that probe either the cold phase (e.g. Mg II) or the hot phase (e.g. O VII) is dependent on the amount of cold or hot gas that the TRML connects. Physically, this corresponds to probing the interior of the cold cloud or the extent of the volume-filling hot gas. Our analytic model is for the TRML at intermediate temperatures only and does not provide a robust prediction for the cold and hot gas content on either side of the TRML. This can also be seen from the fact that Eq. 27 for the column density distribution dN/ dlogT diverges at the cold and hot end of the solution. This is a limitation, but not an inaccuracy, of our model. To understand how this limitation affects our column density calculations, we calculate the column density of a given ion between T = (1 + δT )T cold and T = (1 − δT )T hot , where δT is a small, positive constant. In the bottom panel of Fig. 19, we plot the column densities calculated using δT = 0.01 and use the error bars to indicate column density values at δT = 0 (upper bound) and δT = 0.1 (lower bound). As shown by the error bars, the ions that are the most abundant at either the cold or the hot phase (e.g. Mg II, O VII) are mildly affected by the choice of δT , while most of the other ions are unaffected. We note that the process described above assumes ionization equilibrium, which may be significantly in error in some circumstances (Oppenheimer & Schaye 2013). An additional point of caution in interpreting Fig. 19 is an inconsistency in accounting for selfshielding in our calculations. Physically, as the ionizing photons make their way in from the hot phase towards the center of the cold cloud they will be absorbed. Eventually the photons will be fully absorbed and the assumption of photo-ionization equilibrium will no longer be valid. Where self-shielding becomes an important factor will depend on the details of the mixing layer profile, but will primarily impact only the coldest gas and therefore the ions with the lowest ionization potential. This effect is ignored in our tabulated cooling function, but is approximately taken into account in the Ploeckinger & Schaye (2020) tables we adopt. It will be straightforward in the future to self-consistently account for this by using CLOUDY to calculate the equilibrium ionization state for the specific density and temperature profiles from our TRML calculations.
We have prepared an annotated notebook that produces figures similar to Fig. 19 with any ion of interest. The notebook can be found at this URL. We emphasize that these column densities are for a specific choice of pressure, velocity, metallicity, and radiative background and so serves as a demonstration of the model capabilities rather than a definitive prediction.

Links to Observations -Calculating Surface
Brightness Using Our Analytic Model Although harder to come by than absorption line observations, direct observations of the emission from galactic winds and the CGM provides an extremely powerful probe. Absorption lines generally give us a single pencil beam view per system, but emission observations have the potential to provide valuable insight into the spatial structure of these systems. New instruments are making emission observations much more readily accessible. Specifically, KCWI and MUSE are providing a promising path forward to greatly improve our understanding galactic winds and the CGM as has been recently demonstrated (e.g., Hayes et al. 2016;Rupke et al. 2019;Burchett et al. 2021;Reichardt Chu et al. 2022). As these new observations become available our need for robust models to interpret them becomes paramount.
As discussed in § 4.3, our analytic solution provides us with the temperature, pressure, and density profiles of a mixing layer, which, after being converted to physical units (see § 2.1 for more details), allows us to determine the emissivity (e), and thus the surface brightness of any given emission line of interest. To find the emissivity of a given emission line, we use the properties of shielded gas tabulated by Ploeckinger & Schaye (2020), which uses CLOUDY (Ferland et al. 1998) to calculate the equilibrium cooling and heating rates in the presence of a UV background based on Faucher-Giguère (2020), interstellar radiation and cosmic rays that depend on local gas properties through the Kennicutt-Schmidt (KS) relation (Kennicutt 1998), and dust. The emissivity tables in Ploeckinger & Schaye (2020) are given as a function of redshift, temperature, metallicity, and density. We assume a redshift of 0 and solar metallicity (see § 4.3 for a discussion of these choices), which means the emissivities are functions of temperature and density only (e(T, n)). Given the temperature and density profiles from our analytic solution, we can determine the emissivity of any given emission line as a function of position in the mixing layer. The surface brightness of that emission line can be calculated by SB = e(T (z), n(z))dz, where SB is the surface brightness, e is the emissivity, and the integral runs through the entire mixing layer. In Fig. 20, we plot the cooling flux distribution generated by our analytic model, the emissivities, and the surface brightness of a selected list of emission lines. As in § 4.3, we calculate the surface brightness of a given ion between T = (1 + δT )T cold and T = (1 − δT )T hot , where δT is a small, positive constant. In the bottom panel of Fig. 20, we plot the column densities calculated using δT = 0.01 and use the error bars to indicate column density values at δT = 0 (upper bound) and δT = 0.1 (lower bound). As shown by the error bars, most of the emission lines that are being plotted are unaffected by the choice of δT .
We have also prepared an annotated notebook that produces similar figures with any emission of interest. The notebook can be found at this URL.

Caveats section
In this section, we discuss some physical effects that are not fully captured in our model and identify possible directions of future work.   Realistic Thermal Conduction and Viscosity: In our analytic model we rely solely on the effective turbulent thermal conduction and viscosity, however in reality there will also be real conduction and viscosity that arises from particle collisions (Spitzer 1962). The strength and scaling of the conduction and viscosity coefficients in the tenuous plasmas that we are considering is uncertain. There can be large suppression possible due to the anisotropic motion of charged particles along the magnetic fields and by the excitation of resonant whistler waves (Roberg-Clark et al. 2016;Komarov et al. 2018;Drake et al. 2021).
Nevertheless, the impact of these missing process can be assessed and could be added in the future in a straightforward manner. The standard Spitzer conductivity of, for example, a system with p hot /k B = 1000K cm −3 is κ Sp (10 6 K) = 10 24.5 cm −1 s −1 and κ Sp (10 4 K) = 10 19.5 cm −1 s −1 . We can directly compare this to the effective turbulent conductivity of this system if v rel = 100 km/s, and h = 100 pc, which yields a characteristic value of κ turb ≈ n mix hv rel ≈ 10 25.5 cm −1 s −1 , where n mix = (p hot /k B )/T mix = (p hot /k B )/ √ T hot T cold . Therefore, we expect that the effective turbulent conductivity will dominate even over full strength Spitzer conduction. This is also true for the viscosity since the Spitzer Prandtl number is Pr Sp = m e /m p ≈ 1/42, while the effective turbulent Prandtl number is Pr ≈ 1/10. Nevertheless, there may be appreciable changes to the detailed structure of a TRML when realistic conduction and viscosity are included, and it would be straightforward to add to our existing TRML model.
Fractal Nature of TRMLs: The fractal nature of the thin corrugated surface where all the cooling takes place in 3D simulations of TRMLs identified by Fielding et al. (2020) is not directly captured in our analytic model. Remarkably, missing this physics does not prevent the analytic model from reproducing the mass flux and phase distributions of the 3D simulations. This suggests that our effective turbulent conductivity and viscosity prescription and the introduction of the "model-specific" parameters Pr and f ν can be easily tuned to produce the same cooling volume.
Non-equilibrium Cooling: The cooling curve we are using assumes equilibrium cooling, which can be incorrect when cooling is rapid (e.g., Gnat 2017). However, our analytic model is flexible enough that it should be relatively straightforward to implement a setup with non-equilibrium cooling.
Magnetic Fields and Cosmic Rays: Our analytic model does not include the effects of magnetic fields and cosmic rays. These components can suppress Kelvin-Helmholtz instabilities (Ji et al. 2019) and possibly affect the phase structure and mass transfer between the phases.
In this work, we present an analytic, steady-state, 1.5D model for Turbulent Radiative Mixing Layers (TRMLs) that includes, for the first time, a physically motivated model of turbulent conduction and viscosity that is key for setting the flow of gas through the TRML. We first summarize our trajectory in this paper before reviewing the key results.
In § 2, we start from mass, momentum, and energy conservation to derive a system of second order, coupled differential equations that describe the TRML. We then define the dimensionless numbers, χ, M rel , and τ , that most strongly affect the structure of TRMLs and discuss conversion between code and physical units. Finally, we introduce other features of our model, including a piecewise power-law cooling curve and, most crucially, an effective turbulent conductivity and viscosity that are proportional to the shear velocity gradient ( dv x / dz). We explore two forms for the shear velocity, one motivated by the v x profile in simulations, and the other by the ρv x profile. The introduction of a physical model for the turbulent conduction and viscosity which depends on a prescribed velocity model is a key new element of the model and allows for a much more selfconsistent solution.
We solve this system of equations numerically to obtain solutions, and in § 3, we present those results. We first examine individual solutions to the TRML obtained from the v x profile and identify the viscous heating zone that takes up a tiny fraction of the temperature space near the cold phase but a significant fraction of the position space. Next, we analyze the physics of TRMLs through the integrated energy sources and sinks and argue that the criterion H visc =ṁv 2 rel 2 must be satisfied for energy conservation.
By looking at how the total cooling in the mixing layer (Q cool ) and the enthalpy flux from the hot phase vary as a function of M rel , we conclude that there are two components of viscous dissipation in the TRML that controls the phase structure and mass transfer between the phases: viscous dissipation at intermediate temperatures reduces the enthalpy flux from the hot phase, and viscous dissipation at the cold phase enhances radiative cooling. The energy conservation criterion H visc =ṁv 2 rel 2 reduces one degree of freedom in the parameter space of our analytic model. In other words, Pr and f ν 5 , the two parameters that are introduced by our analytic model and not present in the 3D simulation, can effectively be reduced to one. Comparing with 3D simulation results allows us to identify the pair of Pr and f ν that best matches theṁ and Q cool for each of the two forms of the shear velocity profile (cosine fits to v x and ρv x ). (see Fig. 14 and Fig. 15 for more details). The best fit analytic solutions also produces a remarkably good fit to the phase distributions, especially with the cosine v x profile setup. We show this result in Fig. 18. In general, we expect the best fit analytic solutions produced by the two forms of the shear velocity profile to bound reality.
In § 4, we put this work in the context of the existing literature. We also link our analytic model with observations and discuss the calculation of ion fractions and ion column densities in the TRML. For the convenience of the reader, we have uploaded annotated notebooks for implementing our analytic model and calculating column densities of ions.
Our key findings are as follows: • Our model provides essential insights into the physics that controls the evolution of TRMLs. Specifically, any TRML that satisfies energy conservation must obey H visc =ṁv 2 rel 2, where H visc is the amount of viscous heating in the TRML,ṁ is the mass flux between the hot and cold phases, and v rel is the relative shear velocity between the phases. This viscous dissipation of relative kinetic energy is particular important in balancing radiative cooling as the shear Mach number approaches unity.
• By introducing two "model-specific" parameters, Pr and f ν , and calibrating their values using results from 3D simulations, we are able to reproduce the mass fluxṁ, total cooling Q cool , and phase structure of the computationally expensive 3D simulations (Fielding et al. 2020).
• We outline a quick and accurate way of calculating the column density of ions and surface brightness of emission lines in a given TRML using our model. The flexibility of our model allows for such calculations to be carried out for a wide range of physical conditions (pressure, metallicity, shear velocity), and we provide a demonstration of such calculations using a fiducial set of parameters.

APPENDIX
A. DETAILS OF OUR COOLING CURVE As mentioned in § 2.2, we use a piece-wise power-law cooling curve that is intended to reproduce the overall shape of the realistic cooling curves. Additionally, we include a heating term such that the cooling curve has two equilibria at both the cold and the hot phase temperatures. This heating term improves numerical performance but has no appreciable impact on our results.
Our cooling curve is given byĖ where β = β cold if T ≤ T peak β hot if T > T peak , log 10 ( T cold /T peak) log 10 χ , and α heat = (β cold − β hot ) log 10 ( T cold / T peak ) log 10 χ − β hot Note that τ is the dimensionless number that controls the normalization of the cooling rate, which we introduced in § 2.1.1, and T peak is the temperature where the cooling rate peaks and where the cooling time, defined to be is the shortest. As shown in Eq. A1 and Eq. A2, the dimensionless number τ serves as a normalization to both the cooling rate and the cooling time. In 3D simulations, this normalization is set by the shear time t sh = L/ v rel , where L is the shear layer length, which is well-defined in 3D simulations. However, here in our 1.5D model we do not have a transverse lengthscale, which precludes us from implementing the same setup. Instead, we use the characteristic lengthscale introduced in Eq. 25 to define the dimensionless number τ , which serves as the normalization in our setup.

B. DETAILS OF THE NUMERICAL INTEGRATION
We seek to numerically integrate 2 coupled differential equations to solve for the z-velocity (v z ) and temperature profiles. The equations at hand are second order differential equations in v z and T (Eq. 12 and Eq. 13) that we obtained from the steady state fluid equations. Our general strategy is to use scipy's solve-ivp integrator to perform the numerical integration after defining appropriate initial conditions.
We are primarily interested in condensation solutions that flow from the hot to the cold phase. The natural choice of the direction of integration is to follow the direction of gas flow: from hot to cold. This means we can deal with positive values of v z and the mass fluxṁ and avoid potential numerical difficulties. Now we are left with defining the initial conditions, which consists of T and v z at the hot phase, and the gradients of T and v z at the hot phase. The hot phase temperature and pressure is set to 1 for convenience, and the hot phase density can be calculated through ρ hot = P hot / T hot . Then for a given choice of the mass fluxṁ, the hot phase v z can be calculated through v z,hot =ṁ/ ρ hot . As for the initial gradients, we note that by definition our cooling curve has an equilibrium at the hot phase, so to produce non-trivial solutions that actually depart from the hot phase, we need to give the flow an initial "nudge" to push it away from the equilibrium and evolve towards the cold phase. To do so, we define the initial dT / dz to be a small negative number. The choice of the initial gradients does not affect the solutions we obtain as long as the initial gradients are non-zero and reasonably small (approximately 10 −6 ).
With the initial conditions defined, we proceed to numerically integrate our system of differential equations for a given value of mass fluxṁ. We use scipy's solve-ivp integrator with an rtol of 3×10 −14 and atol of 10 −11 and terminate the integration if T drops below the cold phase T calculated through the density contrast χ or if T exceeds the hot phase T (more on why this can happen in the following paragraphs).
Note that in the above description, the mass fluxṁ seems to be a free parameter. Now we impose the additional constraint that the temperature profile needs to be flat ( dT / dz = 0) at both the hot and cold phases by definition. For any given choice of parameters, there's only one value of the mass fluxṁ that allows for that. We call that the eigenvalue ofṁ (ṁ eigen ).
To find the eigenvalue ofṁ, it is crucial to realize that whenṁ is greater thanṁ eigen , the T profile shoots down to the cold phase with a steep slope, meaning that dT / dz < 0 at the end of the solution; whenṁ is less thanṁ eigen , the T profile barely reaches the cold phase before bouncing back up, and since we terminate the solution when T exceeds the hot phase T , we end up with dT / dz > 0 at the end of the solution in this case. Given the different signs of dT / dz at the end of the solution in these two cases and the fact that we want to find a solution with dT / dz = 0 at the cold phase, we adopt a bisection method in the parameter space ofṁ to pinpointṁ eigen . (see Fig. 5 for more details.) Specifically, we start off from the maximum possible value ofṁ given byṁ crit = P hot / √ T hot = ρ hot √ T hot and record the sign of dT / dz at the end of the solution withṁ =ṁ crit . Then we decrease our guess ofṁ by 10% and record the sign of dT / dz at the end of the new solution. We continue this process until we notice a sign change. Theṅ m eigen must be sandwiched between the two most recent guesses forṁ. Finally, we use scipy's optimize.root-scalar with rtol=10 −14 to find theṁ value that corresponds to the solution with dT / dz = 0 at the cold phase. This is the eigenvalue ofṁ we are looking for.
C. PHASE DISTRIBUTION AS A FUNCTION OF β cold AND β hot Besides serving as a computational efficient counterpart to the 3D simulation, our analytic model also has value as a flexible stand-alone framework. In particular, we can easily test how changing different parameters of our model affects the structure of the phase distributions. In this section, we focus on understanding how β cold and β hot , parameters of the cooling curve (see § 2.2 for more details), affect the shape of the phase distributions. Fig. 21 is made using the cosine v x profile setup and shows how the phase distributions vary as a function of β cold and β hot . Note that the values of Pr and f ν are carefully selected such that all the mixing layers being plotted obey energy conservation.
As we vary β cold , the slopes of the phase distributions at T < T peak changes linearly with it, while the slopes at T > T peak remain unchanged. The opposite holds true when we vary β hot , although the slope of the column density distribution at T > T peak seems to change only very slightly. Fig. 22 shows the same experiment but with the cosine ρv x profile setup. The same conclusions hold. In the top half of the figure we show the cooling flux distributions, column density distributions, and the power-law slope of these phase distributions for a series of mixing layers with different β cold . We use the cosine vx profile setup here, and all other parameters are held constant: M rel = 0.75, τ = 10 −1 , χ = 10 2 , Pr = 0.07, fν = 10 −2 , and β hot = 1. The cooling flux distributions are normalized by F0 = ρ hot c 3 s,hot , and the column density distributions are normalized by N0 = ρ hot L0/ mH . Note that Pr and fν are selected such that energy conservation is satisfied. The color scheme denotes the different choices of β cold . Generally speaking, changing β cold affects the power-law slopes of the cooling flux and column density distribution at T < T peak linearly and has almost no effect on the power-law slopes at T > T peak . The bottom half of the figure is identical to the top half other than that we fix β cold = −2 and vary β hot . Varying β hot affects the power-law slopes of the cooling flux and column density distributions at T > T peak linearly and has almost no effect on the power-law slopes at T < T peak .