MNRAS 525, 3831-3848 (2023) than its globally averaged value. On a similar note, Voit ( 2021 ) proposes that locally t cool / t ff ࣠ 1 leads to condensation but on a global scale the threshold condition depends on the amplitude of entropy fluctuations.

Turbulence plays a critical role in the evolution of the ICM. It is driven on large scales ( ∼100-500 kpc) by galaxy motions during mergers and on smaller scales by AGN ( ∼10-100 kpc). It can transfer the heat from the gas heated by AGN jets to the ambient ICM through turbulent mixing (Banerjee & Sharma 2014 ) and viscous dissipation. Further, Voit ( 2018 ) shows that turbulence can drive buoyancy oscillations that lead to condensation when 10 ࣠ t cool / t ff ࣠ 20. Gaspari et al. ( 2018 ) argue that the turbulent mixing time t mix is a more important time-scale than t ff , and the regions with cold gas are traced better by t cool / t mix ࣠ 1. Mohapatra & Sharma ( 2019 ) show that the onset of multiphase condensation is delayed when one drives turbulence on smaller scales, since t mix is shorter for small-scale driving.

Oli v ares et al. ( 2019) and O'Sulli v an et al. ( 2021) find that t cool / t mix ≈ 1 in regions of clusters, where the cold-phase gas is observed. Ho we ver, it is dif ficult to disentangle the importance of the two ratios ( t cool / t ff ) and t cool / t mix ) from observations, since (1) t cool varies more strongly with radius compared to t ff and t mix in cluster centres, and (2) we do not have many direct observations of turbulent velocities of the hot phase, except by Hitomi for the Perseus cluster (Hitomi Collaboration 2016 ). Hence, we rely on indirect methods of constraining turbulence and t mix (see Simionescu et al. 2019 for a re vie w).

Turbulence plays a dual role in multiphase condensation. On one hand, turbulence drives large density fluctuations on the driving scale in the ICM, leading to multiphase gas condensation. On the other hand, turbulent mixing suppresses the density contrast and multiphase condensation. Baek et al. ( 2022 ) find molecular gas co-spatial with sloshing features seen in the X-ray emission, implying that the velocity field affects condensation locally. Using idealized simulations, Mohapatra, Federrath & Sharma ( 2020, 2021, 2022b ) have shown that the amplitude of turbulence-driven (other sources, e.g. cooling, buoyancy, jet/outflo ws can also dri ve density fluctuations) density fluctuations depends on the degree of stratification of the ICM, the turbulent Mach number and the nature of driving (solenoidal versus compressive modes). Howev er, man y previous theoretical and numerical studies of the ICM initialize seed density fluctuations by hand, independent of the gas turbulence.

In order to better constrain the conditions required for the onset of multiphase condensation and to separate the two proposed threshold ratios of the time-scales, we conduct high-resolution hydrodynamic simulations of turbulence in a stratified medium, including radiative cooling of the gas. In our study, density fluctuations develop naturally due to the large-scale turbulence driving. We vary four main parameters rele v ant to cluster haloes -(1) the strength of stratification, which controls t ff , (2) the strength and (3) the nature of turbulence forcing, which controls t mix and the amplitude of density fluctuations, and (4) the initial gas density, which controls t cool .

This paper is organized as follows. We introduce our model, numerical set-up and tools in Section 2 . Then we present our results and discuss them in the context of galaxy cluster haloes in Section 3 . We summarize our key findings regarding the two time-scale ratios in Section 4 . In Section 5 , we discuss some of the shortcomings of our model and set-up, missing physics and how they might affect our results as well as the future prospects of this work. Finally, we present our concluding remarks in Section 6 .

We use Euler equations to model the ICM, with acceleration due to gravity ( g ) and turbulence ( a ), radiative cooling with a rate density L , and thermal heating with a rate density Q as additional source terms. We assume an ideal gas equation of state with an adiabatic index γ = 5/3. We evolve the following equations:

)

where ρ is the gas mass density, v is the velocity, P = ρk B T /( μm p ) is the thermal pressure, μ is the mean particle weight, m p is the proton mass, k B is the Boltzmann constant, and T is the temperature. In the energy equation (equation 1c ), the total energy density is given by E and the cooling rate density L is given by

where n e and n i are the electron and ion number densities, respectively. We use the temperature-dependent cooling function ( T ).

The time-scales of interest in this study are -the gas cooling time t cool , the isobaric thermal instability growth time t ti , the sound crossing time t cs , the gas free-fall time t ff and the turbulent mixing time on the driving scale t mix . They are defined as follows:

)

where α characterizes the density dependence of the heating rate density Q , with Q ∝ ρ α . The sound speed c s is given by √ γ P /ρ. For a deri v ation of equation ( 3b ) using linear stability analysis, see section 4.1 in McCourt et al. ( 2012 ). The two scales L and ℓ driv denote the size of the system and the driving scale of turbulence, respectively. In our simulations, ℓ driv = L /2 and v ℓ driv ≈ v, so t mix L /(2 v).

We use a modified version of the FLASH code (Fryxell et al. 2000 ;Dubey et al. 2008 ), version 4, to solve equation ( 1a ) to equation ( 1d ) in our simulations. For time integration, we use the MUSCL-Hancock scheme (Van Leer 1984 ;Waagan 2009 ) with the HLL5R approximate Riemann scheme (Waagan, Federrath & Klingenberg 2011 ). We use a second-order reconstruction method that uses primitive variables and ensures that density and internal energy are positive. Our simulation domain size is the same as in Mohapatra et al. ( 2020 ) -we use a cuboidal box with L x = L y = L = 40 kpc and L z = 1.5 L = 60 kpc. The box is centred at the origin (0,0,0). We implement periodic boundary conditions along the x -and y-direction for all variables. In the zdirection, we implement diode boundary conditions for the velocity. For density and pressure, we fix the values in the guard cells to their initial values throughout the duration of the simulation. In addition to using a larger box along the z direction to minimize the effect of the boundaries, we further smoothly decay the source terms -turbulent acceleration a , gas cooling rate density L , and gas heating rate Q for | z| > L /2, where the weighting function w( z) is given by

We analyse the outputs from our simulations only in the central cubical region with

We set-up a gravitationally stratified atmosphere with a constant g oriented along the -ˆ z direction. Pressure and density follow exponential profiles along the z direction at time t = 0 and the gas is at hydrostatic equilibrium, given by

H is the scale height of pressure/density and P 0 , ρ 0 ( = P 0 / gH ) are the initial values of pressure and density at z = 0, respectively. The pseudo-entropy S = P / ρ γ has a scale height H S ( ≡ 1/[dln S /d z]) = H /( γ -1). Since γ = 5/3, H S > 0, and the equilibrium is conv ectiv ely stable. The degree of stratification is denoted by the Froude number Fr on the integral scale ℓ int and is given by

, where (5c)

and N = √ g/ ( γ H S ) is the Brunt-V äis äl ä oscillation frequency, and v is the rms velocity. The quantity E ( k ) denotes the velocity power spectrum.

To force turbulence, we use a spectral forcing method using the stochastic Ornstein-Uhlenbeck (OU) process to model a (Eswaran & Pope 1988 ;Schmidt, Hillebrandt & Niemeyer 2006 ;Federrath et al. 2010 ). 1 The autocorrelation time of the driving is set to roughly match an eddy turno v er time on the driving scale. We drive turbulence only on large scales, corresponding to 1 ≤ k | L /2 π ≤ 3, where k is the magnitude of the wave vector k . The power is a parabolic function of k , peaking at 4 π / L , which corresponds to ℓ driv = L /2. We consider two types of forcing in this study: (1) natural mixture and (2) compressive modes only. For a more detailed description of the turbulence driving, we refer the reader to section 2.2.1 of Mohapatra, Federrath & Sharma ( 2022b ).

We use the temperature-dependent cooling function from Sutherland & Dopita ( 1993 ) corresponding to Z /3 (a third solar) metallicity. To control the code evolution time-step set by t cool , we introduce cut-offs on the cooling rate based on the gas pressure ( P cut-off ) and temperature ( T Cut-off ). We switch off the gas cooling when the gas pressure or temperature drop below these cut-off values. We also set a ceiling on the gas density ( ρ ceiling ) abo v e which we switch off the cooling. The complete cooling function is given by

where H is the Heaviside function. We have set T cut-off = 10 4 K, which is also the lower limit of the cooling function in Sutherland & Dopita ( 1993 ). We fix P cut-off = P 0 /1000 and the ρ ceiling = 500 × ρ 0 . For faster time-steps, we modify the criterion for setting the global time-step of the code d t code , such that d t code = min(0.5 × sub factor × t cool, min , d t CFL ), where t cool, min is the minimum cooling time o v er the domain, d t CFL is the code time-step set by the Courant-Friedrichs-Lewy criterion and sub factor is the subcycling factor which we set to 25. We refer the reader to appendix C of Mohapatra et al. ( 2022b ) for a discussion of this implementation. Note that we resolve cooling at most times when we update the internal energy using subcycling.

To prevent a runaway cooling flow in the simulation, we implement a shell-by shell balance (in constant z shells) between the net energy lost due to cooling and the net energy added by turbulence and thermal energy input. We inject thermal energy into each shell at a rate Q ( z) proportional to the local gas density in each shell ( Q ∝ ρ in equation 1c and α = 1 in equation 3b ). Ho we ver, if the turbulent energy input exceeds the total energy lost in a shell due to cooling, we set Q ( z) = 0 and do not apply any additional cooling. We implement this energy balance at each time-step. Mathematically, the heating rate is given by

We define the turbulent heating fraction f turb as

where we carry out the volume integration over the region defined by

We set our initial conditions to model the dense central regions of CC clusters. We initialize the gas with a constant initial temperature throughout the domain, set to T 0 = 1.07 × 10 7 K, such that the initial sound speed c s0 = 500 km s -1 . We set the gas number density n 0 = 0 . 1 cm -3 , so ρ( t = 0) = n 0 μm p exp (z/ H ) (except for four low-density simulations, where n 0 is 2 times smaller). We drive turbulence on 20 kpc scales, which roughly mimics the size of X-ray cavities seen in the ICM (see e. Note. Column 1 shows the simulation label. The number following H denotes the scale height of the initial pressure/density profile in code-units. We show the type of turbulence driving in column 2. In column 3, we show the average Froude number Fr of the simulations. The fourth column shows the time at which multiphase gas condenses out of the hot phase through thermal instability for a simulation. We denote it as 'NA' if there is no multiphase gas condensation in the particular simulation. In columns 5 and 6, we show the volume-weighted rms Mach number and its compressive component M comp , respectively. In column 7, we show the volume-weighted standard deviations of velocity v. We show the average value of the ratio between the thermal instability time-scale t ti and important dynamical time-scales -the free-fall time-scale t ff and the turbulent mixing time-scale t mix in columns 8 and 9, respectively . Finally , in column 10, we show σ 2 s, hot , the square of the standard deviations of the logarithms of density of the hot phase. All time-averaged statistics in columns 3, 5, 6, 7, 8, 9, and 10 are averaged for t ≤ t mp for runs in which multiphase gas forms. Movies of simulations are available at this playlist.

steady state, the rms velocity of the gas is approximately 250 km s -1 for our fiducial runs, consistent with the observations by Hitomi in the core regions of the Perseus cluster (Hitomi Collaboration 2016 ).

The cooling function ( T ) ∝ T 1/2 for free-free cooling at T ∼ 10 7 K. Since Q ∝ ρ, this gives t ti ≈ (10 / 3) t cool , using γ = 5/3 in equation ( 3b ).

We have conducted a total of 16 simulations in this study, which are listed in Table 1 . By default, our simulations have 512 2 × 768 resolution elements, with 768 cells along the z-axis. Since L z = 1.5 L , the individual resolution elements (or cells) are all cubical, organized in a uniformly spaced Cartesian grid. Since we only use the central cubical region with | x | , | y | , | z| < L /2 for the post-processing of our results, the ef fecti ve resolution is 512 3 .

By default, we drive the natural mixture of turbulent modes (i.e. we do not remo v e either solenoidal or compressiv e components of a ; see Federrath et al. 2010 ). Our fiducial set consists of two simulations with different strengths of gravity/stratification (and different t ff ) labelled H 1.0 and H 4.0 (so the value of g is in the ratio 4:1). The number following H in the label denotes the scale height of pressure/density in the simulation in code units (i.e. with respect to L ). We repeat this fiducial set as we vary other simulation parameters in our set. To check the effect of the nature of turbulence forcing, we keep all other parameters fixed but set ∇ ×a = 0 (compressive forcing; see Federrath et al. 2010 ). These two runs are indicated by ζ 0.0 in the label, where ζ denotes the fraction of solenoidal modes. In order to vary t mix while keeping t ti and t ff constant, we have two sets of simulations with weak driving and strong driving, denoted as 'wdriv' and 'sdriv' in the labels, respectively . Similarly , to check the effect of a longer t ti , we repeat the fiducial set and compressive forcing set of simulations with half the initial density ( n 0 = 0.05 cm -3 ) and pressure, so that the initial temperature still stays the same. This doubles the initial t ti and t cool . These four runs are marked by 'ldens' (low density) in the label. To compare our results directly with previous studies without constant turbulent forcing, we switch off the turbulent forcing and repeat the fiducial set with seed density perturbations at t = 0. These are marked by 'NoTurb' in the run label. Finally, to check the convergence of our results, we have two higher resolution versions of our fiducial simulations with 1024 2 × 1536 resolution elements. These simulations are denoted by 'HR' in the label.

In this section, we present and discuss the results of our simulations. We have run all our simulations till t end = 2.344 Gyr. Thermal instability leads to cold gas condensing out of the hot phase in 8 out of our 14 simulations. For runs that form multiphase gas, we define the time at which cold ( T ࣠ 2 × 10 4 K) gas first forms (when the cold gas mass fraction m cold /m tot > 0 . 01 per cent ) as t mp and list it in column 4 of Table 1 . We have also listed some time and volumeaveraged statistics in Table 1 , such as Fr, the rms Mach number M , the rms velocity v, the average value of the ratio between important time-scales t ti / t ff and t ti / t mix , and the square of logarithmic-density ( s ) fluctuations σfoot_4 s, hot in columns 3, 5, 6, 7, 8, and 9, respectiv ely. F or runs that do not form multiphase gas, these quantities are averaged over the last 120 Myr of the simulation. For runs that form multiphase gas, these averages are calculated in the 120 Myr just before t mp , but after the first 100 Myr, so that there is some time for turbulence to grow. 2 We begin this section by briefly discussing some key statistical properties of the gas in the fiducial set and the compressive forcing set of runs. These are crucial to understanding the second part of our study, where we vary the simulation parameters such as the strength of the turbulence forcing and the cooling rate. In the later subsections, we mo v e our focus to the non-linear evolution of thermal instability in the system and how it is affected by the different parameter choices.

Three of these runs form cold gas through thermal instability, but the H 4.0 run does not. In Fig. 1 , we show the projections of gas density (volume-weighted, first column), temperature (massweighted, second column), and column density fluctuations (after dividing out the xy -averaged density profile) in the hot phase ( T ≥ 10 6 K, third column). These snapshots are plotted when the runs have the maximum mass fraction of cold gas ( m cold / m tot ) and at t = t end for the H 4.0 run. The insets in column 1 show the projections of gas density at t = 0. Clearly, the runs with H = 1.0 have stronger gradients in the initial density than the runs with H = 4.0.

Thermal instability produces large variations in density, with much stronger variations compared to the initial density gradient. In all runs that form multiphase gas, the dense regions correspond to cooler gas and the rarer regions correspond to hotter gas, as expected.

For the H 1.0 run, the cold clouds are misty, i.e. they are small in size and occur throughout the simulation domain. In comparison, the compressive driving runs show many large clouds, with size ∼ℓ driv = 20 kpc. These results are similar to what we observed for different forcing runs in simulations without gravity in Mohapatra et al. ( 2022b ).

For the H 4.0 run, the net variations in density and temperature are much smaller compared to the other runs. Column density fluctuations in the hot phase are also much weaker for this run. For the other runs, we find that the regions with cold gas (in column 2) are associated with strong, positive fluctuations in the column density in the hot phase (in column 3). Such features are also observed in multiwavelength observations of the ICM (see e.g. Werner et al. 2013 ;Anderson & Sunyaev 2018 ;Baek et al. 2022 ). In our simulations, the spatial o v erlap between the different phases could be either due to turbulent mixing with the cooler gas making the hot phase denser or the cold gas could have directly formed from these dense regions of the hot gas, which have shorter cooling time (since t cool ∝ ρ -1 ).

In Fig. 2 , we show the time evolution of the mass fraction of cold gas ( T ࣠ 2 × 10 4 K) in the first row, the v olume-a veraged M in the second row, f turb (defined in equation 6c ) in the third row and the standard deviation of logarithmic density of the hot-phase σ s , hot in the fourth row.

Cold gas forms at different times ( t mp ) for the three different runs. The time t mp is clearly affected by the driving, multiphase gas condensation occurs much earlier for the compressive forcing runs. This is due to the stronger seed density fluctuations generated by the compressive forcing, as seen in the fourth row of Fig. 2 . The ratio m cold / m tot initially increases, reaches a maximum value, and then decreases with time. The rate of decrease in m cold / m tot is much faster for the runs with stronger gravity (i.e. H = 1.0), since the cold clumps being heavier than the ambient hot gas, fall faster to the ne gativ e z boundary.

At initial times, M for all runs reaches values of 0.5-0.7. The turbulent heating fraction f turb is approximately a few × 10 per cent . Ho we ver, for the runs forming multiphase gas, we find that both M and f turb decrease at t = t mp . By design, the turbulent forcing amplitude remains the same throughout the duration of the simulation. Coldgas condensation is associated with the production of fast-cooling dense gas at intermediate temperatures (2 × 10 4 K ࣠ T ࣠ 10 6 K), which increases the cooling rate. This is compensated by an increase in the heating rate since we impose energy balance in z-shells. The rarer hot-phase gas is heated more (because L ∝ ρ 2 , Q ∝ ρ), which increases c s and decreases M .

At late times, the simulation reaches a steady state at a lower M but higher f turb . The atmosphere is hotter and has a smaller net cooling rate, such that f turb increases. For the H 1.0 run, after the removal of extra mass, the turbulent heating alone is sufficient to balance the reduced steady-state cooling rate ( f turb = 1).

Among the two fiducial runs ( H 1.0 and H 4.0), the hot-gas density fluctuations are slightly larger for the H 1.0 run for t < t mp . This happens because the H 1.0 run is more strongly stratified (Fr listed in column 3 of Table 1 ) compared to the H 4.0 run. Mohapatra et al. ( 2020Mohapatra et al. ( , 2021 ) ) showed that for weak and moderate levels of stratification (Fr 1) the density fluctuations increase with increasing stratification (decreasing Fr) for fixed M and driving. These larger seeds lead to multiphase condensation developing in the H 1.0 run (and a slightly shorter cooling time, whose effect we discuss later), whereas they do not develop in the H 4.0 run.

The hot-gas density fluctuations show a sharp increase at t t mp for the H 1.0 run -bringing its value closer to the amplitudes for the compressive forcing runs. Clearly, the density fluctuations due to multiphase condensation are much larger than those due to stratified turbulence at t < t mp . Using unstratified multiphase turbulence simulations in Mohapatra et al. ( 2022b , fig. 6 and section 3.5), we showed that these larger fluctuations are due to the strong compressive velocities during cold-gas condensation and the baroclinicity of a multiphase turbulent system.

In Fig. 3 , we show the mass-weighted probability distribution functions (PDFs) of the Mach number (first column), temperature (second column), and gas density (third column) for our fiducial and compressi ve dri ving sets of runs. The PDFs for the three multiphase runs are averaged from 1.4 to 1.64 Gyr and for the single-phase H 4.0 run, they are averaged from 1.4 Gyr till t end . We show the 1σ spread in PDF values as shaded regions. The runs forming multiphase gas show two strong peaks in all three PDFs, whereas the H 4.0 run shows a single peak. The two peaks correspond to the hot and cold phases.

The hot phase is subsonic ( M hot < 1) for all four runs, as is expected from ICM observations (Hitomi Collaboration 2016 , see

Simionescu et al. 2019 for a re vie w). The high M peak corresponds Downloaded from https://academic.oup.com/mnras/article/525/3/3831/7251496 by guest on 14 September 2024 to the supersonic cold-phase gas, which has much smaller sound speed. Since we use the same forcing scheme to drive turbulence in all four runs, the shapes of the distributions of M are quite similar for M 1. The small offsets can be explained by differences in the temperature/sound speed among the different runs.

In the temperature PDFs, we observe a strong cold-phase peak at T cut-off = 10 4 K and the hot-phase peak at T ∼ 10 7 -10 8 K. The features in the PDF between these two peaks correspond to the shape of the cooling curve that we use. The temperature of the hot-phase peak is higher for the compressive forcing runs.

In the density PDFs, the low-density peak corresponds to the hot phase and the high-density peak to the cold phase. The hot-phase gas has much lower density for the compressive forcing runs, while the density of the cold-phase peak is similar. Thus, the ratio between the densities of the phases χ = ρ cold / ρ hot is much larger for compressive forcing. This is caused by strong converging and diverging motions on the driving scale (Schmidt et al. 2009 ;Federrath et al. 2010 ;Seta & Federrath 2022 ). For the H 4.0 run, the density PDF is lognormal with a power-law tail at low densities. The low-density tail is a known feature of the PDFs when the adiabatic index γ > 1, also reported in Passot & V ázquez-Semadeni ( 1998 ), Federrath & Banerjee ( 2015 ), and Mohapatra et al. ( 2020 ).

In Fig. 4 , we show the joint mass-weighted PDFs of the logarithms of temperature and density, temporally averaged over the same duration as the 1D PDFs in Fig. 3 . The different lines show the nature of fluctuations: adiabatic ( δT / T 0 ∝ ( γ -1) δρ/ ρ 0 ), isothermal ( δT = 0) at 10 5.5 K and T cut-off = 10 4 K, isobaric ( δT / T 0 = -δρ/ ρ 0 ), and isochoric ( δρ/ ρ 0 ). From a theoretical viewpoint, understanding the nature of fluctuations is important to calculate the growth rate of thermal instability through the different fluctuation modes (Das, Choudhury & Sharma 2021 ). They are also useful to compare with observ ations. For instance, Zhuravle v a et al. ( 2018) inferred the mode of perturbations from X-ray observations of the ICM.

In our single-phase H 4.0 run, the fluctuations are composed of isobaric and adiabatic components. This is in agreement with the stratified turbulence simulations (without radiative cooling) of Mohapatra et al. ( 2020 ), where we showed that unstratified turbulence produces adiabatic fluctuations, and the fraction of isobaric fluctuations increases with increasing strength of the stratification.

For the multiphase runs, we observe some clear trends in the PDFs -the hot phase (10 6 -10 8 K) is isobaric, the intermediate temperatures are isochoric, with a drop in temperature around 10 5.5 -10 6 K and the cold phase is approximately isothermal at T cut-off . We reported the same features in the temperature-density joint PDFs in Mohapatra et al. ( 2022b , fig. 5), so they are not strongly affected by the stratification.

The isochoric drop at T ∼ 10 5.5 -10 6 K is associated with the peak of ( T ), where t cool < t cs . The cooling time for the gas at intermediate temperatures is quite short and such gas may not be able to attain pressure equilibrium. Ho we ver, some of this pressure drop could be due to our lack of resolution of the cooling length ( ℓ cool = min ( c s t cool )). Recent high-resolution simulations of multiphase systems such as Fielding et al. ( 2020 ); Abruzzo, Fielding & Bryan ( 2022 ) argue that this could be due to lower spatial resolution in large-scale boxes, which do not resolve ℓ cool . While resolving ℓ cool is important to model the properties of the cold phase after it forms, it is not necessary to determine when or where it forms. In this study, we mainly focus on the latter part, so we do not expect our results to strongly depend on resolution. We have checked our results for convergence in the Appendix. The TNG50 simulations (Nelson et al. 2020 ;Ramesh, Nelson & Pillepich 2023 ), which track the cold gas better than our fixed-grid simulations, do not show this isochoric drop. Ho we ver, this could be partly due to the orders of magnitude variation in halo pressure in TNG50 haloes (therefore the sharp isochoric temperature drop is not as clear), whereas the vertical extent of our simulation box is much smaller to have a large pressure variation.

Theoretical studies such as Voit et al. ( 2017 ) report that the large-scale entropy gradient is important to thermal instability. They propose that haloes in thermal balance (applicable to our set-up) with a shallower entropy gradient are more susceptible to condensation. In Fig. 5 , we show the z-shell averaged entropy profiles ([ S / S 0 ]( z), where S 0 = P 0 /ρ γ 0 ) of the hot gas ( T 10 6 K) for our fiducial and compressive forcing sets of runs at t = 0 and t = t end . For the three runs that form multiphase gas, we also plot the entropy profile at the onset of multiphase condensation ( t mp , denoted in the titles of the respective columns).

For the H 1.0 run, the entropy gradient is steep at t = 0, but it flattens out around the onset of multiphase condensation ( t = t mp ). This is due to turbulent mixing, which mixes the low-and high-entropy regions together and makes the entropy gradient disappear. After cold gas condenses and mo v es out of the box through the bottom z boundary, at t = t end the entropy increases by almost an order of magnitude. We find that the gas has redeveloped a weak entropy gradient at this time.

The single-phase H 4.0 run starts out with a much weaker entropy gradient compared to the H 1.0 run. Despite starting out with a flatter entropy gradient, this run never forms multiphase gas. By t = t end , its entropy gradient also disappears and its entropy value is slightly larger than that for the H 1.0 run just before condensation.

The two compressive forcing runs form multiphase gas fairly quickly. Our snapshots just before thermal condensation show that the initial entropy profiles have large-scale v ariations e ven within the first ∼300 Myr of the simulations. By this time, the turbulence is still developing, such that a large-scale entropy gradient has not been lost to the mixing. By t = t end , the average entropy for both runs increases by an order of magnitude. Unlike the H 1.0 run, we still observe a strong entropy gradient for the ζ 0.0 H 1.0 run. The large-scale entropy profile shows a very disturbed state for the ζ 0.0 H 4.0 run due to strong large-scale perturbations induced by the compressive forcing, which are not mo v ed out of the box by the weaker gravity.

In summary, we find that a smaller initial entropy gradient (larger H ) does not necessarily imply better thermal stability of the halo. The entropy profile can be strongly modified by large-scale turbulence, which can remo v e the initial gradients, giv en enough time ( H 1.0 and H 4.0 runs). Further, the different amplitudes of density fluctuations also play a key role -larger fluctuations can seed multiphase condensation even when the entropy gradient is steep.

Following the discussion on the role played by the entropy profile, we now mo v e our attention to the z shell-averaged values of the three important time-scales of the system t ti , t mix , and t ff (defined in Section 2.2 ). The ratio between these time-scales is expected to play Figure 5. The vertical profiles of entropy of the hot phase ( T 10 6 K, averaged along the x-y plane) at t = 0 (red dotted line) and t = t end (blue dash-dotted line) for our fiducial and compressive driving sets of runs. We also show the entropy profile at t = t mp (green dashed line, when cold gas has just started forming) for runs that form multiphase gas.

The variation of important time-scales for the hot-phase gas ( T 10 6 K)-t ti , t ff and t mix averaged in shells parallel to the z-axis for our fiducial and compressi ve dri ving sets of runs at t = 0. For runs in which multiphase gas forms through thermal instability, we also show t ti and t mix at the onset of multiphase condensation (at t plot = t mp , denoted in the column titles). For the single phase runs, we show these time-scale profiles at t plot = t end .

a key role in the thermal stability of the system and has been studied in theoretical (e.g. McCourt et al. 2012 ;Sharma et al. 2012 ;Gaspari et al. 2018 ), numerical (e.g. Prasad et al. 2015 ;Beckmann et al. 2019 ;Butsky et al. 2020 ), and observational (e.g. Voit & Donahue 2015 ; Oli v ares et al. 2019 ) studies. In Fig. 6 , we show these quantities for the hot phase ( T ≥ 10 6 K) at t = 0 and at the onset of multiphase condensation ( t plot = t mp ). For the runs that do not form multiphase gas, we set t plot = t end = 2.344 Gyr.

We start with an isothermal profile, so at t = 0, t ti ∝ ρ -1 (see equation 3a ). It varies exponentially with z, with a scale height H . The free-fall time t ff is a constant throughout space and time, since we fix g to a constant value.

For the H 1.0 run, the z-gradient of t ti flattens and its value decreases slightly, following the same trend as the evolution of the entropy profile shown in Fig. 5 . Around the time when cold gas starts condensing out of the medium ( t = t mp ), t ti / t ff = 3.87 ± 0.05 and t ti / t mix = 5.92 ± 0.08. This medium satisfies the instability criterion ( t ti / t ff ࣠ 10) proposed by Sharma et al. ( 2012 ) and produces multiphase gas. Ho we ver, Gaspari et al. ( 2018 ) argue that when t ti / t mix > 1, turbulent mixing should be able to stop multiphase gas from de veloping. Ho we ver, this criterion does not correctly predict the outcome of the H 1.0 simulation. By t = t end , cold gas condenses out and falls through the bottom z-boundary. In the new steady state, the hotter and rarer atmosphere has t ti ∼ 10 Gyr, t ti / t ff ≈ 80 (see movie Downloaded from https://academic.oup.com/mnras/article/525/3/3831/7251496 by guest on 14 September 2024 of time-scale profiles evolution in supplementary material or at this link) and is stable against undergoing further thermal condensation.

For the single-phase H 4.0 run, the evolution of t ti is similar to that of the H 1.0 run, but its average value is slightly larger. The ratio t ti / t ff = 2.19 ± 0.02 and t ti / t mix = 6.92 ± 0.04. For this run, the criterion by Gaspari et al. ( 2018 ) correctly predicts that multiphase condensation does not occur in this system, while the Sharma et al. ( 2012 ) prediction does not hold true.

The amplitude of seed density fluctuations plays a key role in determining whether the systems undergo condensation. The H 4.0 run has weaker seed density perturbations compared to the H 1.0 run (see row 4 in Fig. 2 ) and a slightly larger t ti / t mix . The relatively faster mixing of the weaker seeds successfully prevents cold gas from condensing out. The two compressive forcing runs have much larger seed density perturbations. Despite having t ti / t mix = 6.2 ± 0.5 and 5.4 ± 0.4 at t = t mp for the ζ 0.0 H 1.0 and ζ 0.0 H 4.0 runs, respectively, they both form multiphase gas. At t = t end , the ζ 0.0 H 1.0 run has a similar value of t ti as the H 1.0 run, albeit with larger variations due to the compressive forcing. In comparison, the ζ 0.0 H 4.0 run reaches a larger t ti in steady state, but a similar t ti / t ff ≈ 100.

Considering the importance of the turbulence driving for the formation of multiphase gas seen in the previous subsections, here we analyse four more runs, where we vary the strength of the turbulence forcing. In steady state, v ∼ 20-40 km s -1 for the two 'wdriv' runs and ∼ 400 km s -1 for the two 'sdriv' runs. Similar to Fig. 2 , in Fig. 7 , we show the time evolution of the m cold / m tot , M , f turb and σ s,hot . We present the z shell-averaged profiles of important time-scales (for the hot phase) in Fig. 8 .

Out of the four runs, H 4.0wdriv and H 1.0sdriv form multiphase gas, whereas H 1.0wdriv and H 4.0sdriv do not. First, we focus our discussion here on the 'wdriv' set of runs. Due to the weak forcing, these two runs are the most comparable to thermal instability studies that do not explicitly drive turbulence (such as Sharma et al. 2012 ;Choudhury et al. 2019 ). 3The turbulent eddy turno v er time for these two runs is around 0.5-0.7 Gyr. Due to the weaker forcing, turbulence is strongly stratified, with Fr 1. In this regime, Mohapatra, Federrath & Sharma ( 2021 , fig. 5) showed that density fluctuations decrease with increasing stratification, due to strong buoyancy forces limiting motions in the z-direction. This is clearly observed in our simulations (fourth row of Fig. 7 ) as the density fluctuations are smaller for the H 1.0wdriv run compared to those for the H 4.0wdriv run (for t 0.8 Gyr). The weaker seed fluctuations are thus unable to induce multiphase condensation in the H 1.0wdri v run, e ven though t ti / t ff = 6.2 ± 0.1. In Fig. 8 , we find that the weak forcing is unable to significantly modify the initial profile of t ti by t = t end , unlike the fiducial set, which flattened the z-profiles of t ti (and entropy).

For the H 4.0wdriv run, t mix ∼ t ti around 1.316 Gyr, when the driven turbulence is expected to reach a steady state. Due to the weak turbulent mixing between the z-shells, most of the cold gas condensation occurs from the lower half of the box, which has a smaller initial t movies of simulation in supplementary material or at playlist link). Compared to the ζ 0.0 H 4.0 t ti 2-5 Gyr at t = t end , which is an order of magnitude smaller. Thus, for weaker driving, the system does not lose as much to condensation during the simulation period of 2.344 Gyr.

The trend in the two 'sdriv' runs are similar to what we observe for the fiducial set -out of the two, the more strongly stratified H 1.0sdriv run forms multiphase gas, while the weakly stratified H 4.0sdriv run does not. There are a few differences -the initial density fluctuations are larger for H 1.0sdriv so the multiphase gas forms much earlier compared to the H 1.0 run from the fiducial set even before the z-profile of t ti is flattened by turbulent mixing.

Before the onset of multiphase condensation, the amplitude of fluctuations in the H 1.0sdriv and H 4.0sdriv runs around t = 0.2 Gyr are similar (in agreement with expectations from Mohapatra et al. 2021 , for M ∼ The k ey difference between the tw o is the shorter average t ti in H 1.0sdriv. Although t / t mix = 9.3 ± 0.2, it is still unable to stop multiphase gas from developing. In H 4.0sdriv run, the turbulent heating due to the strong driving ( v = 410 ± 20 km s -1 ) is more than to offset the cooling ( f turb 1). The gas heats up with time, showing a gradual decrease in M and a larger value of t ti at t = t end .

8. Similar to Fig. 6 , but for our for our weak (wdriv) and strong (sdriv) driving sets of runs. For weak driving, the weaker stratification run forms multiphase gas, while for strong driving, the stronger stratification run shows multiphase gas.

For the runs described in this subsection, we lower ρ 0 and P 0 by half compared to the fiducial set (so initial T is fixed). This doubles t ti , while t ff and t mix are unaffected. We show the time evolution of rele v ant quantities in Fig. 9 and the z shell-averaged time-scale profiles in Fig. 10 . These are low-density (or longer t ti ) counterparts to Figs 2 and 6 for the fiducial set.

We find that only the two compressive forcing runs form multiphase gas, while the natural forcing runs do not. Since t cool and t ti are doubled, t mp ∼ 500 Myr is also doubled for these runs compared to ∼250-300 Myr for the fiducial compressive set with the same parameters. These two runs show a clear decrease in M around t mp associated with the hot phase becoming hotter. Since the cooling is weaker, f turb is larger, roughly by a factor of 2 for all the low-density runs compared to their fiducial counterparts. The fraction f turb ≈ 30 per cent for the natural forcing runs and 50-100 per cent for the compressive forcing runs for t < t mp . For t > t mp , f turb decreases, similar to what we observe for the fiducial set.

In Fig. 10 , we find that turbulent mixing flattens the z profiles of t ti for both the natural driving runs. The average t ti / t ff = 6.6 ± 0.2, t ti / t mix = 10.4 ± 0.1 for H 1.0ldens and t ti / t ff = 4.52 ± 0.06, t ti / t mix = 14.6 ± 0.2 for H 4.0ldens run. The larger value of these ratios compared to the fiducial set, ensures that multiphase condensation does not occur in either of these runs.

F or the compressiv e forcing runs, the av erage values of t ti / t ff = 13.5 ± 0.6, t ti / t mix = 16 ± 3 for ζ 0.0 H 1.0ldens and t ti / t ff = 6.4 ± 0.3, t ti / t mix = 16 ± 2 for ζ 0.0 H 4.0ldens. Both of these ratios are much larger than 1. Both Sharma et al. ( 2012 ) and Gaspari et al. ( 2018 ) models would predict the ζ 0.0 H 1.0ldens run to not produce multiphase gas, contrary to what we find. 4 Ho we ver, the large density fluctuations due to the compressive forcing grow before either mixing or buoyancy can prevent them from becoming multiphase. By t = t end , t ti ∼ 10-30 Gyr similar to that of their fiducial counterparts, despite their longer initial t ti . Thus, σ s , t ti / t ff and t ti / t mix determine the final value of t ti rather than the initial value of t ti .

Here, we summarize our results from all our simulations and discuss them in the broader context of the conditions that lead to multiphase condensation in the halo gas. In Fig. 11 , we show the time taken to form multiphase gas normalized by the thermal instability time-scale ( t mp / t ti ) (first row), minimum values of the ratios t ti / t ff (second row), t ti / t mix (third row), and t ti /min ( t ff , t mix ) (fourth row)foot_9 as a function of the standard deviation of logarithmic density (normalized) for all of our 16 simulations. For runs that form multiphase gas, we show these values just before t mp and plot them as filled data points. For the runs that do not form multiphase gas, we plot the ratios at t = t end using unfilled data points. The coloured dashed lines show the time evolution of these quantities as a function of σ s prior to multiphase condensation (or the end of the simulation).

Out of our 16 simulations, 9 form multiphase gas. For the seven simulations that remain single phase till t = t end , we plot t end /avg( t ti ) as a lower limit to t mp /avg( t ti ), in the first row of Fig. 11 . The singlephase simulations are generally concentrated to the upper left part of the figure, whereas the multiphase simulations are to the bottom right. This denotes that larger density fluctuations aid the formation of multiphase gas. Among the runs that form multiphase gas, we find that we can further divide them into three subgroups. The forcing in the four compressive driving runs and the strong driving H 1.0sdriv generates large density fluctuations ( σ s 0.3) and the gas forms localised high-density pockets with a short cooling time. The multiphase gas forms in t mp ࣠ 0.5 t ti for these simulations. The remaining four multiphase runs form cold gas at t mp t ti . We note that the runs with stronger turbulence ( H 1.0 and H 1.0HR) have stronger density fluctuations but form multiphase gas later compared to the runs with weak or no turbulent forcing ( H 4.0wdriv and H 4.0NoTurb). This highlights that turbulence driving generates stronger density fluctuations b ut turb ulence mixing slows the onset of multiphase condensation. On the other hand, in the absence of mixing the amplitude of density fluctuations keeps growing with time for the H 4.0wdriv and H 4.0NoTurb runs till t = t mp (see the fourth panel of fig. 7).

In this subsection, we first discuss how the predictions of thermal instability criteria proposed by Sharma et al. ( 2012 ) and Gaspari et al. ( 2018 ) hold for our set of simulations. We also attempt to construct a modified condensation curve based on these two criteria for our simulations, taking into account the local variation in t ti due to density fluctuations, as well as the lognormal shape of the density distribution (and consequently t cool , since t cool ∝ ρ -1 ) before multiphase condensation occurs (e.g. see the density PDF for the H 4.0 run in fig. 3 ). Since condensation is a local phenomenon, i.e. dense pockets of gas with a short ratio of the time-scales can condense out even when the atmosphere is globally stable (also seen in Choudhury et al. 2019 ), we consider the minimum value of these time-scales in our criterion. The densest regions would have gas density ρ max ∼ ρ exp ( c 1 σ s ), where c 1 is a positive constant. As t cool ∝ ρ -1 , min ( t cool ) ∼ t cool × exp (c 1 σ s ). Similar to Voit ( 2021 ), we use an exponential condensation curve that depends on σ s , and which takes into account these local variations in t ti (or t cool ) due to density fluctuations.

Figure 11. First row: Scatter plot of the time taken to form multiphase gas normalized by the z shell-averaged thermal instability time-scale ( t mp / t ti ) versus the standard deviation in the logarithm of gas density ( σ s ) for all our runs. The filled points show runs that form multiphase gas, while the unfilled points show runs that remain single phase till t = t end . For the latter set of runs, we show the lower limits to the ratio, denoted by the upward facing arrows in the symbols. Second row: The minimum value of the ratio of t ti to the z shell-averaged free-fall time-scale ( t ff ) t ti / t ff with the same x -axis. Third row: Similar to the upper panel, but we show min ( t ti / t mix ), the minimum value of the ratio between the z shell-averaged t ti and the turbulent mixing time-scale ( t mix ) instead along the y -axis. Fourth row: Here, we show min ( t ti /min ( t mix , t ff )), using the minimum of t ff and t mix in the denominator instead. The black line corresponds to the condensation curve described in equation ( 7a ). For the third and fourth rows, the black dashed line is given by equation ( 7b ). It clearly separates between the single phase and multiphase runs in the fourth row. The coloured dashed lines show the time evolution of these ratios as a function of σ s till t = min ( t mp , t end ).

The importance of t ti / t ff Sharma et al. ( 2012 ) propose the criterion t ti / t ff ࣠ 10 for the onset of multiphase condensation. This is satisfied in all our simulations, barring the ζ 0.0 H 1.0ldens run. Yet 8 out of the 15 simulations do not form multiphase gas, indicating that turbulent mixing has a significant effect on the conditions required for multiphase condensation (also discussed in Banerjee & Sharma 2014 ;Voit 2018 ). We find that the simulations that form multiphase gas are concentrated to the bottom right part of the figure, where either σ s is large or t ti / t ff is short. This is in agreement with the findings of Choudhury et al. ( 2019 ), who showed that the min ( t ti / t ff ) required for cold gas to condense out depends on the amplitude of density fluctuations. They also showed that the min ( t ti / t ff ) for which the gas becomes multiphase for a given σ s (or amplitude of density fluctuations) rises steeply once σ s 0.5. This effect is seen for our compressive driving run ζ 0.0 H 1.0ldens which has t ti / t ff > 10 but still undergoes multiphase condensation.

We attempt to construct a condensation curve like in Voit ( 2021 , see their section 4) with the functional form

to separate between the single phase and multiphase runs. We choose c 1 = 6 from an empirical fit to our data. Ho we v er, we hav e two outlier runs, H 1.0 and its high-resolution counterpart H 1.0HR that have t ti / t ff ∼ 2 but still do not form multiphase gas. Since this curve ignores the importance of turbulent mixing of fluctuations, it is unable to predict the occurrence of multiphase condensation correctly for runs with strong turbulent mixing.

The importance of t ti / t mix Now, we discuss the effects of the ratio t ti / t mix on the multiphase condensation. As discussed earlier, Gaspari et al. ( 2018 ) propose that gaseous haloes become multiphase if t ti / t mix ࣠ 1 and remain stable otherwise. This criterion does not correctly predict the outcomes of our simulations, since 7 out of the 15 haloes with t ti / t mix > 1 form multiphase gas. We think this discrepancy may partly arise because Gaspari et al. ( 2018 ) use δρ/ρ ∝ M (or σ s ∝ M ) to derive the amplitude of density fluctuations in their study (based on the results from cluster-scale simulations in Gaspari & Churazov 2013 ), which would make the density fluctuations directly related to t mix . This is not in agreement with our results. Recent studies have shown that σ s depends on M , the degree of stratification (denoted by Fr or H S ) (Mohapatra et al. 2020(Mohapatra et al. , 2021 ) ) and the Mach number of the compressive component of the velocities (Konstandin et al. 2012 ;Mohapatra et al. 2022b ), which correctly predict the amplitude of σ s in our simulations. Thus, understanding density fluctuations in cluster environments is key to predicting the thermal stability of the halo gas. Similar to Section 4.2.1 , we attempt to construct a condensation curve of the form min ( t ti / t mix ) = c 2 exp ( c 1 σ s ). We set c 1 = 6 and c 2 = 1.8 empirically. This curve correctly predicts the outcome of simulations with σ s 0.1. Ho we ver, this criterion ignores the importance of t ff . Thus, it fails to predict the outcome of the two runs with weak/no driving and strong gravity ( H 1.0wdriv and H 1.0NoTurb), where min ( t ti / t mix ) 1 but min ( t ti / t ff ) is much larger.

Instead of using the two ratios t ti / t ff and t ti / t mix separately, we construct a new ratio t ti /min ( t mix , t ff ) by taking the minimum of the two time-scales in the denominator. Our new condensation curve is given by min

where c 1 = 6 and c 2 = 1.8 are empirically determined from fitting our data. As discussed in earlier works and in previous sections of this study, multiphase condensation is inhibited when either of these time-scales are short enough. We plot the minimum value of this new ratio against σ s in the third row of Fig. 11 . This new condensation curve clearly separates all the simulations into subsets of single phase (unshaded region) and multiphase (grey-shaded region). In the limit of weakly forced turbulence with a long t mix , multiphase condensation is predicted well by the t ti / t ff ratio. Similarly in the limit of weak stratification, the ratio t ti / t mix predicts whether multiphase condensation occurs. Our new combined criterion co v ers both of these cases.

Although the behaviour of the condensation curve in our study is similar to that of Choudhury et al. ( 2019 ) ( t mix t ff in their study), we find that our curve flattens to a smaller threshold min ( t ti / t ff ) in the limit σ s → 0. We think this difference arises because they plot min ( t ti / t ff ) and density fluctuations δρ at t = 0 in their condensation curve, whereas we show these values just before multiphase condensation occurs. We expect δρ to grow (for e.g. see H 1.0wdriv run in the fourth panel of fig. 7) and min ( t ti / t ff ) to decrease by t = t mp , which would make the results consistent with each other.

Predictability of the outcome of a simulation : Here, we discuss whether one can predict the occurrence of multiphase condensation for a given set of simulation parameters -namely Fr, M , M comp , and the ratio of pressure and entropy scale heights R PS . The dashed lines in the second, third, and fourth rows of Fig. 11 show the co-evolution of the corresponding ratios and σ s . Except for the H 4.0sdriv run, these ratios do not show significant variation with time (after turbulence reaches a roughly steady state). Hence, if one can determine the value of σ s using the simulation parameters, then one can predict whether multiphase condensation occurs. We find two expressions for σ 2 s in the literature rele v ant to the turbulence parameters in our simulations:

from Mohapatra et al. ( 2021 ) for subsonic stratified turbulence (where R PS = H P / H S = 0.67 for our simulations) and

from Konstandin et al. ( 2012 ) for compressively forced subsonic turbulence. As we show in Fig. 12 , equation ( 7c ) agrees well with the the measured value of σ s in our natural driving simulations (left column), except the 'wdriv' runs. Similarly, equation ( 7d ) accurately predicts the scaling with M comp for our compressively driven turbulence simulations. The 'wdriv' (where turbulence may not have saturated yet) and 'NoTurb' runs (where we seed initial density fluctuations by hand) do not show good agreement with either scaling relation. Importance of f turb: Among the simulations that do not form multiphase gas, most reach a steady state where the thermal energy lost due to radiative cooling is replenished by turbulence dissipation and thermal heating. The steady-state value of σ s varies only by a fe w per cent. Ho we ver, as seen in the third ro w of Fig. 7 , f turb > 7d). The measured σ s shows a remarkable agreement with equation ( 7c ) predicted values for the natural driving runs, except weak turbulent forcing ('wdriv' runs, which may not have reached a turbulent steady state yet). On the other hand, the compressive forcing ( ζ 0.0) runs agree well with the equation ( 7d ). The runs without driven turbulence ('NoTurb' runs) do not agree well with either of the scaling relations.

1 for the H 4.0sdriv run. Thus, the heating rate due to turbulence exceeds the net cooling rate (thermal heating is switched off to prevent further overheating). Initially, the strong turbulence drives large density fluctuations and the pink dashed line initially crosses o v er to the multiphase side of the condensation curve (in the fourth row of Fig. 11 ). However, within a few t mix , the gas is o v erheated, which increases the temperature, decreases M and σ s , and raises the value of min t ti . When f turb > 1, even when the gas properties instantaneously satisfy the condensation criterion, the gas can be heated up on time-scales t < t ti , and multiphase condensation is prevented.

Here, we discuss some of the shortcomings of our study and possible ways to address them. We also outline some future prospects of this work.

In this set of simulations, all our standard set of runs use 512 3 × 768 resolution elements to resolve the domain of size 40 2 × 60 kpc 3 . So the minimum length that we can resolve is ∼80 pc. In order to capture the turbulent mixing layers between the hot-and coldphase gas, as well as to reproduce the pressure-temperature phase diagrams, one needs to resolve the cooling length ℓ cool , which is orders of magnitude below our resolution limit. In particular, the clear evidence for isochoric cooling in Fig. 4 is an indication that cold gas has collapsed to the grid scale. At that point, the gas cannot be compressed anymore because of insufficient resolution, pressure equilibrium cannot be maintained, and the gas cools isochorically.

Further, the small-scale turbulence is also not well-resolved in this study. Hence, we have not analysed the scale-by-scale kinematics of the hot and cold phases here and leave it to a follow-up study.

We conduct two high-resolution simulations -H 1.0HR and H 4.0HR with 1024 2 × 1536 resolution elements. We present these in the Appendix. The results of the higher resolution simulations are similar to those presented in the main text. Ho we ver, our resolution is still far from what is required to resolve the cooling length ℓ cool , so although the convergence in the Appendix is encouraging it is far from a guarantee that the results would be the same if our resolution were sufficient to resolve all the key length scales in the problem.

Throughout the duration of the simulation, we constantly force turbulence on large scales. Further, to prevent the model from undergoing a global runaway cooling flow, we have applied a shellby-shell energy balance at all times. Instead of such a fine-tuned balance at all times, clusters are rather expected to undergo cycles of heating and cooling, where a cooling episode triggers strong feedback, heats the gas, and prevents it from further cooling (as seen in simulations, such as Prasad et al. 2015 ;Beckmann et al. 2019 ). In a future study, we plan to explore the effect of episodic turbulence driving and decay, to mimic AGN on-off scenarios.

The density-dependent heating model that we use in our simulations (defined in Section 2.4.4 ) is quite idealized. We have ignored other possible heating sources such as cosmic rays (Butsky et al. 2020 ;Kempski & Quataert 2020 ;Su et al. 2020 ), thermal conduction (Br üggen & Scannapieco 2016 ;Jennings et al. 2023 ), mixing of hot bubbles with the surrounding ICM (Banerjee & Sharma 2014 ;Hillel & Soker 2017 ), etc. We have also ignored the effect of magnetic fields in this study. Ji, Oh & McCourt ( 2018 ) have shown that magnetic fields, independent of orientation can destabilize buoyant oscillations and modify both the amplitude and morphology of density fluctuations, which are critical to understanding the onset of multiphase condensation. Wang et al. ( 2021 ) and Mohapatra et al. ( 2022a ) show that magnetic fields can modify the kinematics of both the hot and cold phases. We plan to conduct follow-up studies exploring the effects of some of these physical elements.

We have modelled the ICM as a plane-parallel atmosphere with constant acceleration due to gravity. Ho we ver, cluster atmospheres are expected to be spherical/elliptical. Choudhury & Sharma ( 2016 ) showed that the amount of cold gas condensing depends on the variation of g (or t cool / t ff ) along the radial separation from the cluster centre. The energy and mass budgets are also expected to be different in a spherical atmosphere, since the denser central gas has a smaller mass fraction. The hot gas would be able to expand and cool more easily compared to the plane-parallel atmosphere. We plan to look into the effects of the cluster geometry in a future study.

In this work, we have explored the conditions that lead to cold gas condensation from the thermally unstable hot phase in the intracluster medium. We have conducted 16 idealized simulations of a local box of size (40 2 × 60) kpc 3 including radiative cooling, density-dependent thermal heating and turbulent driving (in 14 out of 16 simulations). The important time-scales that go v ern multiphase condensation in such a system are (1) thermal instability time t ti ( ∝ t cool , the cooling time); (2) gravitational free-fall time ( t ff ); and (3) turbulent mixing time ( t mix ). A short t ti makes condensation more likely, whereas shorter t ff and t mix are expected to prevent condensation. Since t cool ∝ ρ -1 (gas density), the amplitude of logarithmic density fluctuations σ s is also an important parameter to determine local variations in t ti . The ratios between the aforementioned timescales of the systemt ti / t ff and t ti / t mix are important to predict the occurrence of multi-phase condensation. Here, we summarize the main tak eaw ay points of this w ork, focusing on the importance of these ratios: (i) In the limit of weak stratification, the ratio t ti / t mix predicts the occurrence of multiphase condensation. We find that turbulent mixing suppresses multiphase gas condensation even for runs with min ( t ti / t ff ) 2 (see H 4.0 run in Figs 2 and 6 ). This result is further corroborated by our findings in our strong turbulent driving set of runs (labelled 'sdriv', see Figs 7 and 8 ).

(ii) In our weak turbulence driving simulations (labelled 'wdriv') and simulations without constantly driven turbulence (labelled 'NoTurb'), we find the occurrence of multiphase condensation is predicted well by the t ti / t ff ratio (see Figs 7 and 8 ). Strong stratification suppresses multiphase condensation even when min ( t ti / t mix ) 1 in our H 1.0wdriv and H 1.0NoTurb runs.

(iii) Large density fluctuations al w ays increase the likelihood of multiphase condensation. Cold gas forms in our simulations with min ( t ti / t mix ) 1 and min ( t ti / t ff ) 10, if the turbulence driving promotes strong density fluctuations, such as for compressive driving (see ζ 0.0 runs in Figs 2 , 6 , 9 , and 10 ). This happens due to the formation of dense pockets of cold gas with short t ti . The dependence of multiphase condensation on σ s is clearly seen in Fig. 11 .

(iv) Thus the two ratios min ( t ti / t ff ) and min ( t ti / t mix ) collectively predict whether multiphase condensation occurs. In the limit that one of these ratios is much larger than the other, the larger of the two determines whether multiphase gas forms. Taking into account our findings abo v e, we propose a new condensation criterion that considers the importance of both t ff and t mix as well as the variability in t ti due to large density fluctuations, which we parameterize using σ s . Our new multiphase condensation criterion is given by min ( t ti /min ( t mix , t ff )) = c 2 × exp ( c 1 σ s ) with c 1 = 6 and c 2 = 1.8, empirically determined and shown in the bottom panel of Fig. 11 . When the minimum value of the ratio t ti /min ( t mix , t ff ) falls below this threshold, multiphase condensation occurs in our simulations.

(v) Unlike previous studies, we find that the entropy scale height does not al w ays play a significant role in determining whether or not a system forms multiphase gas. Turbulent mixing flattens the entropy gradient on scales smaller than the driving scale in a few mixing time-scales. Ho we ver, in the limit of weak or no turbulence, simulations with a steeper entropy gradient are more stable against thermal condensation.

(vi) Our simulations that form multiphase gas reach a second steady state after most of the condensed cold gas rains down through the bottom z-boundary. In this state, we find the value of min ( t ti /min ( t mix , t ff )) to be independent of the initial value of min ( t ti /( t mix , t ff )) (before the condensation begins). Instead, its steady-state value increases with the amplitude of turbulent density fluctuations.

Movies of projected density and temperature as well as timeev olution of z-a veraged time-scale profiles of different simulations Downloaded from https://academic.oup.com/mnras/article/525/3/3831/7251496 by guest on 14 September 2024

MNRAS 525, 3831-3848 (2023) are available as online supplementary material, as well as at the following links:

(i) Playlist of fiducial runs and compressi ve dri ving ( ζ 0.0) sets of runs;

(ii) Playlist of weak ('wdriv') and strong ('sdriv') driving sets of runs;

(iii) Playlist of low-density ('ldens') sets of runs;

(iv) Playlist of high-resolution ('HR') runs.

(v) Playlist of runs without external driving 'NoTurb'.