# Non-linear saturation and energy transport in global simulations of magneto-thermal turbulence in the stratified intracluster medium

Jean M Kempf, François Rincon

► To cite this version:

Jean M Kempf, François Rincon. Non-linear saturation and energy transport in global simulations of magneto-thermal turbulence in the stratified intracluster medium. Astronomy & Astrophysics - A&A, 2025, 694, pp.A25. 10.1051/0004-6361/202451788 . hal-04820267

HAL Id: hal-04820267

https://hal.science/hal-04820267v1

Submitted on 9 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Non-linear saturation and energy transport in global simulations of magneto-thermal turbulence in the stratified intracluster medium

J. M. Kempf* and F. Rincon

Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, Toulouse, France

Received 4 August 2024 / Accepted 19 November 2024

## ABSTRACT

**Context.** The magneto-thermal instability (MTI) is one of many possible drivers of stratified turbulence in the intracluster medium (ICM) outskirts of galaxy clusters, where the background temperature gradient is most likely aligned with the gravity. This instability occurs because of the fast anisotropic conduction of heat along magnetic field lines. However, the extent to which it impacts the ICM dynamics, energetics, and overall equilibrium is still a matter of debate.

**Aims.** This work aims to understand MTI turbulence in an astrophysically stratified ICM atmosphere, its underlying saturation mechanism, and its ability to carry energy and to provide non-thermal pressure support.

**Methods.** We performed a series of 2D and 3D numerical simulations of the MTI in global spherical models of a stratified ICM thanks to the finite-volume Godunov-type code IDEFIX and using Braginskii magnetohydrodynamics. We used well-controlled volume-averaged, shell-averaged, and spectral diagnostics to study the saturation mechanism of the MTI and its radial transport energy budget.

**Results.** The MTI is found to saturate through a dominant balance between injection and dissipation of available potential energy, which amounts to marginalising the Braginskii heat flux but not the background temperature gradient itself. Accordingly, the strength and injection length of MTI-driven turbulence exhibit clear dependencies on the thermal diffusivity. With realistic Spitzer conductivity, the MTI drives cluster-size motions with Mach numbers up to $M \sim 0.3$, even in the presence of strong stable entropy stratification. We show that such mildly compressible flows can provide about ~15% of the non-thermal pressure support in the outermost ICM regions close to the cluster accretion shock and that the convective transport itself is much less efficient (a few percent only) than conduction at radially transporting energy. Finally, we show that the MTI saturation can be described by a diffusive mixing-length theory, shedding light on the diffusive buoyant nature, rather than the adiabatic convective nature, of the instability.

**Conclusions.** The MTI seems relevant to both the dynamics and energetics of the ICM through radially biased magnetic fields that enhance the background Braginskii heat flux. Further work including externally forced turbulence, for instance, mimicking accretion-induced turbulence, is needed to assess its overall relative importance in comparison to other drivers of ICM turbulence.

**Key words.** instabilities – magnetohydrodynamics (MHD) – turbulence – galaxies: clusters: intracluster medium

## 1. Introduction

The hot ($k_B T \sim 5$ keV) and diffuse ($\rho \sim 10^{-27}$ g/cm$^3$) intracluster medium (ICM) plasma is susceptible to the magneto-thermal instability (MTI; Balbus 2000, 2001; Parrish & Stone 2005) in the periphery of galaxy clusters. This instability originates from a combination of features conducive to ICM outskirts: the dilute nature of its constitutive plasma and the inward alignment of the background temperature gradient with gravity.

Despite the seemingly weak magnetic field filling galaxy clusters ($B \sim 1-10$ µG; Govoni & Feretti 2004; Ferrari et al. 2008), ICM plasma is strongly magnetised. The Coulomb mean-free path of charged particles is 14 orders of magnitude larger than their gyration radii around the local magnetic field. This peculiarity makes the macroscopic transport of heat and momentum anisotropic with respect to the magnetic field lines (Braginskii 1965), which in turn strongly alters the buoyant response to perturbations of a stratified plasma layer threaded by weak magnetic fields. Balbus (2000, 2001), Quataert (2008) showed that the buoyant stability of such an atmosphere is no longer controlled by the sign of the entropy gradient ($N^2 > 0$), with N the Brunt-Väisälä frequency, as given by Schwarzschild's

criterion; Schwarzschild 1906) but by that of the temperature gradient. In particular, despite the stable entropy stratification of galaxy clusters (Cavagnolo et al. 2009), the ICM is unstable either to the MTI when the background temperature gradient is negative with respect to the radial direction (Balbus 2000, 2001) or to its sister instability, the heat-flux-driven buoyancy instability (HBI; Quataert 2008), in the opposite case, both in the limit of fast parallel conduction and weak magnetic field. The MTI is therefore relevant to the periphery of galaxy clusters, while the HBI is more likely to develop in the cool-core of relaxed clusters, which exhibits an increasing temperature with radius (Cavagnolo et al. 2009). Negative radial temperature gradients in galaxy clusters are actually ubiquitous, as observed in X-rays for example (Pratt et al. 2007; Leccardi & Molendi 2008; Ghirardini et al. 2019), reflecting the hierarchical history of dynamical structure formation through gravitational accretion (Peterson & Fabian 2006). As shown by McCourt et al. (2013), temperature profiles observed in the Perseus cluster (Simionescu et al. 2011), or in a larger sample of intermediate redshift clusters (Leccardi & Molendi 2008) are a natural consequence of the cluster's dynamical accretion history and only modestly affected by the thermal conduction. The gravitational assembly process appears to be the main driver of the temperature distribution in galaxy clusters.

<footnote>* Corresponding author; jean.kempf@irap.omp.eu</footnote>
Additionally, ICM turbulent fluid motions induced by such diffusive magneto-buoyant instabilities may induce a hydrostatic mass bias. The square of the Mach number $M$, defined as the ratio between the fluid bulk velocity $v$ and the speed of sound, $c_s = (\gamma p/\rho)^{1/2}$, with $\gamma$ as the adiabatic index and $p$ as the pressure, provides an estimate of the non-thermal turbulent pressure support found in galaxy clusters because it is proportional to the ratio of the dynamical to thermal pressures. Parrish et al. (2012b) argued that the MTI can reach Mach numbers as high as $M \sim 0.3$ and similar subsequent non-thermal to thermal pressure ratios at saturation in the extended outskirts of galaxy clusters. Such mildly compressible flows, if present in the outer ICM, would to some extent bias the X-ray and Sunyaev-Zeldovich (SZ) mass determination of galaxy clusters because these methods rely on the assumption of the ICM being at hydrostatic equilibrium (HSE; Vazza et al. 2018; Eckert et al. 2019; Angelinelli et al. 2020). These possible mass biases are of significant cosmological importance since precise measurements of galaxy cluster mass distributions could help alleviate, or at least elucidate, the current tension around the Hubble $H_0$ and the $\sigma_8$ cosmological parameters (Allen et al. 2011; Pratt et al. 2019).

Anisotropic thermal conduction, which is essential to the MTI, was first implemented in a large-scale structure formation simulation by Ruszkowski & Oh (2010). Since then, only Ruszkowski et al. (2011) specifically looked for any compelling signatures of the MTI in a simulation of structure formation through gravitational collapse, but they could not find any evidence of MTI-induced turbulence. This could be due either to the MTI being significantly inhibited by (viscous and/or perpendicular thermal) numerical diffusion or to structure formation flows energetically outperforming the MTI. As a third option, Ruszkowski et al. (2011) suggested that any MTI flows, if present, would be hard to characterise given the very complex structure of the velocity field induced by gravitational collapse. A more detailed characterisation of MTI turbulence in high-resolution idealised global simulations, as first attempted by Parrish et al. (2008, 2012b), is needed to bridge the gap between such models and large-scale structure formation simulations including anisotropic thermal conduction (Ruszkowski et al. 2011), as well as to clarify the apparently conflicting conclusions reached by these different approaches.

The issue of the MTI saturation level is all the more puzzling since two apparently different descriptions of the MTI turbulence at saturation have emerged during the past decade in works by Parrish et al. (2012b) and McCourt et al. (2013) and Perrone & Latter (2022a,b). The first two works claim that MTI-driven turbulence can be described by a standard mixing-length theory (MLT) involving the pressure scale height $H_p = (\partial_r \log p)^{-1}$; while Perrone & Latter (2022a,b) invok a locally balanced mechanism between the energy harvested from the background temperature gradient and that dissipated by parallel thermal conduction. Kempf et al. (2023) suggest that which of these descriptions correctly depicts MTI-induced turbulence might depend on the typical injection length of the turbulence $\ell_i$ and more specifically on how it compares with the typical scale height $H_p$ of the cluster. A precise knowledge of the MTI turbulence strength and structure at saturation would help better assess its detectability in large-scale simulations of structure formation (Ruszkowski et al. 2011) or in the context of future X-ray observations of ICM outskirts dynamics (Kempf et al. 2023). With the current study, we therefore also aim to clarify this issue through global models of MTI in a stratified ICM atmosphere and using phenomenological arguments.

The outline of the paper is as follows. In Sect. 2, we first recap the basic physics driving the instability. We then review the main descriptions proposed for the non-linear saturation of the MTI. In Sect. 3, we present the numerical model, diagnostics, and set of global simulations we performed to study MTI-driven turbulence in the presence of strong stable entropy stratification. These numerical simulations are then analysed in detail in Sect. 4, where we bring to light the saturation mechanism of the MTI in our models. We also document the dependency of the MTI structure and strength on the level of thermal conductivity at saturation. In Sect. 5, we introduce a diffusive MLT of the MTI, and show how the two previously introduced descriptions of MTI-driven turbulence blend into a common framework. We then make use of this theory and of the numerical results to deduce the levels of convective flux and non-thermal turbulent pressure support that can be expected from MTI-driven turbulence in the outermost ICM regions. Finally, the main conclusions, caveats, and perspectives of our study are presented in Sect. 6.

## 2. Phenomenology of the non-linear MTI

Before reviewing previously suggested descriptions of the MTI dynamics at saturation, we first recall the schematic linear mechanics of the instability. We picture a dynamically negligible horizontal magnetic field that threads a plasma layer with temperature gradient and gravity both pointing downwards. The initial magnetic field lines are therefore isothermal. When a fluid element is subsonically displaced from its initial equilibrium position, it remains thermally connected to its initial neighbours thanks to the dragging of the frozen-in magnetic field lines, in the limit of highly electrically conducting fluids. Provided that thermal conduction parallel to the magnetic field is fast enough, the perturbation is then isothermal. A rising blob will be hotter (and less dense) than the neighbouring fluid and keep rising through buoyancy, triggering the MTI. In the ICM, the latter can fully develop, with a maximal growth rate $\omega_T$ of about 0.5 to 1 Gyr$^{-1}$, up to a state of saturated turbulence, where $\omega_T = (-g_0/H_T)^{1/2}$ is the MTI frequency, $g_0$ is the gravitational acceleration, and $H_T = (\partial_r \log T_0)^{-1}$ is the scale height of the background temperature $T_0$.

As discussed by Kempf et al. (2023), different physical mechanisms have been invoked to describe the turbulent saturation regime of the MTI. We chronologically review the two main schools of thought. On the one hand, Parrish et al. (2012b), McCourt et al. (2013) proposed to apply mixing-length arguments borrowed from the solar convection community (Vitense 1953; Böhm-Vitense 1958; Hurlburt et al. 1984; Chan & Sofia 1989; Canuto 1990; Cattaneo et al. 1991; Abbett et al. 1997; Porter 2000) to predict both the intensity of temperature fluctuations (and the subsequent level of convective heat flux; McCourt et al. 2013) and the kinetic energy (Parrish et al. 2012b) induced by MTI-driven turbulence at saturation. The standard MLT, as applied to thermal convection, states that the temperature fluctuations induced by convective motions scale as (Vitense 1953)

$$ \delta T \approx (\partial_r T_0 - \partial_r T|_{\text{ad}}) \ell_m, \quad (1) $$

where $\ell_m$ is the mixing length and $\partial_r T|_{\text{ad}}$ the adiabatic gradient, that is the temperature change of a fluid element due to adiabatic expansion or compression as it vertically moves in a stratified atmosphere.

In this theory, the mixing length $\ell_m$ is taken to be the typical vertical length wandered by a blob of fluid before it loses its
coherence and mixes with the background on account of diffusive effects. The incompleteness of the MLT lies in the choice of the free parameter $\ell_m$, which is typically taken to be a significant fraction $\alpha \sim 0.1$ of the pressure scale height $H_p$ in stellar evolution models. Great efforts have been undertaken in the solar convection community to better constrain the mixing length thanks to numerical simulations (Chan & Sofia 1989; Abbett et al. 1997; Porter 2000). In the MLT as applied to thermal convection, temperature perturbations scale as $\propto (\partial_r T_0 - \partial_r T_{\text{ad}})$ (rather than simply $\propto \partial_r T_0$) because fluid elements are brought out of their equilibrium position adiabatically (rather than isothermally). In the case of the MTI, however, initial perturbations are isothermal rather than adiabatic on account of the fast parallel conduction along horizontal magnetic field lines:

$$ \frac{\delta T}{T_0} \approx \frac{\ell_m}{H_T}. \qquad (2) $$

A similar relation can be obtained for the velocity fluctuations, assuming that all the work performed by the buoyancy is converted into kinetic energy (Vitense 1953; Abbett et al. 1997; Porter 2000):

$$ v \approx \ell_m \omega_T. \qquad (3) $$

Therefore, $v$ must vary with $H_T$ (the precise power-law depending on how $\ell_m$ scales itself with $H_T$). This residual dependency was seen in the global numerical simulations of ICM from Parrish et al. (2012b), which certainly feature cluster-size motions $\ell_m \sim H_p$. Absent a better physically motivated prescription for the mixing length $\ell_m$, these authors were not able to fully exploit Eqs. (2) and (3). This limitation will be further explored, and overcome, in Sect. 5.1.

On the other hand, Perrone & Latter (2022a,b) argued that the MTI saturates once the energy locally injected by the instability (which critically depends on the thermal diffusion) into temperature fluctuations at a length-scale $\ell_i \ll H_T$, is dissipated by the parallel thermal diffusion itself. Perrone & Latter (2022a) showed that this balance can be formally rewritten as

$$ \frac{\nabla \delta T}{T_0} \approx H_T^{-1} \qquad (4) $$

and used to further deduce scaling laws for the MTI injection length and velocity fluctuations:

$$ \ell_i \approx (\chi \omega_T)^{1/2} / N, \qquad (5) $$

$$ v \approx (\chi \omega_T^3)^{1/2} / N. \qquad (6) $$

At first glance, the latter equation appears to be very different from Eq. (3), especially so if the mixing length $\ell_m$ is taken to be some fraction of the pressure scale height $H_p$, as suggested by McCourt et al. (2013). We emphasise that the scenario described by Perrone & Latter (2022a) strongly relies on the energetics of the Boussinesq approximation, which nonetheless looks quite different from that of the compressible Braginskii magnetohydrodynamic (B-MHD) equations. Therefore, its leading energetic balance remains to be further theoretically understood and numerically tested in the context of MTI turbulence driven at global scales in a highly stratified atmosphere representative of any astrophysical ICM. The resolution of this twofold problem is one of the main objectives of the present study.

In both cases, however, Eqs. (2) and (4) encapsulate the idea that steady-state turbulence is reached when temperature fluctuations rearrange so as to fight back the background temperature

gradient, that is the source of free energy. Conversely, they are also the only dimensionally consistent way to translate this physical behaviour. The two differences between these theories are the privileged length scale and mode of heat transport invoked: the MLT from Parrish et al. (2012b), McCourt et al. (2013) is based on the presumed convective aspect of the MTI, and has been deduced from global simulations of MTI-driven turbulence, in which $\ell_m \sim H_p$; while the theory of Perrone & Latter (2022a,b) relies on the diffusive nature of the instability, and has been checked with Boussinesq simulations that inherently assume $\ell_i \ll H_T$.

Accordingly, the first simulations of MTI-driven turbulence performed by Parrish & Stone (2005, 2007), Parrish et al. (2008) showed that the initial background temperature can be brought back to isothermality on cosmological time scales, at least when it is not continuously supplied in energy by boundary heat fluxes. In these simulations, however, the respective contributions of convection and conduction in this isothermalisation process were not fully quantified. Under the homogeneous Boussinesq approximation used by Perrone & Latter (2022a,b), no similar conclusions could be reached because the evolution of the background temperature profile was not self-consistently modelled. Consequently, ambiguous conclusions have been drawn regarding the ability of the MTI to efficiently carry heat through convective motions in the ICM. We will aim at clarifying this point in this paper too. In the same manner, which of Eqs. (3)–(6) correctly describes the MTI velocity intensity at saturation may, critically or not, impact the level of non-thermal pressure support, and the subsequent estimates of hydrostatic mass bias, that can be expected from the turbulence induced by the instability. An objective of the current study is to provide further answers to this question. In the next section we describe the numerical tools used to achieve these goals.

### 3. Physical and numerical models

In this section we present the set of Braginskii magnetohydrodynamic (B-MHD) equations used in this study, and describe how they are normalised in the initial HSE model of a global stratified ICM common to all our simulations. Finally, we introduce the numerical methods used to solve these equations, and the energetics diagnostics used to analyse the simulations.

#### 3.1. The Braginskii magnetohydrodynamic equations

The compressible B-MHD equations are a system of non-linear partial differential equations describing the evolution of the density, the momentum and the energy of a perfectly ionised and dilute collisional plasma subject to anisotropic transport of heat and/or momentum (in a given spherical solution domain $\mathcal{V} = (e_r, e_\theta, e_\varphi)$ in the present study; see Sect. 3.3). In conservative form, the set of equations is

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (7) $$

$$ \frac{\partial(\rho\mathbf{v})}{\partial t} + \nabla \cdot \left(\rho\mathbf{v}\mathbf{v} + \mathcal{P}\mathbf{I} - \frac{\mathbf{B}\mathbf{B}}{\mu_0} - \mathbf{\Sigma}\right) = -\rho\nabla\Phi, \qquad (8) $$

$$ \frac{\partial E}{\partial t} + \nabla \cdot \left( (E + \mathcal{P}) \mathbf{v} - \frac{(\mathbf{B} \cdot \mathbf{v}) \mathbf{B}}{\mu_0} - \frac{\eta \mathbf{B} \times \mathbf{j}}{\mu_0} - \mathbf{\Sigma} \cdot \mathbf{v} + \mathbf{q}_B \right) = -\rho \mathbf{v} \cdot \nabla \Phi, \qquad (9) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{v} \times \mathbf{B}) + \eta \Delta \mathbf{B}. \qquad (10) $$
<table><caption>**Table 1.** Set of typical cluster scales used to non-dimensionalise the physical quantities of the simulations, computed for the Perseus cluster (Simionescu et al. 2011).</caption><thead><tr><th>Physical quantity</th><th>Normalisation</th><th>Reference value</th></tr></thead><tbody><tr><td>Length</td><td><i>R</i><sub>vir</sub></td><td>1.8 Mpc</td></tr><tr><td>Time</td><td><i>t</i><sub>dyn</sub></td><td>1.4 Gyr</td></tr><tr><td>Density</td><td><i>ρ</i><sub>vir</sub>/3</td><td>1.04 × 10<sup>-28</sup> g/cm<sup>3</sup></td></tr><tr><td>Velocity</td><td>√<i>G</i><i>M</i><sub>vir</sub>/<i>R</i><sub>vir</sub></td><td>1260 km/s</td></tr><tr><td>Thermal energy</td><td>μ<i>m</i><sub>H</sub><i>v</i><sub>vir</sub><sup>2</sup></td><td>10.3 keV</td></tr><tr><td>Conductivity</td><td><i>R</i><i>ρ</i><sub>vir</sub><i>v</i><sub>vir</sub><i>R</i><sub>vir</sub></td><td>9.70 × 10<sup>12</sup> g cm/s<sup>3</sup>/K</td></tr><tr><td>Magnetic field</td><td>√<i>ρ</i><sub>vir</sub>μ<sub>0</sub><i>v</i><sub>vir</sub></td><td>45.6 μG</td></tr></tbody></table>

Here, $\rho$ is the mass density, $v$ is the fluid velocity (and $v$ its norm), and $p$ is the thermal pressure (i.e. the trace of the usual stress tensor). The vector $B$ stands for the magnetic field (and the letter $B$ for its strength). The total pressure $\mathcal{P} = p + B^2/(2\mu_0)$ is the sum of the hydrostatic and magnetic pressures (and $\mu_0$ is the vacuum permeability). The gravitational acceleration field $g$, whose norm is $g_0$, equals the opposite gradient of the gravitational potential $\Phi$. The total energy volume density $E = E_K + E_M + E_I$ is the sum of the kinetic $E_K = \rho v^2/2$, magnetic $E_M = B^2/(2\mu_0)$, and internal $E_I = \mathcal{R}\rho T/(\gamma - 1)$ energy volume densities, where $T$ is the temperature of the gas and $\mathcal{R}$ is the specific gas constant (i.e. the ratio between the Boltzmann constant and the mean molecular weight). We close the system with the usual equation of state for a perfect gas $p = \mathcal{R}\rho T$. The magnetic diffusivity is represented by the letter $\eta$ and the current density by the vector $j = \nabla \times B/\mu_0$.

The heat flux $q_B$ in a dilute, collisional, plasma is given by Braginskii's closure (Braginskii 1965):

$$q_B = -\kappa (\hat{b} \cdot \nabla T) \hat{b}, \quad (11)$$

where $\kappa$ is the plasma thermal conductivity, and $\hat{b}$ is the unit vector in the direction of the magnetic field. In this study, we set the tensor $\Sigma$ to be the isotropic viscous stress tensor:

$$\Sigma = \mu (\nabla v + \nabla v^T) - \frac{2}{3} \mu (\nabla \cdot v) \mathbf{I}, \quad (12)$$

with $\mu$ as the plasma dynamic viscosity. A complete description of the transport processes occurring in such plasmas would require the use of Braginskii viscosity (Parrish et al. 2012a; Kunz et al. 2012) rather than the isotropic strain-rate tensor Eq. (12). Because the dynamics at the viscous scale are absent from both MTI theories reviewed in Sect. 2, the dynamics at larger scales, in which we are mostly interested in this paper, should, however, be unaffected by this choice. Accordingly, we neglect the effects of kinetic micro-instabilities possibly triggered by Braginskii viscosity (which is a disguised form of anisotropic pressure; Schekochihin et al. 2005, 2010), on the macroscopic heat transport (Riquelme et al. 2016; Komarov et al. 2018; Drake et al. 2021), as they are not well described within B-MHD. In the dilute ICM, the plasma transport coefficients $\kappa$ and $\mu$ are respectively set by the electrons and the protons (Kunz et al. 2012). The plasma thermal conductivity is given by (Spitzer 1962; Chandran & Maron 2004)

$$\kappa = 6.87 \times 10^{12} \left( \frac{k_B T}{5 \text{ keV}} \right)^{5/2} \text{ g cm s}^{-3} \text{ K}^{-1}. \quad (13)$$

In the ICM, the Prandtl number – defined as the ratio between the plasma kinematic viscosity, $\nu = \mu/(m_H n_i)$, and thermal diffusivity, $\chi = \kappa/(n_e k_B)$ – is – is Pr $\approx 0.01$ (Kunz et al. 2012), where $m_H$ is the proton mass and $n_e, n_i$ are the electron and ion number densities, respectively.

## 3.2. Normalisation, hydrostatic equilibrium, and parameters

In our simulations, lengths, times and densities are respectively normalised by the cluster virial radius $R_{vir}$, the dynamical freefall time $t_{dyn} = (GM_{vir}/R_{vir}^3)^{-0.5}$, and the third of the virial baryonic density $\rho_{vir}/3 = f_b M_{vir}/(4\pi R_{vir}^3)$, where $f_b = 0.17$ is the cosmic baryon fraction. To clarify expectations, the numerical values of these astrophysical scales are given in Table 1 for the Perseus cluster (Simionescu et al. 2011). With such a choice, physical quantities other than lengths, times and densities are self-consistently normalised (through the B-MHD equations) by other typical cluster scales, defined and gathered in Table 1. For instance, the velocity is non-dimensionalised to the virial velocity $v_{vir} = (GM_{vir}/R_{vir})^{0.5}$ and the temperature to twice the virial temperature $2k_B T_{vir} = \mu m_H v_{vir}^2$.

The 1D spherical hydrostatic ICM model of McCourt et al. (2013) is used as a common initial condition for all the simulations of the paper<fnref n="2" />. The corresponding density $\rho_0$ and temperature $T_0$ profiles are shown in blue in Fig. 1. This spherically symmetric model efficiently reproduces some key thermodynamical features of the ICM as natural consequences of the cluster accretion history, such as the non-zero temperature gradient despite the possible erasing effect of thermal conduction. More specifically, this model assumes fast relaxation of the ICM towards HSE between successive accretion events of baryonic matter shells (for a type III accretion history; McBride et al. 2009) in a pre-settled and independent Navarro-Frenk-White (NFW; Navarro et al. 1997) gravitational potential of dark matter (we use a concentration parameter of 5), without any effect of the thermal conduction on the entropy evolution. Although thermal conduction can be accounted for in this 1D spherical model, the resulting temperature profiles are never at thermal equilibrium (TE) in presence of isotropic conductivity (or, equivalently, Braginskii with purely radial field lines), that is they do not satisfy $\partial_r(r^2\kappa\partial_r T_0) = 0$ (for some given boundary conditions), because this condition is never imposed *per se* in this model. For further reference though, an ICM atmosphere at both hydrostatic and thermal equilibria with a Spitzer conductivity, Eq. (13), in a NFW potential well is shown in red in Fig. 1 too. This ICM atmosphere at both HSE and TE is qualitatively very similar to that of the McCourt et al. (2013) model when a conductive heat flux is taken into account (see the red and yellow curves in their Fig. 7 and 9). In our simulations (which initially have $\partial_r(r^2\kappa\partial_r T_0) \neq 0$), a modification of the background thermodynamic profile must be expected as soon as a non-negligible radial background heat flux operates. In Sect. 4.3 we will see that this adjustment manifests itself in simulations, after the opening of radial magnetic field lines by the onset of the MTI, through the development of a thermal wind, that remains even when a quasistationary state is reached, as previously noted by Zingale et al. (2002) in the case of an isotropic heat flux. In the linear phase though, no such dynamics are occurring, and the initial back-

<footnote n="1">¹ In the present paper, $\kappa$ (resp. $\chi$) is a thermal conductivity (resp. diffusivity) in units of g cm/s³/K (resp. cm²/s): their respective meanings are reversed with respect to Kunz et al. (2012), Berlok et al. (2021). Our definition of the thermal diffusivity also differs from that of the previous authors by a factor $n_e/(n_e + n_i)$, hence the different Pr.</footnote>

<footnote n="2">² We have successfully reproduced this model, which has already been used and partially described in Kempf et al. (2023).</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Stratified ICM atmosphere. The density $ρ_0$ (in g/cm$^3$, left); temperature $T_0$ (in keV, centre); and entropy $S_0$ (right) profiles are shown. In blue: Initial profiles common to all simulations. In yellow: Same profiles but averaged over time for the 3D run F0 with $κ_{eff} = 0.075$. In green: Same but for the 2D run S0 with a realistic Spitzer conduction (see discussion in Sect. 4.3). In red: ICM atmosphere at both hydrostatic and thermal equilibria in the presence of a hypothetical background heat flux with a Spitzer conductivity, Eq. (13), subject to the same boundary conditions on the temperature (for comparison only).</figcaption></figure>

ground profiles are at rest. In agreement with observations, the atmosphere predicted by this model is stably stratified (with respect to the Schwarzschild criterion, $N^2 = g_0/(γH_S) > 0$) in entropy $S_0 = \frac{R}{γ-1} \log(T_0ρ_0^{1-γ})$, with the ratio $N/ω_T$ of the Brunt-Väisälä to the MTI frequencies of order unity, though slightly radii-dependent (Kempf et al. 2023).

The MTI is excited through small initial white-noise velocity field perturbations, modelled by a Gaussian random field with an amplitude of $10^{-4}$ for each component. The initial magnetic intensity is weak with $B_0 = 10^{-4}$; and the corresponding volume-averaged plasma beta parameter, defined as $\langle β \rangle_V = 2\langle p_0/B_0^2 \rangle_V$, is ~$10^{10}$. The associated magnetic configuration is a divergence-free, radially modulated and purely azimuthal (i.e. horizontal) field, which is the magnetic geometry most unstable to the MTI:

$$B = B_0 \sin \left( \pi \frac{r-r_{\text{in}}}{r_{\text{out}}-r_{\text{in}}} \right) e_{\varphi}, \quad (14)$$

with $r_{in} = 0.3$ and $r_{out} = 1.5$ (so that the magnetic field is initially zero at the inner and outer boundaries). In doing so, the background heat flux is initially shut down and the initial ICM atmosphere is at HSE. Then, we expect the background thermodynamic profile to stay still until an appreciable background heat flux appears due to the opening of the magnetic field lines by the MTI.

The thermal diffusivity $\chi$ is a key quantity of MTI-driven turbulence. Once the initial background density profile $ρ_0$ is chosen, $\chi$ only depends on the thermal conductivity $κ$, and on the density fluctuations. To ensure that our conclusions are generic, two types of conduction profiles are explored in this paper: Spitzer (i.e. $κ \propto T^{\frac{5}{2}}$, noted Spit), for it is the most consistent with plasma physics, at least when kinetic effects arising in dilute plasma are not accounted for, and a type with constant diffusivity across radii (i.e. $κ \propto ρ$, noted Cdiff) for its simplicity. In both cases, we must specify the proportionality constant $κ_0$. For a Spitzer conduction profile, a natural choice appears: $κ_0 = 4.2$, which is the numerical value given by Eq. (13) for $T = 2T_{vir}$, further divided by $κ_{vir}$. For the constant diffusivity profile, there is no such obvious option. For further comparison, we introduced an additional dimensionless parameter $κ_{eff}$ that we define as the ratio between the volume-averaged conductivity in the solution

<table><caption>**Table 2.** Simulations in 2D and 3D with their varying parameters.</caption><thead><tr><th>Name</th><th>Resolution</th><th>κ-profile</th><th>κ<sub>0</sub></th><th>κ<sub>eff</sub></th></tr></thead><tbody><tr><td>3dSf4e0</td><td>256×128×512</td><td>Spit.</td><td>4.2</td><td>1</td></tr><tr><td>2dSf4e0/S0</td><td>256 × 512</td><td>Spit.</td><td>4.2</td><td>1</td></tr><tr><td>3dDf5e-2</td><td>256×128×512</td><td>Cdiff.</td><td>0.05</td><td>0.38</td></tr><tr><td>2dDf5e-2</td><td>256 × 512</td><td>Cdiff.</td><td>0.05</td><td>0.38</td></tr><tr><td>2dSf5e-1Mc</td><td>256 × 512</td><td>Spit.</td><td>0.5</td><td>0.12</td></tr><tr><td>3dDf1e-2/F0</td><td>512×256×1024</td><td>Cdiff.</td><td>0.01</td><td>0.075</td></tr><tr><td>2dDf1e-2</td><td>512 × 1024</td><td>Cdiff.</td><td>0.01</td><td>0.075</td></tr><tr><td>2dSf6e-2</td><td>512 × 1024</td><td>Spit.</td><td>0.06</td><td>0.014</td></tr><tr><td>2dSf3e-2Mc/S1</td><td>1024 × 2048</td><td>Spit.</td><td>0.03</td><td>0.0072</td></tr></tbody></table>

**Notes.** All runs start with a purely azimuthal magnetic field $B_0 = 10^{-4}$, and have $Pr = 0.01$, $Pm = 1$. Spit (resp. Cdiff) stands for conduction profile $κ \propto T^{5/2}$ (resp. $κ \propto ρ$).

domain (for any conduction type) and the volume average of the (full) Spitzer conductivity:

$$\kappa_{\text{eff}} = \frac{\langle \kappa \rangle_V}{\langle 4.2 T_0^{\frac{5}{2}} \rangle_V} \quad (15)$$

The Prandtl number $Pr = v/\chi$, and its magnetic counterpart $Pm = v/\eta$, are set to 0.01 and 1, respectively, in all simulations. The Brunt-Väisälä and the MTI frequencies, $N$ and $ω_T$ (two critical MTI timescales in the ICM), are enforced through the choice of the initial HSE, which is common to the main simulations of this study and reflects astrophysical ICM stratification regimes. We subsequently carry out an extensive parameter study, summarised in Table 2, with respect to the profile and intensity of thermal conductivity only.

### 3.3. Numerical methods, domains, and boundary conditions

We solve the system of Eqs. (7–10), normalised as described in Sect. 3.2, using IDEFIX, a new Godunov-type finite-volume code for astrophysical flows (Lesur et al. 2023). At each cell interface, the in-and-out hyperbolic fluxes are computed through the approximate resolution of the current Riemann problem
with the HLLD solver (which is a well-known compromise between numerical diffusion and stability; Miyoshi & Kusano 2005) after the respective left and right states have been systematically reconstructed with a second-order scheme (associated with a slope limiter in order for the numerical scheme to mimic the total-variation-diminishing physical nature of the B-MHD equations). The Braginskii heat flux is computed with a finite-difference-like scheme (Parrish & Stone 2005). We found that using a slope limiter on the anisotropic heat flux is mandatory to avoid episodic anti-diffusion (Sharma & Hammett 2007) in a few runs, especially those in 2D in which the thermal conductivity has a Spitzer dependency. In the runs that have the suffix Mc in Table 2, a monotonised central scheme ensures the monotonicity-preserving nature of the numerical operator accounting for anisotropic diffusion. Regarding the time-dependent part of the B-MHD equations, the hyperbolic terms are integrated thanks to a three-stage (third-order) Runge-Kutta algorithm, while the integration of parabolic operators is sped up using a Runge-Kutta-Legendre super time-stepping scheme (Meyer et al. 2014; Vaidya et al. 2017). The magnetic field is evolved thanks to a constrained transport algorithm (Fromang et al. 2006), associated with the HLLD upwind constrained transport averaging scheme for the electromotive force (Londrillo & Del Zanna 2004).

Next, we discuss the solution domain, alongside the boundary conditions and the numerical resolutions used in the simulations. We performed simulations in both two and three dimensions. In the former case, the B-MHD equations are solved in a curved domain $V^2$, which is a donut-like equatorial plane extending from 0.3 to 1.5 $R_{vir}$ (the usually expected position of the virial shock) in radius $r$, and from 0 to $2\pi$ in azimuthal angle $\varphi$; while in three dimensions, $V^3$ is a spherical wedge extending over the same range of radii and azimuthal angles, and from $\theta = \pi/4$ to $3\pi/4$ in the polar direction.

In two and three dimensions, the boundary conditions (BC) in the azimuthal direction are naturally periodic; while in three dimensions only, additional BCs are needed in the polar direction; we chose periodic BCs for all fields, making our 3D model semi-global. In both cases, however, BCs are still needed in the radial direction. As mentioned in Sect. 1, a numerical ICM atmosphere in which efficient MTI-driven turbulence is present can be quickly brought back to isothermality. Motivated by the non-marginal temperature gradients observed in astrophysical galaxy clusters, we indirectly impose Dirichlet BC on the temperature in the main simulations of this paper, by enforcing the pressure and the density to keep their initial HSE values in the ghost cells, so as to avoid this pathological behaviour. The consequence, on the temperature profile, of choosing inhomogeneous Dirichlet, rather than homogeneous Neumann, BCs is explored in detail in Sect. 4.3. Regarding the velocity, we use stress-free penetrative BC (i.e. homogeneous Neumann BC for both the normal and tangential components of the velocity field) at the inner edge of the domain, and a stress-free non-penetrative BC (i.e. homogeneous Dirichlet on the normal component) at the outer edge; such BCs seem more appropriate than penetrative BCs with respect to the physics of the accretion shock in astrophysical galaxy clusters. The impact of the velocity field BC on the overall dynamics of the simulation is quickly discussed in Sect. 4.3. At the inner and outer edges of the solution region, we impose that the tangential components of the magnetic field are zero (that is homogeneous Dirichlet BCs on the tangential components), while the radial component is automatically computed so as to satisfy the divergence-free nature of the field. Therefore, the magnetic field is purely radial (or null) at the inner and outer boundaries of

the solution region, consistent with the initial magnetic geometry Eq. (14). The uniform, non-adaptive resolutions of our 2D simulations range from 256 × 512 (i.e. 8.4 kpc in the radial direction) to 1024 × 2048, to better resolve the runs with the lowest values of thermal diffusivity, as expected from the scaling of the MTI injection length, Eq. (5). In 3D, the resolution of our simulations varies between 256 × 128 × 512 and 512 × 256 × 1024.

### 3.4. Energy budgets

We finally introduce the diagnostics used to characterise the turbulence generated by the MTI at saturation. The first quantities analysed are the averages of the Mach number $\langle M \rangle_V$ and those of the available potential $\langle E_A \rangle_V$, kinetic $\langle E_K \rangle_V$, and magnetic $\langle E_M \rangle_V$, volume energy densities in the solution domain $V$. The available potential energy (APE) is defined as the amount of internal and gravitational potential energy which can be reversibly converted into kinetic energy; its budget will turn out to be particularly useful in Sect. 4.1 to understand how the MTI saturates. The use of this physical quantity is borrowed from environmental fluid dynamics, where it was first introduced by Lorenz (1955), further developed by Holliday & McIntyre (1981), Andrews (1981), Winters et al. (1995), and is currently still intensively investigated by Tailleux (2009, 2013, 2018). Here, we simply define the APE volume density as

$$E_A = \frac{g_0^2}{2N^2} \frac{\delta \rho^2}{\langle \rho \rangle_{\Omega}}, \quad (16)$$

where we have used the local density fluctuation $\delta\rho = \rho - \langle\rho\rangle_{\Omega}$ with respect to the spherically averaged density $\langle\rho\rangle_{\Omega}$ over solid angles $\Omega = (\theta, \varphi)$ at a given radius and time. From now on, spherical averages, denoted $\langle\cdot\rangle_{\Omega}$, of the thermodynamic quantities will formally replace the initial background profiles, denoted $g_0$, in any qualitative discussion or quantitative diagnostic to come, because the background thermodynamic profiles are dynamically evolving through the course of the different simulations. To our knowledge, the notion of APE has never been directly mentioned in an astrophysical context, while for example a similar concept of available energy has been recently studied in the fusion plasma community (Helander 2017, 2020). Nevertheless, we stress that a quadratic APE, as given by the right-hand side of Eq. (16), has often appeared (but has never been formally identified as such) in the context of astrophysical stratified fluid dynamics, such as in studies of galaxy cluster turbulence (Perrone & Latter 2022a,b; Mohapatra et al. 2020, 2021), of wave physics in stellar radiative zones (Prat & Lignières 2014; Ahuir et al. 2021), or in planetary atmospheric dynamics (Dhouib et al. 2024); precisely because the APE is exactly the part of the internal and gravitational potential energy that is modelled and accessible to a stratified fluid in the Boussinesq and anelastic approximations (through the buoyancy equation). Therefore, the leading energetic balance at the core of Perrone & Latter (2022a,b) phenomenological picture of the MTI saturation mechanism is an APE budget in disguise. Additional details about the APE are presented in Appendix B. Here we just give its volume-averaged evolution equation in compressible B-MHD, assuming that no APE enters or leaves the domain and that the background profiles are constant in time:

$$\frac{d\langle E_A \rangle_V}{dt} = - \langle \mathcal{A} \rangle_V - \langle C \rangle_V + \langle \epsilon_A \rangle_V + \left\langle \frac{\delta T}{T} \Sigma : \nabla \mathbf{v} \right\rangle_V + \left\langle \frac{\delta T}{T} \eta \mu_0 j^2 \right\rangle_V, \quad (17)$$
where $\delta p = p - \langle p \rangle_{\Omega}$ and $\delta T = T - \langle T \rangle_{\Omega}$ are the local pressure and temperature fluctuations. On the right-hand side of Eq. (17), the first two terms quantify the amount of buoyant and compressible works, $\mathcal{A} = -g_0 \delta \rho \delta v_r$ and $C = -v \cdot \nabla \delta p$ respectively, that can reversibly convert APE into kinetic energy. The third term, $\varepsilon_A = q_B \cdot \nabla \left(\frac{\delta T}{T}\right)$, is the conversion rate between the background potential energy (BPE; see Appendix B for an introduction) and the APE; the latter is created when $\varepsilon_A \ge 0$. It can be further expanded into two contributions, $\varepsilon_i$ and $\varepsilon_k$, by decomposing the temperature into its spherical average and fluctuating parts:

$$ \varepsilon_i = -\kappa b_r \partial_r \langle T \rangle_{\Omega} \hat{b} \cdot \nabla \left( \frac{\delta T}{T} \right), \quad (18) $$

$$ \varepsilon_{\kappa} = -\kappa (\hat{b} \cdot \nabla \delta T) (\hat{b} \cdot \nabla \frac{\delta T}{T}). \quad (19) $$

The nature of the fourth and fifth viscous and resistive terms on the right-hand side of Eq. (17) are further discussed in Appendix B. They should be negligible on account of the regime of Prandtl and magnetic Prandtl numbers chosen in the simulations (Pr = 0.01, Pm = 1), and even more so for the resistive term in a realistic ICM since Pm $\gg$ 1 there (Schekochihin et al. 2005).

Next, to better characterise the energy transport within the simulation domain $\mathcal{V}$, we introduced the evolution equation for the spherically averaged total energy $E_I + E_K + E_M + E_G$, where $E_G = \rho\Phi$ is the gravitational potential energy:

$$ \frac{\partial}{\partial t} \langle E + E_G \rangle_{\Omega} + \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 [Q_r + \mathcal{H}_r + \mathcal{G}_r + \mathcal{K}_r + \mathcal{D}_r + \Pi_r] \right) = 0. \quad (20) $$

We refer the reader to Appendix A for a description of the spherically averaged kinetic $\mathcal{K}_r$, viscous $\mathcal{D}_r$, and Poynting $\Pi_r$ fluxes. We argue that because of the very hot and gravitationally bounded nature of the ICM plasma (implying $E_M, E_K \ll E_I, |E_G|$), these fluxes are subdominant with respect to the spherically averaged radial Braginskii heat flux $Q_r$, enthalpy flux $\mathcal{H}_r$, and gravitational flux $\mathcal{G}_r$:

$$ Q_r = \langle -\kappa (\hat{b} \cdot \nabla T) b_r \rangle_{\Omega}, \quad (21) $$

$$ \mathcal{H}_r = \frac{\mathcal{R}\gamma}{\gamma - 1} \langle \rho T \rangle_{\Omega}, \quad (22) $$

$$ \mathcal{G}_r = \Phi \langle \rho v_r \rangle_{\Omega}. \quad (23) $$

The radial Braginskii heat flux is the sum of the contributions from the average background temperature and from the temperature fluctuations:

$$ Q_{r,0} = \langle -\kappa b_r^2 \rangle_{\Omega} \partial_r \langle T \rangle_{\Omega}, \quad (24) $$

$$ Q_{r,1} = \langle -\kappa (\hat{b} \cdot \nabla \delta T) b_r \rangle_{\Omega}. \quad (25) $$

The enthalpy and gravitational fluxes can similarly be subdivided into advective and convective contributions:

$$ \mathcal{H}_{r,0} = \frac{\mathcal{R}\gamma}{\gamma - 1} \langle \rho T \rangle_{\Omega} \langle v_r \rangle_{\Omega}, \quad (26) $$

$$ \mathcal{H}_{r,1} = \frac{\mathcal{R}\gamma}{\gamma - 1} \left( \langle \rho \rangle_{\Omega} \langle \delta T \delta v_r \rangle_{\Omega} + \langle T \rangle_{\Omega} \langle \delta \rho \delta v_r \rangle_{\Omega} \right), \quad (27) $$

$$ \mathcal{G}_{r,0} = \Phi \langle \rho \rangle_{\Omega} \langle v_r \rangle_{\Omega}, \quad (28) $$

$$ \mathcal{G}_{r,1} = \Phi \langle \delta \rho \delta v_r \rangle_{\Omega}. \quad (29) $$

The advective fluxes $\mathcal{H}_{r,0}$, $\mathcal{G}_{r,0}$ are related to a global radial mass transfer (which is zero in absence of any bulk motion $\langle v_r \rangle_{\Omega}$),

whereas the convective fluxes $\mathcal{H}_{r,1} = \mathcal{H}_{r,\delta T} + \mathcal{H}_{r,\delta \rho}$, $\mathcal{G}_{r,1}$ are due to a positive correlation between the turbulent fluctuations of temperature, density, and radial velocity. Accordingly, we also monitor the volume-averaged level of the following cross-correlation function at saturation:

$$ \alpha (\delta v_r, \delta T) = \frac{\langle \delta v_r \delta T \rangle_{\nu}}{\sqrt{\langle \delta v_r^2 \rangle_{\nu} \langle \delta T^2 \rangle_{\nu}}}. \quad (30) $$

Finally, spectral diagnostics are used to assess both the convergence of the simulations and the spectral dynamics of MTI-driven turbulence. In 2D and 3D, we Fourier-analyse the azimuthal-$\varphi$ dependence, since spherical fields are naturally periodic in this direction, calling $m$ the associated angular degree. Additionally in 3D, we also do so for the polar-$\theta$ dependence ($l$ is the corresponding angular degree), as permitted by quasi-global $\theta$-periodic model used in this work. Therefore in 3D, the different fields are expanded as 2D plane waves $\propto \exp(i\ell\theta + im\varphi)$ (and 1D plane waves $\propto \exp(im\varphi)$ in 2D runs). We stress that, in this work, $l$ is the angular degree of a plane wave and not the order of a spherical harmonic. We do not expand fields on a spectral basis in the radial direction, but we analyse the spectra of the radially averaged angularly expanded fields instead. Given the highly stochastic nature of turbulence, we systematically accumulate as much statistics as possible for the sake of diagnostic convergence. So, in the next section where simulation results are presented, all diagnostics are time-averaged (during the saturation phase of developed turbulence) when not plotted as a function of time. We note, however, that such time averages are never made explicit in notations. All simulations are integrated during 200 dynamical times $t_{dyn}$, knowing that the turbulent turnover time at saturation is $\sim 10t_{dyn}$ in all simulations.

# 4. Simulation results

In this section we present the results of the 2D and 3D simulations of MTI-driven turbulence whose parameters are summarised in Table 2. We first focus on the time evolution of the volume-averaged energy densities and APE conversion rates of the fiducial 3D run F0. Then, we assess the structure and the intensity of the MTI flows, in all 2D and 3D simulations, using direct visualisation, global averages, and spectral diagnostics. Finally, we analyse the energy flux budget in the quasi-stationary saturated state of MTI-driven turbulence in the 3D run F0.

## 4.1. Evolution of the fiducial 3D MTI simulation

We first review the basic properties of the 3D MTI simulation F0, which has $\kappa_{eff} = 0.075$. We choose it as our fiducial run because of its intermediate regime of thermal conductivity, and the simple shape of its diffusivity profile (that is constant across radii, similar to all Cdiff runs). It will be later compared to the two 2D simulations, S0 and S1, with extreme values of thermal conductivity, and with a realistic Spitzer dependency, to emphasise the similar physics governing the dynamics of all of these runs, no matter the regime, and dependency, of the thermal conductivity. The time evolution of the different component contributions to the average kinetic $\langle E_K \rangle_{\nu}$ and magnetic $\langle E_M \rangle_{\nu}$ volume energy densities, and available potential energy (APE) $\langle E_A \rangle_{\nu}$, are shown on the top panel in Fig. 2.

All of these quantities quickly increase at the start of the simulation, as a result of MTI exponential growth. All along the simulation, the contribution of radial motions to the kinetic energy
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Key quantities of magneto-thermal turbulence as a function of time during the 3D run F0. Top: Volume-averaged time evolution of the available potential energy density $\langle E_A \rangle_V$ and of the different contributions to the kinetic $\langle E_K \rangle_V$ and magnetic $\langle E_M \rangle_V$ volume energy densities. Bottom: Volume-averaged time evolution of the APE injection rate $\langle \epsilon_i \rangle_V$, APE thermal dissipation rate $-\langle \epsilon_k \rangle_V$, and the (opposite amount of) reversible buoyancy work $-\langle \mathcal{A} \rangle_V$. Inset: Zoom-in of the highlighted region in which the additional quantity $-\langle \epsilon_k \rangle_V - \langle \mathcal{A} \rangle_V$ is plotted. Dashed lines are indicative of negative quantities whose sign has been switched for the purpose of visualisation.</figcaption></figure>

outweighs that of their horizontal counterparts, on account of the buoyant nature of the instability. Therefore, the MTI is able to drive energetic vertical motions despite the strong stable entropy stratification that tends to impede adiabatic fluid motions in the ICM (Mohapatra et al. 2020, 2021). As highlighted in the first column in Table 3, the volume-averaged Reynolds number of the different simulations, including the fiducial 3D F0 run under consideration here, defined as $Re = \ell_i \langle v_{rms} \rangle_V \langle \rho \rangle_V / \langle \mu \rangle_V$ (the injection length $\ell_i$ is quantified is the next subsection), is only moderate. Accordingly, the magnetic Reynolds number, defined as $Rm = RePm$, is always lower than the critical threshold above which the MTI-driven turbulence can operate a dynamo<fnref n="3" /> (Perrone & Latter 2022b). This is because of the $Pm = 1$ regime chosen here to avoid too constraining resolutions. Therefore,

we do not expect significant feedback from the magnetic field through the (bulk) Lorentz force in these simulations, given the very high plasma beta of the initial state and the absence of dynamo, even in our 3D runs. Consequently, despite its initial exponential growth due to the MTI field-line stretching, the magnetic energy is several orders of magnitude smaller than the kinetic energy after the linear phase. The anisotropy between the radial and horizontal contributions to the magnetic energy is even stronger than for the kinetic energy. The slow, constant decrease of magnetic energy is due to magnetic field being advected out of the domain (through the Poynting flux) by a thermal wind that is further described in Sect. 4.3. Such dynamics, however, only modestly affect the turbulence here since all components are decaying at the same rate, and in the B-MHD regime we considered, only the direction (and not the magnitude) of the magnetic field (through anisotropic heat flux) matters when the magnetic field stays dynamically negligible (which is the case here).

To understand MTI saturation, we focused on the time evolution of the different sink and source terms of APE, Eq. (17). In the bottom panel of Fig. 2, we display the time evolution of the average conversion rates, $\langle \epsilon_i \rangle_V$, $-\langle \epsilon_k \rangle_V$, along with the (opposite amount of) reversible buoyancy work $-\langle \mathcal{A} \rangle_V$. The first striking feature is the positive sign of $\langle \epsilon_i \rangle_V$. Middleton & Taylor (2020), Tailleux (2024) recently showed that, in oceanic physics, conversion from background potential energy (BPE) into APE can be associated with the growth of double-diffusive Ledoux-like instabilities. On the contrary, in incompressible Boussinesq single-component stratified fluid dynamics, APE can only be dissipated into BPE, and backward transferring from BPE to APE is prevented (Winters et al. 1995; Middleton & Taylor 2020). Therefore, such a positive conversion rate seems indicative of a thermodynamic process that harvests energy from the BPE and invests it into APE in the form of temperature and density fluctuations. The injected APE $\langle \epsilon_i \rangle_V$ is then distributed between the amount of APE dissipated back to the BPE $\langle \epsilon_k \rangle_V$ and the kinetic energy, through the reversible work of buoyancy $\langle \mathcal{A} \rangle_V$. Quantitatively speaking, the relation $\epsilon_i + \epsilon_k + \mathcal{A} = 0$ does not hold exactly (see the inset in Fig. 2), most certainly because of small APE boundary fluxes not accounted for in the current budget (see Appendix B). This suggests that the saturation mechanism proposed by Perrone & Latter (2022a,b) holds for global stratified MTI-driven turbulence developing in an astrophysically stratified ICM atmosphere as well, as this is precisely the balance on which their theory relies ($\epsilon_i \approx -\epsilon_k - \mathcal{A}$). We therefore expect the scaling laws for the MTI kinetic energy Eq. (6) and injection length Eq. (5) to remain valid in global MTI-driven turbulence. As a final note, we highlight that $\epsilon_i \approx -\epsilon_k$ is just a restatement of $q_B \approx 0$. Thus, we found that the MTI saturates by marginalising the Braginskii heat flux $q_B$, whereas Parrish & Stone (2005, 2007) argued that the MTI would saturate by shutting itself off through the marginalisation of the background temperature gradient $\partial_r \langle T \rangle_\Omega$.

## 4.2. Structure and strength of the flow

In this section we characterise the turbulent properties of MTI-induced turbulence, both in terms of strength and structure. On the top-left donut in Fig. 3, we show the velocity field induced by the MTI at saturation in the 3D simulation F0, where buoyant plumes characteristic of the 3D MTI (Perrone & Latter 2022b) are discernible. These turbulent eddies, which are preferentially elongated in the direction of gravity, reflect the anisotropy in the partition of kinetic energy between vertical and horizontal

<footnote n="3">³ Perrone & Latter (2022b) found a critical $Rm$ of about 35. Their definition, however, includes an additional $2\pi$ factor with respect to ours.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Visualisation of magneto-thermal turbulence in the simulated outskirts of galaxy clusters. Top left: Velocity field in the equatorial plane $(r, \varphi)$ of the 3D run F0 at saturation. Top right: Mach number (left half) and magnetic strength (right half) snapshots of the 2D run S0 at saturation. Bottom right: Same but for the 2D run S1. Bottom left: Temperature fluctuation maps and magnetic field lines (in black). Different 2D runs are shown in each sixth of the pie and are ranked (counter-clockwise, starting from the top-left part) by increasing effective thermal conductivity $\kappa_{\text{eff}}$, as in Table 2. In the latter plot, the colour map is scaled independently for each sixth of the pie. Reddish blobs are hotter than the spherical average, and blueish fluid elements are colder. The beige dashed line at $0.6R_{\text{vir}}$ on the top-right donut indicates the limit beyond which the ATHENA/X-IFU would collect less than $1 \times 10^6$ photons in the 0.2–12 keV range during a 1 Ms observation (Kempf et al. 2023, see the discussion later in Sect. 6.1), as computed prior to the recent X-IFU design-to-cost exercise (Barret et al. 2023).</figcaption></figure>

components, that was highlighted in the previous subsection. In all 2D and 3D simulation, the typical phenomenology of the MTI is respected: hot fluid elements rise while colder blobs sink, as quantified by the strong positive cross-correlation $\alpha \sim 0.5$ between the radial velocity and temperature fluctuations (last entry in Table 3). Additionally, the bottom-left donut in Fig. 3 shows that strong gradients of temperature fluctuations, typical of the MTI, are developing perpendicular to the magnetic field lines (that are superimposed as black lines on slices of different 2D runs).

As shown in Sect. 4.1, the MTI saturation mechanism of Perrone & Latter (2022a,b) theory, $\epsilon_i \approx -\epsilon_k \approx -\mathcal{A}$, essentially seems to hold in our global stratified models of MTI-driven tur-

bulence. A key prediction of the theory is then the scaling of the injection length, Eq. (5), and of the kinetic energy, Eq. (6), with the thermal diffusivity. Such a scaling could have important implications for the ICM given the large value of Spitzer conductivity in such plasma (if unsuppressed by kinetic instabilities). We therefore check both scaling laws, along with that for the APE, in global MTI-driven and stratified turbulence, on a volume-averaged basis for all runs in Table 2. To better assess the dependency of the flow strength and structure on the thermal diffusivity, we mostly rely on the 2D runs, which allow for a broader parameter sweep than in 3D. We start with both 2D runs S0 and S1, which only differ in the intensity of their thermal conductivities: they have respectively the highest and the lowest
<table><caption>**Table 3.** Key volume-averaged physical quantities at saturation for the 2D and 3D simulations of MTI-driven turbulence listed in Table 2.</caption><thead><tr><th>Run</th><th>Re</th><th>⟨M⟩<sub>V</sub></th><th>⟨E<sub>K,r</sub>⟩<sub>V</sub></th><th>⟨E<sub>K,h</sub>⟩<sub>V</sub></th><th>⟨E<sub>M,r</sub>⟩<sub>V</sub></th><th>⟨E<sub>M,h</sub>⟩<sub>V</sub></th><th>⟨v<sub>i</sub>⟩<sub>V</sub> / ⟨v<sub>i,rms</sub>⟩<sub>V</sub></th><th>⟨ε<sub>i</sub>⟩<sub>V</sub></th><th>−⟨ε<sub>k</sub>⟩<sub>V</sub></th><th>m<sub>i</sub></th><th>α (δv<sub>i</sub>, δT)</th></tr></thead><tbody><tr><td>3dSf4e0</td><td>100</td><td>0.14</td><td>0.0011</td><td>0.00058</td><td>4.1e−11</td><td>1.7e−11</td><td>0.61</td><td>0.0037</td><td>0.0028</td><td>3.1</td><td>0.59</td></tr><tr><td>2dSf4e0/S0</td><td>67</td><td>0.17</td><td>0.0013</td><td>0.00091</td><td>9.9e−12</td><td>5.2e−12</td><td>0.59</td><td>0.0056</td><td>0.0047</td><td>4.1</td><td>0.54</td></tr><tr><td>3dDf5e−2</td><td>170</td><td>0.13</td><td>0.001</td><td>0.0005</td><td>1.2e−09</td><td>5.7e−10</td><td>0.59</td><td>0.0017</td><td>0.0012</td><td>5.3</td><td>0.55</td></tr><tr><td>2dDf5e−2</td><td>140</td><td>0.15</td><td>0.0011</td><td>0.00071</td><td>7.2e−11</td><td>5.8e−11</td><td>0.6</td><td>0.002</td><td>0.0016</td><td>5.9</td><td>0.43</td></tr><tr><td>2dSf5e−1Mc</td><td>94</td><td>0.078</td><td>0.00021</td><td>0.00025</td><td>9.3e−10</td><td>1e−09</td><td>0.2</td><td>0.00075</td><td>0.00061</td><td>11.4</td><td>0.36</td></tr><tr><td>3dDf1e−2/F0</td><td>111</td><td>0.05</td><td>0.00017</td><td>0.0001</td><td>1.3e−07</td><td>6.8e−08</td><td>0.35</td><td>0.00034</td><td>0.00024</td><td>15.4</td><td>0.5</td></tr><tr><td>2dDf1e−2</td><td>119</td><td>0.065</td><td>0.00019</td><td>0.00026</td><td>3.4e−09</td><td>3.6e−09</td><td>0.27</td><td>0.0005</td><td>0.00041</td><td>12.8</td><td>0.29</td></tr><tr><td>2dSf6e−2</td><td>126</td><td>0.032</td><td>2.9e−05</td><td>3.5e−05</td><td>1.1e−08</td><td>1.3e−08</td><td>0.065</td><td>8.4e−05</td><td>6.6e−05</td><td>30.1</td><td>0.28</td></tr><tr><td>2dSf3e−2Mc/S1</td><td>123</td><td>0.023</td><td>1.6e−05</td><td>2e−05</td><td>2.1e−08</td><td>2.5e−08</td><td>0.048</td><td>4.4e−05</td><td>3.5e−05</td><td>43.4</td><td>0.25</td></tr></tbody></table>

**Notes.** The angular degree $m_i$, representative of the injection scale, is defined in Eq. (31). In 2D, $E_{K,h} = E_{K,φ}$ while in 3D, $E_{K,h} = E_{K,θ} + E_{K,φ}$, and similarly for the magnetic energy $E_M$. The Reynolds number is Re = $\ell_i \langle v_{rms} \rangle_V \langle \rho \rangle_V / \langle \mu \rangle_V$, from which we further define the magnetic Reynolds number $Rm = PmRe = Re$ and the Peclet number $Pe = PrRe = 0.01Re$.

ones. Despite sharing a similar phenomenology, MTI flows in S0 and S1 exhibit very different driving scales at saturation (which is easily seen when comparing both right donuts in Fig. 3). For instance, cluster-size flows are present in S0 (top-right donut), which corresponds to a realistic Spitzer conductivity case; while the typical turbulent length scale is much smaller in S1 (bottom-right donut), and progressively increases with the level of thermal conductivity (counter-clockwise bottom-left donut).

We can better quantify this diffusivity dependence using the spectral diagnostics defined in Sect. 3.4. In Fig. 4, we show the different radially averaged spectral energy densities $E_K(m)$, $E_M(m)$ and $E_A(m)$ for all runs in Table 2. In 3D, an additional average across the polar angular degree $l$ is performed: $E_K(m) = \sum_l E_K(m, l)/n_l$, where $n_l$ is the number of polar modes. These plots illustrate the phenomenological picture of turbulence, in which energy is distributed over a range of scales, but most of it resides at large scales. The slopes are quite different from classical $k^{-5/3}$ (or $k^{-3}$ in 2D) Kolmogorov laws. This is expected for MTI turbulence in the ICM since it is neither homogeneous nor isotropic. We further note that the kinetic and magnetic energies tend to reach smaller scales than the APE (i.e. the thermodynamic fluctuations), on account of the low thermal Prandtl number which makes the parallel thermal dissipation scale larger than the (isotropic) resistive and viscous scales.

We are now in a position to constrain the MTI integral scale $\ell_i$ using the azimuthal angular degree $m_i$ (in both 2D and 3D) deduced from the radially averaged volume kinetic energy spectral density, $E_K(m)$ in 2D and $E_K(l, m)$ in 3D:

$$m_i = \frac{\sum_m m E_K}{\sum_m E_K}, \quad \text{in 2D,} \qquad (31)$$

$$m_i = \frac{\sum_m m (\sum_l E_K)}{\sum_m (\sum_l E_K)}, \quad \text{in 3D.} \qquad (32)$$

This quantity roughly matches the order $m$ at which the kinetic spectrum peaks. In the top panel of Fig. 5, we show the relation between the global horizontal injection length, defined as $\ell_i = \pi (r_{in} + r_{out})/m_i$ (in both 2D and 3D) and the volume-averaged thermal diffusivity, for all simulations in Table 2. The scaling Eq. (5) is very well verified across two decades of variation in thermal diffusivity. From an astrophysical viewpoint, this means that the MTI alone drives cluster-size flows in the ICM, with unsuppressed realistic Spitzer conductivity ($\kappa_{eff} = 1$), as illustrated in the top-right donut in Fig. 3.

Next, we focus on the strength of the flow induced by MTI turbulence and compare it to the scaling Eq. (6). Again, this com-

parison is made in a global sense. On the middle panel in Fig. 5, we show the evolution of the volume-averaged squared velocity fluctuations as a function of the volume-averaged thermal diffusivity, for all simulations in Table 2. The scaling Eq. (6) is extremely well verified across two decades of thermal diffusivity, with a slight drop in the case of S0 and its 3D counterpart. The latter might be attributed to box-size effects that prevent the MTI from fully developing, as previously shown by McCourt et al. (2013) using local models. The present study shows that such effects could be astrophysically relevant. The validity of Eq. (6) suggests that the MTI could drive turbulent motions of about a few hundreds of kilometers per second in the outermost regions of galaxy clusters, which correspond to Mach numbers $M \sim 0.3$, as illustrated on the top-right donut in Fig. 3. Despite their potential strength, such MTI motions would be hardly detectable at the edges of galaxy clusters, on account on the ICM low emissivity there. The latter drops as $\propto \rho^2$. Beyond the beige dashed line at $0.6R_{vir}$ on the top-right donut, the future ATHENA/X-IFU observatory would collect less than $1 \times 10^6$ photons in the 0.2–12 keV range during a 1-Ms observation (Kempf et al. 2023), given the X-IFU specifications before the recent design-to-cost exercise (Barret et al. 2023). In Sect. 5.3, the associated induced level of non-thermal pressure support is further discussed, along with the possible caveats of this determination with our modelling approach.

The same exercise is then performed, on the bottom panel in Fig. 5, to assess the dependency of the APE (i.e. thermodynamic fluctuations) on the volume-averaged thermal diffusivity. Similarly to the case of the turbulent velocity, the APE seems to follow a power-law $E_A \propto \chi$ over two decades in thermal conductivity. This scaling slightly differs from that proposed by Perrone & Latter (2022b) ($E_A \propto \chi^{3/4}$). Our scaling is used in Sect. 5.2 to predict the level of convective flux that should be expected from MTI-driven turbulence.

Finally, directly comparing the global averages of the 3D simulations to their 2D counterparts (red and blue markers at the same $\langle \chi \rangle_V$ in Fig. 5) reveals a systematic bias, which is particularly visible in the case of the APE (which is lower in 3D than in 2D). Similarly to the APE, the kinetic energy of the 3D runs is always below than that of their 2D counterparts, while their typical injection length is usually slightly above that of the 2D runs. We think that this effect is purely geometric: the geometry and the dimensionality of the magnetic field vector must somehow differ between 2D and 3D simulations. This could for example explain the lesser efficiency of the Braginskii heat flux to enforce the isothermality required by the MTI to drive buoyancy
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Radially averaged spectra plotted against the azimuthal order *m* for all simulations in Table 2. Darker colours correspond to a higher *κ*<sub>eff</sub>. Dashed lines are for 3D runs, and regular lines are for 2D runs. Top row: Kinetic energy spectral density. Middle row: Magnetic energy spectral density. Bottom row: APE spectral density.</figcaption></figure>

(see the discussion in Sect. 5.1). Mathematically speaking, the augmented dimensionality simply induces different prefactors in the relation εᵢ ∝ -εₖ ∝ -A (through their dependence on the unit magnetic field vector **b**), which may explain the systematic biases observed in Fig. 5 between the 3D and 2D simulations.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** MTI scaling laws. Volume average of the turbulent injection length (top), the velocity fluctuation strength (middle), and specific APE (bottom), as a function of the volume-averaged thermal diffusivity, for all simulations in Table 2. Diamonds (respectively, circles) correspond to runs with Spitzer (respectively, constant diffusivity) conduction profile; while the colour blue (resp. red) stands for 2D (resp. 3D) runs.</figcaption></figure>

4.3. Boundary conditions, thermodynamic profiles, thermal
wind, and energy flux budget at saturation

Parrish & Stone (2005, 2007) first highlighted that the MTI has
the apparent ability to bring the initial temperature gradient
close to isothermality, especially when adiabatic boundaries are imposed (homogeneous Neumann BC on the temperature). In the design of our simulations, we thus tested several types of BC while paying special attention to the evolution of the background thermodynamic profiles with time. Accordingly, we found that imposing Neumann (zero-flux) BC on the temperature field leads to almost flat temperature distributions at the end of the simulations with the highest conductivities, or partial and still ongoing flattening in the case of the runs with the lowest conductivities. We argue, however, that this feature is only an indirect consequence of the MTI; the latter only creates a background diffusive heat flux by radially re-orienting the magnetic field lines, and does not marginalise the profile by its own convective heat transport. We note that isotropic turbulence would have the same effect, expect perhaps less efficiently so. Indeed, we found that the timescale of this isothermalisation process is simply governed by the thermal diffusion time. For instance, we performed similar runs in which we replaced Braginskii by isotropic heat flux. Such simulations exhibit a similar flattening of the temperature profile, and most importantly, with very similar timescales, while no turbulence is occurring. This behaviour was first seen in simulations by Parrish et al. (2008), but only later elucidated by Parrish et al. (2012b). Therefore, the flattening of the ICM temperature profile should not be seen as a specific property of the MTI saturation mechanism, since any isotropically tangled magnetic fields (which are most likely found in astrophysical galaxy clusters, contrary to purely isothermal horizontal magnetic field lines) of any turbulent origin would induce similar behaviours, to within a geometric suppression factor. As pointed out by Parrish et al. (2012b), the genuine problem of these kinds of simulations lies in the fact that they do not self-consistently resolve the physics of both the MTI and the accretion shock, which acts as the main driver of the temperature gradients in the ICM (McCourt et al. 2013). Imposing inhomogeneous Dirichlet BC on the temperature then provides a convenient and consistent way to provide our numerical ICM atmosphere access to a large energy reservoir that has been filled, all along the cluster dynamical hierarchical accretion history (which is not self-consistently modelled here), thanks to heating by the accretion shock. This is the path we took in the present work.

Nevertheless, even with such a BC, non-trivial saturation of the thermodynamic profiles still occurs in our simulations. In the top panel of Fig. 6, we show the leading radial fluxes transporting energy through the solution region of F0 after a quasi-stationary regime of MTI saturation is reached. The presence of a strong advective flux $H_{r,0}$, Eq. (26), can be somewhat surprising, but is actually easily explained. Our initial HSE state is not at thermal equilibrium in presence of any heat flux in the radial direction, but since the initial magnetic field lines are purely horizontal and the radial heat flux is shut off, it is nevertheless initially at rest. At the end of the linear phase, the opening of the magnetic field lines by the MTI in the radial direction induces a Braginskii heat flux in this direction. However, such a non-zero radial heat flux destabilises the initial profile by bringing it out of HSE and causes a subsequent uniform radial bulk motion, which we refer to as a thermal wind. Such dynamics are quite common in planetary atmospheric physics: upward thermal conduction of energy deposited deep into the atmosphere can drive a wind that tends to escape the planet's gravitational field (Johnson 2010). This model of atmospheric escape is usually referred to as the slow hydrodynamic escape (Volkov et al. 2011), and it is straightforward to explain, at least qualitatively. If a shell of stratified fluid heated from below is receiving more energy than it can thermally diffuse to cool down, it tends to expand. However, in

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Energy transport in the 3D simulation F0. Top: Dominant equilibrium between Braginskii heat, enthalpy, and gravitational energy fluxes, Eqs. (24), (26), and (28), in the energy flux budget, Eq. (20), at the quasi steady state. Bottom: Subdominant perturbed energy fluxes, Eqs. (25), (27), and (29). Dashed lines are indicative of negative quantities whose sign has been switched for the purpose of visualisation.</figcaption></figure>

doing so, the pressure stratification induces asymmetry and prevents the shell from expanding towards the inner region. Such an expansion takes then the form of an outwards wind, escaping the gravitational potential. In our simulations, the relative strength of the thermal wind (with respect to the MTI turbulence) is quantified through the ratio $\langle v_r \rangle_\nu / \langle v_{r,\text{rms}} \rangle_\nu$ in Table 3. The latter ratio increases with the thermal diffusivity because $\langle v_r \rangle_\nu \propto \chi$ (this result is empirical but physically in line with our interpretation of the wind origin), while for the MTI, $\langle v_{r,\text{rms}} \rangle_\nu \propto \chi^{1/2}$.

This bulk advective heat flux was already seen in the simulations of Parrish et al. (2008), and deemed to be astrophysically irrelevant in that same work on the basis that the dark matter profile should somehow be able to adjust to the change in the ICM structure so as to limit this fluid advection. We do agree with the fact that this thermal wind might only bear a limited astrophysical relevance, but not with the physical mechanism invoked to justify this: the dynamics of the dark matter could
hardly be set by that of the baryons, while the inverse is certainly
true on account to their respective mass budgets (e.g. gas slosh-
ing in galaxy clusters; Johnson et al. 2012; ZuHone et al. 2013).
Consequently and similarly to the possible flattening of the tem-
perature profile, we rather see this wind as a consequence of not
resolving, in isolated cluster simulations, the whole gravitational
structure formation flows, whose interplay with other physical
mechanisms would otherwise self-consistently set the thermo-
dynamic ICM gradients without necessarily resorting to a ther-
mal wind. However, we do not rule out possible short episodes
of thermal wind in the ICM, for example when AGN is feed-
ing back more energy in the cluster core than the ICM plasma
is able to thermally diffuse, although Bremsstrahlung cooling
might also be at play in this process. Additionally, we found the
structure of this wind to be critical for the shape of the thermo-
dynamic profiles at saturation, in the runs with the highest con-
ductivities κ<sub>eff</sub> (see the green thermodynamic profiles in Fig. 1
for the 2D run S0), but not that much for runs with moderate to
low thermal diffusivities (yellow curves on the same plot for the
3D run F0). The wind structure is also very dependent on the
BC imposed on the stratified thermodynamic quantities (either
inhomogeneous Dirichlet or a local HSE BC, as described by
Zingale et al. 2002), and on the velocity field (either homoge-
neous Dirichlet or Neumann), which may be indicative of the
limited relevance of this effect from an astrophysical viewpoint.
We provide an additional study of the influence of those BCs on
the wind, and on the MTI-driven turbulence itself more gener-
ally, in Appendix C.

In addition to heat, the mass transfer induced by this ther-
mal wind advects gravitational energy. The ordering E<sub>K</sub>, E<sub>M</sub> ≪
E<sub>I</sub>, |E<sub>G</sub>| found in the ICM, imposes the predominance of these
three fluxes (heat, enthalpy and gravitational) over the others
(kinetic, viscous, Poynting; see Appendix A). Interestingly, the
top panel in Fig. 6 shows that, in our 3D simulation F0, the
evolution of the background thermodynamic profiles saturates
by reaching a non-trivial dominant balance between the lead-
ing fluxes (purple dashed line). The same goes for the 2D run
S0 (not shown). In the latter case, the new equilibrium profile
(green curves in Fig. 1) looks similar to an intermediate state
between the initial atmosphere (blue curves) and a hypotheti-
cal ICM atmosphere that would be both at hydrostatic and ther-
mal equilibria (red curves) in presence of a radial background
heat flux with Spitzer conductivity. Accordingly, the outermost
shell of fluid exhibits a decreasing entropy profile, which vio-
lates Schwarzschild's criterion for adiabatic convective stabil-
ity. However, this fluid shell is very thin, and subject to both
anisotropic heat flux and viscosity. It is therefore very unclear
whether such a region would actually be unstable to standard
thermal convection, in addition to the MTI, or not, since nega-
tive entropy gradient is only a necessary condition for the onset
of standard thermal convection.

In the fiducial 3D run F0, a subdominant balance between
the perturbed fluxes induced by MTI-driven turbulence is par-
tially achieved. This balance is visible in the bottom panel in
Fig. 6. The negative perturbed conductive flux $Q_{r,1}$, Eq. (25), and
the convective gravitational flux $G_{r,1}$, Eq. (29), partially com-
pensate each other, while the perturbed convective heat fluxes
$H_{r,\delta T}$ and $H_{r,\delta\rho}$ cancel each other at each radii, on account of
the weakly compressible nature of the flow. The negative sign
of $Q_{r,1}$ translates the fact that, in MTI-driven turbulence, the
temperature fluctuations tends to fight the background tempera-
ture gradient, by rearranging according to Eqs. (2)–(4). We note,
however, that the convective energy fluxes, $H_{r,\delta T}$, $H_{r,\delta\rho}$, and $G_{r,1}$
are small with respect to the background conductive heat flux

$Q_{r,0}$. These results reinforce our argument that the isothermalisa-
tion of the temperature profile is not the dominant, specific, pro-
cess saturating the MTI; and that the latter is rather the marginal-
isation of the total Braginskii flux, as diagnosed through the APE
budget in Sect. 4.1. In Sect. 5.2, we provide an explanation for
these weak levels of convective energy fluxes, and discuss the
implication of this result for galaxy clusters.

**5. Implications for ICM energetics and dynamics**

In this section we show how the two different descriptions of the MTI at saturation, presented in Sect. 2, can be reconciled into a single diffusive MLT (DMLT). We then make use of this theory to predict the levels of convective fluxes and turbulent pressure support produced by the MTI in the ICM.

5.1. A diffusive mixing-length theory of the MTI

Here, we attempt to unify both MTI saturation descriptions that
have been previously reviewed in Sect. 2. A striking feature is
the close resemblance of Eqs. (2)–(4), especially if we assim-
ilate $\nabla$ to $\ell_i^{-1}$, the injection length of turbulence in the theory
of Perrone & Latter (2022a). Taking this path strongly suggests
that, for the MLT to correctly describe the MTI at saturation,
the mixing length $\ell_m$ should be chosen as the injection length $\ell_i$.
When doing so, Eq. (3) is necessarily equivalent to Eq. (6), since
$\omega_T$ is taken to be the natural time scale of MTI-driven turbulence
in the theory of Perrone & Latter (2022a,b). Whether $\omega_T$ and $\ell_i$
(and not the same quantities multiplied by some power of $N/\omega_T$
for instance) are exactly the right quantities to be regarded as
the time scale and the mixing length of such turbulence cannot
be decided within the current work, because no parameter study
varying independently $\omega_T$ and $N$ was carried out. However, this
caveat is only a theoretical limitation, which has very limited
to no impact on any astrophysical conclusions since those time
scales are very similar (and thus $N/\omega_T \sim 1$) in galaxy clusters
and for our initial HSE.

Actually, another and perhaps better physically motivated
reason why $\ell_i$, to within any reasonable power of $N/\omega_T$, should
be taken as the mixing length $\ell_m$ goes as follows. Using a Taylor
expansion to write Eq. (2), we implicitly assumed that the pertur-
bation was perfectly isothermal. In this framework, the mixing
length $\ell_m$ should therefore be understood as the typical length
over which any perturbation can stay isothermal with respect to
the temperature of its initial magnetic field line, in the cartoon
picture of the MTI at least. Given some typical time scale $\omega^{-1}$ of
MTI-driven turbulence (whose precise value does not matter for
the current discussion, it should be quite close to $\omega_T$ and $N$ in the
ICM anyway), the typical length over which a fluid element stays
thermally connected to the radius at which it started moving, dur-
ing the time scale $\omega^{-1}$ of the motion, must be proportional to the
conduction length, defined as $\ell_\chi = (\chi/\omega)^{1/2}$, to within possible
geometrical prefactors. Such a choice for the mixing length $\ell_m$
embeds the diffusive nature of the MTI in a MLT, which is oth-
erwise diffusion free and would therefore misrepresent the basic
physics of the instability. We refer to this approach as a diffu-
sive mixing-length theory (DMLT), which, satisfactorily, leads
to similar dependencies of the MTI injection length and veloc-
ity fluctuation strength on the thermal diffusivity (by comparing
Eq. (3) to Eq. (6) when $\ell_m \propto \ell_\chi$).

Similarly, we can compare the scaling of the temperature
fluctuation intensity derived from both the DMLT and the work
of Perrone & Latter (2022a,b). In the case of the DMLT, we find
$(\delta T/T_0)^2 \propto \chi$ rather than $(\delta T/T_0)^2 \propto \chi^{3/4}$. The former scaling
is also seemingly preferred in our simulations (bottom panel in Fig. 5). This might be so far the only example of divergent result between our work and Perrone & Latter (2022a,b), in which only the gradient of the APE is constrained thanks to phenomenological arguments: in their theory, the scaling of the APE is therefore numerically measured and not theoretically derived. Similarly to the issue of the precise MTI-driven time scale $\omega^{-1}$ and mixing length $\ell_m$, this discrepancy with our work is only of secondary importance regarding astrophysical conclusions. We therefore propose the three following scalings:

$$ \ell_i \approx 4.6 \left(\frac{\chi}{\omega}\right)^{0.5}, \qquad (33) $$

$$ v^2 \approx 0.13 \omega \chi, \qquad (34) $$

$$ \left( g_0 \frac{\delta T}{T_0} \right)^2 \approx 0.03 \omega^3 \chi, \qquad (35) $$

where $\omega$ is some typical frequency of MTI turbulent motions (taken equal to 1 Gyr$^{-1}$ here), which is likely to be a reasonable mix of $\omega_T$ and $N$. In the next sections, we apply these equations derived from the DMLT of MTI saturation to assess both the amount of convective fluxes and non-thermal turbulent pressure support that can be expected from the MTI in the outskirts of galaxy clusters, thereby revisiting the estimates of Parrish et al. (2012b) and McCourt et al. (2013).

## 5.2. Convective fluxes carried by the MTI

Before assessing the convective efficiency of the MTI, a specific issue first needs to be tackled. When previous authors referred to the convective flux induced by the MTI, they systematically meant the quantity named $H_{r,\delta T} \propto \langle \delta T \delta v_r \rangle_{\Omega}$ in the current work, which is the expression of the convective flux usually found in stellar convection (Hurlburt et al. 1984; Cattaneo et al. 1991; Porter 2000). But in all generality, two other convective fluxes arise from the strong cross-correlation between temperature, density and radial velocity fluctuations, namely $H_{r,\delta\rho} \propto \langle \delta\rho\delta v_r \rangle_{\Omega}$ and $G_{r,1}$, Eq. (29), and those must also be accounted for in any energy flux budget. The latter two fluxes are obviously proportional to each other, but they are related through another and more subtle relation too, and only cancel each other for adiabatic stratification. Indeed, the combination of the first and second principles of thermodynamics for an atmosphere at HSE can be rewritten (Hurlburt et al. 1984; Rieutord 2015) as

$$ \frac{\mathcal{R}\gamma}{\gamma-1} \frac{d\langle T \rangle_{\Omega}}{dr} + \frac{d\Phi}{dr} = \langle T \rangle_{\Omega} \frac{d\langle S \rangle_{\Omega}}{dr}, \qquad (36) $$

and therefore

$$ H_{r,\delta\rho} + G_{r,1} = \left[ \left( \frac{\mathcal{R}\gamma}{\gamma-1} \langle T \rangle_{\Omega}^* + \Phi^* \right) + \int_{r^*}^{r} \langle T \rangle_{\Omega} \frac{d\langle S \rangle_{\Omega}}{dr} \right] \langle \delta\rho \delta v_r \rangle_{\Omega}, \qquad (37) $$

where $r^*$ is some reference radius (and $\langle T \rangle_{\Omega}^*, \Phi^*$ are the background temperature and gravitational potential at this radius). A natural choice for the latter radius is that at which the constant term $(\frac{\mathcal{R}\gamma}{\gamma-1}\langle T \rangle^* + \Phi^*)$ cancels (which is always possible to find because the gravitational potential is only defined up to a constant). From Eq. (37), it is clear that $H_{r,\delta\rho} \approx -G_{r,1}$ in the limit of an isentropic background stratification (which is the natural endpoint of saturation in standard thermal convection such as that encountered in stellar convection). In galaxy clusters, however, the background entropy profile is not adiabatic because it is

imposed, at first order, by the dynamical hierarchical accretion history, and not by standard thermal convection. Then, there is no simple relation between $H_{r,\delta\rho}$ and $G_{r,1}$. Nevertheless, in the limit of mildly compressible MTI flows, we still have $H_{r,\delta T} \approx -H_{r,\delta\rho}$ since pressure fluctuations are small (see for example the bottom panel in Fig. 6). The appropriate turbulent energy transport term to consider when referring to the convective flux induced by mildly compressible MTI turbulence in this context is therefore $G_{r,1}$ rather than $H_{r,\delta T}$. Next, we assess its strength with respect to the background heat flux.

By essence, the DMLT introduced in Sect. 5.1 to describe the large-scale saturation of the MTI is well-suited to calculate the induced level of convective flux. Indeed, convection-driven turbulence (to which the MLT usually applies) exhibits a significant cross-correlation between temperature and velocity fluctuations, and so does the MTI for which we numerically measured $\alpha(\delta v_r, \delta T) \approx 0.5$ (see Table 3), on account of its buoyant nature. This exercise was previously performed by McCourt et al. (2013), but with an incomplete understanding of the convective fluxes and MLT associated with MTI-driven turbulence. Using Eqs. (34)–(35), we rederive the ratio of the convective flux of gravitational energy, Eq. (29), to the background Braginskii conductive heat flux, Eq. (24), as

$$ \frac{G_{r,1}}{Q_{r,0}} = 7\% \left( \frac{\alpha(\delta v_r, \delta T)}{0.5} \right) \left( \frac{b_r^2}{0.5} \right)^{-1} \left( \frac{\delta v_r^2 / \delta v^2}{0.5} \right)^{1/2} \left( \frac{n_0 / n_e}{0.5} \right) \left( \frac{\omega^2}{\omega_T^2} \right) \times \left( \frac{k_B T}{2 \text{ keV}} \right)^{-1} \left( \frac{\mu m_H \Phi}{2 \text{ keV}} \right). \qquad (38) $$

This result needs to be considered with care. Indeed, Eq. (38) only suggests that the convective gravitational flux induced by the MTI in the ICM is very weak with respect to the background heat flux, $Q_{r,0}$, contrary to previous claims regarding the convective heat flux (Parrish & Stone 2005, 2007; McCourt et al. 2011). However, MTI-driven turbulence is known to significantly open the magnetic field lines in the radial direction. This radially biased magnetic field offers a stronger channel for the background heat flux to flow within the ICM periphery, in comparison to an isotropically tangled magnetic field, and must to some extent be attributed to the MTI itself. This ordering between the convective gravitational and conductive heat fluxes is specific to the MTI: classical convection is usually more efficient at transporting energy than conduction. In the Sun for instance, the convective envelope is unstable to thermal convection precisely because the effective conduction of heat<fnref n="4" /> is too low to efficiently diffuse all the luminosity coming from the radiative zone, and other non-linear modes of energy transfer such as convection naturally arise. We stress that this is definitely not the situation of the MTI in the ICM. Unlike standard thermal convection, the MTI needs fast parallel heat conduction to develop in the first place and to thrive at saturation, as shown by Perrone & Latter (2022a,b). With this in mind, Eq. (38) is only a natural consequence of the MTI diffusive nature and should thus not be surprising.

## 5.3. Non-thermal pressure support generated by the MTI

We turn our attention to the non-thermal pressure support induced by the MTI in the outermost ICM. The scaling Eq. (34),

<footnote n="4">⁴ In the Sun, an effective thermal conduction is due to the radiative diffusion.</footnote>
equivalently described by the DMLT or by the balance mechanism between the APE leading injection and dissipation rates (Perrone & Latter 2022a,b), offers a proxy to assess the ratio of non-thermal to thermal pressures, given thermodynamic profiles, for typical MTI flows in the ICM. Using typical values for the outermost regions, $r \sim 1.5R_{vir}$, leads to

$$ \frac{p_{\text{turb}}}{p_0} = 14\% \left( \frac{p_0}{1 \times 10^{-14} \text{ g/cm/s}^2} \right)^{-1} \left( \frac{\omega}{1 \text{ Gyr}^{-1}} \right) \times \left( \frac{\kappa}{1 \times 10^{11} \text{ g/cm/s}^3/\text{K}} \right) \left( \frac{n_0/n_e}{0.5} \right). \quad (39) $$

Interestingly, this numerical value is non-negligible and is an indicator of moderate Mach numbers possibly induced by the MTI close to the virial shock. Nevertheless, we caution that this ratio is an increasing function of the radius. For example, it quickly drops to 5% for ICM values taken at $r \sim R_{vir}$, and it keeps doing so towards the innermost and denser regions of ICM haloes. The radial dependence of Eq. (39) should therefore carefully be accounted for before drawing definitive conclusions about the ability of the MTI to significantly bias the mass inferred from X-ray and SZ observations combined with the HSE assumption; especially since the highest values of this ratio are found in the least dense ICM regions, which weigh little in the mass budget of galaxy clusters. Another caveat to this result is the weakly to almost non-collisional nature of the ICM plasma at such large radii. For instance, well-known kinetic instabilities (Schekochihin et al. 2005, 2010; Drake et al. 2021) could efficiently suppress the thermal conductivity in weakly collisional plasmas and therefore also the strength of MTI-induced turbulence according to Eqs. (34)–(35). Their interplay with the MTI is briefly discussed in Sect. 6.2.

## 6. Summary and conclusions

In this paper, we have studied turbulence induced in a global spherical stratified ICM atmosphere by the saturation of the magneto-thermal instability (MTI), a magneto-buoyant instability likely to impact the internal dynamics and energetics of galaxy cluster outskirts. This work is specifically aimed at bridging the gap between idealised simulations of MTI-driven turbulence in local and global ICM models (Parrish & Stone 2007; Parrish et al. 2008; McCourt et al. 2011; Parrish et al. 2012b; Perrone & Latter 2022a,b) and large-scale simulations of structure formation through gravitational collapse including anisotropic thermal conductivity (Ruszkowski et al. 2011). The main results are threefold:

- The MTI saturation mechanism introduced by Perrone & Latter (2022a,b) under the Boussinesq approximation is found to result from a balance between the APE leading injection and dissipation rates. This balance still holds for strong stable entropy stratification and collisional regimes of unsuppressed parallel thermal conduction. A positive APE injection rate is only possible because the MTI saturates by marginalising the anisotropic Braginskii heat flux. We emphasise that this saturation mechanism does not imply the marginalisation of the background temperature gradient itself by the MTI, which was previously suggested to be the natural endpoint of this instability (Parrish & Stone 2005, 2007).

- The mixing-length theory (MLT) proposed by Parrish et al. (2012b), McCourt et al. (2013) to describe the large-scale dynamics of the MTI seems unrelated to the previous local

saturation mechanism of Perrone & Latter (2022a,b). Both can be made to be consistent with each other by choosing the mixing length to be proportional to the conduction length $\ell_{\chi} \propto \sqrt{\chi}$. This is a way to re-introduce the diffusive nature of this instability in a MLT, which remains otherwise essentially diffusion-free and thus cannot account for the physics of the MTI. We refer to this approach as a diffusive MLT (DMLT). However, since $\ell_{\chi} \sim H_T \sim H_p$ in the ICM anyway for a complete Spitzer conductivity, choosing $\ell_m \sim H_p$ over $\ell_{\chi}$ does not significantly impact the numerical predictions. This is the reason why the previous estimate of the convective over conductive energy fluxes ratio of a few percent by McCourt et al. (2013) is well in line with our estimate, Eq. (38), though we showed that the right convective flux to consider is $G_{r,1}$, Eq. (29), rather than $H_{r,\delta T}$.

- Finally, we confirmed a clear dependency of the turbulent injection length and strength on the thermal conductivity (over more than two decades of conductivity) that can be inferred from the saturation mechanism of Perrone & Latter (2022a,b). This dependency is well incorporated by our DMLT too. For instance, in the 2D simulation with a complete Spitzer conductivity, the MTI drives cluster-size motions at Mach numbers $M \sim 0.3$.

We summarise the astrophysical implications of these results in the next subsection.

### 6.1. Implications for the ICM

We made use of these results and the scaling Eq. (34) especially in order to deduce the levels of convective flux and non-thermal turbulent pressure support associated with MTI-driven turbulence in the ICM periphery. We found, both theoretically and numerically, that in standard B-MHD with no suppression of conduction, the MTI can sustain mildly compressible turbulence with a typical length scale of the size of the cluster, with Mach numbers as high as $M \sim 0.3$ in the outermost regions. In this case, the ratio of the non-thermal to thermal pressures $p_{\text{turb}}/p_0 = \gamma M^2$ can reach values up to 15% in regions close to the virial shock (see Eq. (39) for a precise estimate). We stress, however, that this ratio is a quickly increasing function of the radius mostly because the background pressure $p_0$ is exponentially decreasing. The radial variation of the MTI dynamical properties was discussed by Kempf et al. (2023), who assessed the possible detectability of this instability in a Perseus-like cluster at $r \sim 0.25R_{vir}$ with the future X-ray spectrometer X-IFU (Barret et al. 2023) on board ATHENA. In Kempf et al. (2023), we adopted a local radius-by-radius approach to the MTI dynamical properties. However, the present global study shows that in case of a realistic Spitzer conductivity, the MTI can drive buoyant plumes extending over the full ICM outskirts. This could possibly slightly boost the estimates of the velocity signal and of the injection length expected from the innermost ICM regions that are at the same time unstable to the MTI and available for X-ray observation with the ATHENA/X-IFU (inside the area delimited by the beige dashed line on the top-right donut in Fig. 3). This is why directly applying Eq. (39) locally at a given radius may not be justified with the highest values of thermal conductivity as far as the assessment of the turbulent pressure support is concerned.

Another application of Eq. (35) is for the assessment of the convective energy transport induced by the MTI. We showed that the ratio between the convective gravitational and the background Braginskii heat fluxes hardly exceeds 7%, Eq. (38), which is pretty low for a physical process that is often referred to as convection, or as a convective instability. Whether the lat-
ter denomination is justified depends on the meaning we give it. Does it refer to the ability of the induced-turbulence to carry significant convective flux (and if so, are we speaking about transport of heat or gravitational energy, or both?) or to the efficient coupling between temperature and radial velocity fluctuations? In standard thermal convection, such as in the Sun, both are true, and the energy carried by the convective flux is a mix between the gravitational and thermal components. However, because of the diffusive nature of the MTI, which needs fast parallel conduction to trigger turbulence, MTI-driven convection is not as efficient as conduction at transporting heat across the ICM. Referring to the MTI as a convective instability might therefore be misleading in this respect, and emphasising the buoyant nature of the instability instead seems more appropriate.

## 6.2. Caveats and future perspectives for ICM studies

In this section, we bring to light other possible applications of the MTI in an astrophysical context along with the associated caveats. First and most importantly, similar further studies about the interplay between the MTI and other forms of turbulence driven by different processes are needed to draw definitive conclusions about the actual impact of the MTI on the dynamics and physics of galaxy clusters. Second, there is also a strong need to understand the interaction of this fluid instability with possible harmful effects of kinetic micro-instabilities occurring in high-beta plasmas, which we completely overlooked in this paper, in order to assess whether this instability is relevant at all in ICM outskirts. This is a research area under intense investigation. For example, Berlok et al. (2021) demonstrated the robustness of the MTI against mirror modes suppressing the heat conductivity, while Perrone et al. (2024a,b) recently shed light on the possibly deadly non-linear feedback of whistler modes on the MTI, which strongly depends on the ability of the latter to sustain sufficient levels of magnetic field.

Speaking of which, a promising area of application for MTI studies of the kind conducted here lies in the possible dynamo effect that the MTI may drive. Whether the MTI can, through a small-scale fluctuation dynamo process, bring the very weak primordial magnetic fields up to equipartition with the kinetic energy, as observed in close clusters, is an outstanding question that we plan to investigate further in future work. Here, a thorough study of realistic magnetic field amplification by MTI-driven turbulence was prevented because of the $Pm = 1$ regime that was chosen and the subsequent moderate magnetic Reynolds numbers $Rm = O(100)$ of the intermediate resolutions of the few 3D runs presented. We further note that the choice of isotropic over Braginskii viscosity must certainly modify the turbulent dynamics at the viscous scale, on which fluctuation dynamo critically depends (Rincon 2019). Braginskii viscosity should therefore be carefully included in future dynamo works.

## 6.3. Possible astrophysical applications of the APE

In this paper, the concept of available potential energy (APE) has been introduced as a useful tool to understand the physics and saturation of the MTI in the ICM. While a complete and accurate description of this quantity in galaxy clusters is far beyond the scope of the current study (further details are still given in Appendix B), we argue that the APE is of primary importance not only for the ICM internal dynamics but also for the observation of these dynamics in X-ray. Indeed, pending the future generation of X-ray spectrometers such as Resolve (Ishisaki et al. 2018) or the X-IFU (Barret et al. 2023), indirect measurements

of the ICM kinematics are only achievable through the observation of X-ray surface brightness fluctuations (Zhuravleva et al. 2015; Dupourqué et al. 2023). In the ICM, density perturbations with respect to the background profile are usually assumed to be directly related to the velocity fluctuations through their respective spectra (Zhuravleva et al. 2014). In this methodology, the interplay between kinetic and available potential energies is clearly at stake and is actually taken to be scale independent. This assumption has been checked to some extent in hydrodynamic simulations of forced stratified turbulence with conditions representative of the ICM (Gaspari et al. 2014; Mohapatra et al. 2021). However, it remains to be tested in the presence of magneto-buoyant instabilities such as the MTI or the HBI since we demonstrated that in MTI-induced turbulence, kinetic energy is driven by, rather than driving, density fluctuations (i.e. APE).

We argue that this framework could also be applied to stellar, solar, and planetary atmospheric dynamics to further understand the energy transport through, possibly double-diffusive, convective fluxes within those astrophysical objects thanks to the fully compressible equations of fluid dynamics. As in the case of the MTI, the saturation of double-diffusive instabilities might be similarly controlled by the injection rate of APE (more details are given in Appendix B). In the same manner, the energetics of the $\kappa$-mechanism, which is responsible for the pulsation of Cepheid stars (Eddington 1917), that are used as standard candidates to build the cosmic distance ladder (Riess et al. 2019) seem to be ruled by the conversion of BPE into APE (though not formally interpreted as such; Gastine & Dintrans 2008). Besides, Tailleux & Dubos (2024) recently proposed static energy asymptotics as a method to better constrain the energetics of the different classes of Boussinesq and anelastic approximations that are used in the ocean and atmospheric modelling community. We argue that similar models used in astrophysics, such as the adiabatic anelastic approximation used in stellar physics (Brun et al. 2004; Miesch 2005), in protoneutron star dynamics (Reboul-Salze et al. 2022), and in planetary interiors (Gastine et al. 2012; Gastine & Wicht 2012, 2021), could benefit from the same treatment to improve our understanding of their energetics, thus building on the early insights from Braginsky & Roberts (1995).

**Acknowledgements.** The authors thank the anonymous referee for valuable comments which helped improve the clarity of this paper. The work presented here benefited from fruitful exchanges with Nicolas Clerc, Jean-Baptiste Durive, Thomas Guillet, Christopher Howland, Laurène Jouve, Henrik Latter, Geofroy Lesur, Leo Middleton, Yves Morel, Lorenzo Perrone, Michel Rieutord, and Remi Tailleux. The authors are grateful to R. Tailleux for his careful reading of Appendix B, and for bringing the relevant notion of extrinsic state function to their attention. J.M.K. has been partially supported through the grant EUR TESS N°ANR-18-EURE-0018 in the framework of the Programme des Investissements d'Avenir. This work was granted access to the HPC resources of CALMIP under the allocation 2023/2024-P16006, and to those of CINES under the allocations 2024-A0160415090 made by GENCI. We are grateful for the work of the associated support teams. J.M.K: Ce travail est dédié à la mémoire de mon défunt professeur Jean-Pierre Simond, dont le sens physique a beaucoup apporté à ma formation scientifique.

## References

* Abbett, W. P., Beaver, M., Davids, B., et al. 1997, ApJ, 480, 395

* Ahuir, J., Mathis, S., & Amard, L. 2021, A&A, 651, A3

* Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, A&ARv, 49, 409

* Andrews, D. G. 1981, J. Fluid Mech., 107, 227

* Angelinelli, M., Vazza, F., Giocoli, C., et al. 2020, MNRAS, 495, 864

* Balbus, S. A. 2000, ApJ, 534, 420

* Balbus, S. A. 2001, ApJ, 562, 909

* Barret, D., Albouys, V., den Herder, J.-W., et al. 2023, Exp. Astron., 55, 373
Berlok, T., Quataert, E., Pessah, M. E., & Pfrommer, C. 2021, MNRAS, 504, 3435

Böhm-Vitense, E. 1958, Z. Astrophys., 46, 108

Braginskii, S. I. 1965, Rev. Plasma Phys., 1, 205

Braginsky, S. I., & Roberts, P. H. 1995, Geo. Astro. Fluid Dyn., 79, 1

Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073

Canuto, V. M. 1990, A&A, 227, 282

Cattaneo, F., Brummell, N. H., Toomre, J., Malagoli, A., & Hurlburt, N. E. 1991, ApJ, 370, 282

Cavagnolo, K. W., Donahue, M., Voit, G. M., & Sun, M. 2009, ApJS, 182, 12

Cengel, Y. A., Boles, M. A., & Kanoğlu, M. 2011, Thermodynamics: An Engineering Approach (New York: McGraw-hill), 5

Chan, K. L., & Sofia, S. 1989, ApJ, 336, 1022

Chandran, B. D. G., & Maron, J. L. 2004, ApJ, 602, 170

Dhouib, H., Baruteau, C., Mathis, S., et al. 2024, A&A, 682, A85

Drake, J. F., Pfrommer, C., Reynolds, C. S., et al. 2021, ApJ, 923, 245

Dupourqué, S., Clerc, N., Pointecouteau, E., et al. 2023, A&A, 673, A91

Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40

Eddington, A. S. 1917, The Observatory, 40, 290

Ferrari, C., Govoni, F., Schindler, S., Bykov, A. M., & Rephaeli, Y. 2008, Space Sci. Rev., 134, 93

Fraser, A. E., Reifenstein, S. A., & Garaud, P. 2024, ApJ, 964, 184

Fromang, S., Hennebelle, P., & Teissier, R. 2006, A&A, 457, 371

Garaud, P. 2018, Annu. Rev. Fluid Mech., 50, 275

Gaspari, M., Churazov, E., Nagai, D., Lau, E. T., & Zhuravleva, I. 2014, A&A, 569, A67

Gastine, T., & Dintrans, B. 2008, A&A, 484, 29

Gastine, T., & Wicht, J. 2012, Icarus, 219, 428

Gastine, T., & Wicht, J. 2021, Icarus, 368, 114514

Gastine, T., Duarte, L., & Wicht, J. 2012, A&A, 546, A19

Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, A&A, 621, A41

Govoni, F., & Feretti, L. 2004, Int. J. Mod. Phys. D, 13, 1549

Helander, P. 2017, J. Plasma Phys., 83, 715830401

Helander, P. 2020, J. Plasma Phys., 86, 905860201

Holliday, D., & Mcintyre, M. E. 1981, J. Fluid Mech., 107, 221

Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1984, ApJ, 282, 557

Ishisaki, Y., Ezoe, Y., Yamada, S., et al. 2018, J. Low Temp. Phys., 193, 991

Johnson, R. E. 2010, ApJ, 716, 1573

Johnson, R. E., Zuhone, J., Jones, C., Forman, W. R., & Markevitch, M. 2012, ApJ, 751, 95

Kempf, J. M., Rincon, F., & Clerc, N. 2023, A&A, 680, A24

Komarov, S., Schekochihin, A. A., Churazov, E., & Spitkovsky, A. 2018, J. Plasma Phys., 84, 905840305

Kucharski, F. 1997, Quarterly J. R. Meteorol. Soc., 123, 2141

Kunz, M. W., Bogdanović, T., Reynolds, C. S., & Stone, J. M. 2012, ApJ, 754, 122

Leccardi, A., & Molendi, S. 2008, A&A, 486, 359

Lesur, G. R. J., Baghdadi, S., Wafflard-Fernandez, G., et al. 2023, A&A, 677, A9

Londrillo, P., & Del Zanna, L. 2004, J. Comp. Phys., 195, 17

Lorenz, E. N. 1955, Tellus A, 7, 157

McBride, J., Fakhouri, O., & Ma, C.-P. 2009, MNRAS, 398, 1858

McCourt, M., Parrish, I. J., Sharma, P., & Quataert, E. 2011, MNRAS, 413, 1295

McCourt, M., Quataert, E., & Parrish, I. J. 2013, MNRAS, 432, 404

Meyer, C. D., Balsara, D. S., & Aslam, T. D. 2014, J. Comp. Phys., 257, 594

Middleton, L., & Taylor, J. R. 2020, J. Fluid Mech., 893, R3

Middleton, L., Fine, E. C., MacKinnon, J. A., Alford, M. H., & Taylor, J. R. 2021, Geo. Res. Lett., 48, e2021GL092779

Miesch, M. S. 2005, Liv. Rev. Solar Phys., 2, 1614

Miyoshi, T., & Kusano, K. 2005, J. Comp. Phys., 208, 315

Mohapatra, R., Federrath, C., & Sharma, P. 2020, MNRAS, 493, 5838

Mohapatra, R., Federrath, C., & Sharma, P. 2021, MNRAS, 500, 5072

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Parrish, I. J., & Stone, J. M. 2005, ApJ, 633, 334

Parrish, I. J., & Stone, J. M. 2007, ApJ, 664, 135

Parrish, I. J., Stone, J. M., & Lemaster, N. 2008, ApJ, 688, 905

Parrish, I. J., McCourt, M., Quataert, E., & Sharma, P. 2012a, MNRAS, 422, 704

Parrish, I. J., McCourt, M., Quataert, E., & Sharma, P. 2012b, MNRAS, 419, L29

Perrone, L. M., & Latter, H. 2022a, MNRAS, 513, 4605

Perrone, L. M., & Latter, H. 2022b, MNRAS, 513, 4625

Perrone, L. M., Berlok, T., & Pfrommer, C. 2024a, A&A, 682, A125

Perrone, L. M., Berlok, T., & Pfrommer, C. 2024b, A&A, 690, A292

Peterson, J., & Fabian, A. 2006, Phys. Rep., 427, 1

Porter, D. H. 2000, ApJS, 127, 159

Prat, V., & Lignières, F. 2014, A&A, 566, A110

Pratt, G. W., Böhringer, H., Croston, J. H., et al. 2007, A&A, 461, 71

Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, Space Sci. Rev., 215, 25

Quataert, E. 2008, ApJ, 673, 758

Reboul-Salze, A., Guilet, J., Raynaud, R., & Bugli, M. 2022, A&A, 667, A94

Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, ApJ, 876, 85

Rieutord, M. 2015, Fluid Dynamics: An Introduction (Cham: Springer International Publishing)

Rincon, F. 2019, J. Plasma Phys., 85, 205850401

Riquelme, M. A., Quataert, E., & Verscharen, D. 2016, ApJ, 824, 123

Roulet, G., & Klein, P. 2009, J. Fluid Mech., 624, 45

Ruszkowski, M., & Oh, S. P. 2010, ApJ, 713, 1332

Ruszkowski, M., Lee, D., Brüggen, M., Parrish, I., & Oh, S. P. 2011, ApJ, 740, 81

Schekochihin, A. A., Cowley, S. C., Kulsrud, R. M., Hammett, G. W., & Sharma, P. 2005, ApJ, 629, 139

Schekochihin, A. A., Cowley, S. C., Rincon, F., & Rosin, M. S. 2010, MNRAS, 405, 291

Schwarzschild, K. 1906, Nachrichten von der Königlichen Gesellschaft der Wissenschaften zu Göttingen, 195, 41

Sharma, P., & Hammett, G. W. 2007, J. Comp. Phys., 227, 123

Simionescu, A., Allen, S. W., Mantz, A., et al. 2011, Science, 331, 1576

Spitzer, L. 1962, Physics of Fully Ionized Gases (New York: Interscience)

Tailleux, R. 2009, J. Fluid Mech., 638, 339

Tailleux, R. 2013, Phys. Scr., T155, 014033

Tailleux, R. 2018, J. Fluid Mech., 842, R1

Tailleux, R. 2024, J. Fluid Mech., 994, A5

Tailleux, R., & Dubos, T. 2024, Ocean Model., 188, 102339

Traxler, A., Garaud, P., & Stellmach, S. 2011, ApJ, 728, L29

Vaidya, B., Prasad, D., Mignone, A., Sharma, P., & Rickler, L. 2017, MNRAS, 472, 3147

Vazza, F., Angelinelli, M., Jones, T. W., et al. 2018, MNRAS, 481, L120

Vitense, E. 1953, Z. Astrophys., 32, 135

Volkov, A. N., Johnson, R. E., Tucker, O. J., & Erwin, J. T. 2011, ApJ, 729, L24

Winters, K. B., Lombard, P. N., Riley, J. J., & D'Asaro, E. A. 1995, J. Fluid Mech., 289, 115

Zhuravleva, I., Churazov, E. M., Schekochihin, A. A., et al. 2014, ApJ, 788, L13

Zhuravleva, I., Churazov, E., Arévalo, P., et al. 2015, MNRAS, 450, 4184

Zingale, M., Dursi, L. J., ZuHone, J., et al. 2002, ApJS, 143, 539

ZuHone, J. A., Markevitch, M., Ruszkowski, M., & Lee, D. 2013, ApJ, 762, 69
## Appendix A: Energy flux budget in B-MHD

The gravitational work on the right-hand-side of Eq. (9) is related to the gravitational energy $E_G = \rho\Phi$ according to

$$ \frac{\partial E_G}{\partial t} + \nabla \cdot (E_G \mathbf{v}) = \rho \frac{\partial \Phi}{\partial t} + \rho \mathbf{v} \cdot \nabla \Phi, \qquad (A.1) $$

where the first term on the right-hand-side is zero in our case, since the NFW gravitational potential of dark matter is independent of time in our models. A total budget equation of radial energy fluxes also including the gravitational energy can thus be obtained by spherically averaging the sum of Eqs. (A.1) and (9):

$$ \frac{\partial}{\partial t} \langle E + E_G \rangle_{\Omega} + \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 [Q_r + \mathcal{H}_r + \mathcal{G}_r + \mathcal{K}_r + \mathcal{D}_r + \Pi_r] \right) = 0, \qquad (A.2) $$

where

$$ Q_r = \langle -\kappa (\hat{\mathbf{b}} \cdot \nabla T) b_r \rangle_{\Omega}, \qquad (A.3) $$

$$ \mathcal{H}_r = \frac{\mathcal{R}\gamma}{\gamma - 1} \langle pv_r \rangle_{\Omega}, \qquad (A.4) $$

$$ \mathcal{G}_r = \Phi \langle \rho v_r \rangle_{\Omega}, \qquad (A.5) $$

$$ \mathcal{K}_r = \frac{1}{2} \langle \rho v^2 v_r \rangle_{\Omega}, \qquad (A.6) $$

$$ \mathcal{D}_r = \langle \mathbf{\Sigma} \cdot \mathbf{v} \rangle_{\Omega} \cdot \mathbf{e}_r, \qquad (A.7) $$

$$ \Pi_r = \frac{1}{\mu_0} \langle \mathcal{E} \times \mathbf{B} \rangle_{\Omega} \cdot \mathbf{e}_r, \qquad (A.8) $$

are the spherically averaged radial heat, enthalpy, gravitational, kinetic, viscous and Poynting fluxes, respectively. In Eq. (A.8), $\mathcal{E} = -\mathbf{v} \times \mathbf{B} + \eta \mathbf{j}$ is the local electric field of the plasma. In the ICM, the kinetic and magnetic energy volume densities are comparable, and several orders of magnitude smaller than the internal and gravitational energy densities. We therefore expect the kinetic, viscous and Poynting fluxes to contribute very little to the total energy flux budget; so we discard them from the current discussion. Focusing on the remaining heat, enthalpy, and gravitational fluxes, we note that they are themselves a sum of very different contributions, that can easily be highlighted when all physical quantities $X$ (other than the unit vector $\hat{\mathbf{b}}$) are further decomposed into their spherically averaged and fluctuating$^5$ parts $\langle X \rangle_{\Omega} + \delta X$. More specifically, the radial heat flux $Q_r$ can be rewritten as the sum of background $Q_{r,0}$ and perturbed $Q_{r,1}$ heat fluxes:

$$ Q_r = -\langle \kappa b_r^2 \rangle_{\Omega} \partial_r \langle T \rangle_{\Omega} - \langle \kappa (\hat{\mathbf{b}} \cdot \nabla \delta T) b_r \rangle_{\Omega}. \qquad (A.9) $$

The radial enthalpy flux $\mathcal{H}_r$ is the sum of an advective $\mathcal{H}_{r,0}$ and two convective $\mathcal{H}_{r,1} = \mathcal{H}_{r,\delta T} + \mathcal{H}_{r,\delta \rho}$ heat fluxes:

$$ \mathcal{H}_r = \frac{\mathcal{R}\gamma}{\gamma - 1} \left( \langle \rho T \rangle_{\Omega} \langle v_r \rangle_{\Omega} + \langle \rho \rangle_{\Omega} \langle \delta T \delta v_r \rangle_{\Omega} + \langle T \rangle_{\Omega} \langle \delta \rho \delta v_r \rangle_{\Omega} \right), \qquad (A.10) $$

and similarly for the gravitational flux, which is the sum of an advective $\mathcal{G}_{r,0}$ and a convective $\mathcal{G}_{r,1}$ flux:

$$ \mathcal{G}_r = \Phi \langle \rho \rangle_{\Omega} \langle v_r \rangle_{\Omega} + \langle \delta \rho \delta v_r \rangle_{\Omega}. \qquad (A.11) $$

## Appendix B: Introduction to the concept of available potential energy

In this work, we introduced the notion of available potential energy (APE) to understand the non-linear saturation of the MTI. The goal of this section is to better describe this physical quantity, since it might look, at first glance, unfamiliar in astrophysics.

Given the close relation, Eq. (36), between the specific internal energy (or enthalpy) and specific gravitational energy for an atmosphere at HSE in a gravitational field $\Phi$, Lorenz (1955) proposed to process both internal and gravitational energies as a single type of energy in a stratified atmosphere. The net gain of kinetic energy through adiabatic buoyancy motions comes at the expense of both energies. Therefore, Lorenz (1955) defined the sum of the gravitational potential and internal energies as the total potential energy (TPE) $E_P = E_G + E_I$. However, by definition of a stably stratified atmosphere at HSE (with respect to thermal convection), it seems that none of the gravitational energy $E_G$ of such an equilibrium state can be converted into kinetic energy, since small velocity perturbations will never grow but rather be systematically damped. This illustrative example shows that, in stratified fluid dynamics, not all the gravitational energy (and a fortiori not all the TPE) is available for conversion into kinetic energy. This motivated Lorenz (1955) to introduce the notion of available potential energy (APE) $E_A$, which quantifies the part of the TPE that is available for reversible conversions into kinetic energy. Formally, the APE is the difference between the TPE and the background potential energy (BPE) $E_B$, which is defined as the minimum TPE that can be obtained from an adiabatic (and thus reversible) redistribution of masses at HSE in the initial gravitational field $\Phi$ (Winters et al. 1995). With such a definition, the APE is nothing else than the amount of buoyant, and adiabatic compressive, work that must be exerted against the background stratification to bring a fluid element from its current physical condition (i.e. its height and pressure) to the position and pressure it would have in the reference state, that realises the BPE (that is the minimum TPE; Tailleux 2018). Guided by such considerations, Andrews (1981) first derived local evolution equations for both the APE and BPE, absent any non-ideal effect. Tailleux (2018) later generalised these equations for a viscous and thermally conducting fluid. Here, we extend it again, to include Ohmic dissipation:

$$ \frac{\partial E_A}{\partial t} + \nabla \cdot \left( E_A \mathbf{v} + \delta p \mathbf{v} + \frac{\delta T}{T} \mathbf{q} \right) = \delta \rho \mathbf{v} \cdot \nabla \Phi + \mathbf{v} \cdot \nabla \delta p + \mathbf{q} \cdot \nabla \left( \frac{\delta T}{T} \right) + \frac{\delta T}{T} \left( \mathbf{\Sigma} : \nabla \mathbf{v} + \eta \mu_0 j^2 \right) + \mathcal{F}, \qquad (B.1) $$

$$ \frac{\partial E_B}{\partial t} + \nabla \cdot \left( E_B \mathbf{v} + p_* \mathbf{v} + \frac{T_*}{T} \mathbf{q} \right) = - \mathbf{q} \cdot \nabla \left( \frac{\delta T}{T} \right) + \frac{T_*}{T} \left( \mathbf{\Sigma} : \nabla \mathbf{v} + \eta \mu_0 j^2 \right) - \mathcal{F}, \qquad (B.2) $$

with $\mathbf{q}$ as some heat flux and where the fluctuating quantities here are defined with respect to the reference state $p_*$ that realises the BPE. In Sect. 3.4 however, these fluctuations are defined locally, with respect to the spherically averaged quantities $\langle \cdot \rangle_{\Omega}$. Therefore, we implicitly chose the background profiles as the reference state, which is, however, not guaranteed to be that of minimum TPE in this case. What exact information (and if any, which amount) is lost when choosing an arbitrary reference state, rather than the sorted state, in such studies is still to be elucidated. In the same manner, it must be noted that the quadratic APE introduced in Eq. (16) is only an approximation valid

<footnote>$^5$ By definition, $\langle \delta X \rangle_{\Omega}$ is zero.</footnote>
in the limit of small displacements and density perturbations (Roullet & Klein 2009; Tailleux 2013), and does certainly not hold for arbitrarily large displacements (Holliday & McIntyre 1981; Andrews 1981), despite having been used as such in simulations of compressible stratified ICM turbulence here or in Mohapatra et al. (2020). The APE and the BPE are both extrinsic state functions, that is a joint function depending on the state variables describing the local fluid properties, as well as the properties of the environment within which it is embedded. Therefore, an additional exchange term $F$ between APE and BPE accounting for the temporal variation of the environment itself, that is the fluid as a whole, appears in both Eqs. (B.1) and (B.2) (Tailleux 2018). Though we acknowledge that a complete volume-averaged budget of APE would require taking into account the diverse boundary fluxes and the non-local exchange term too, we just used the latter equation to illustrate the basic phenomenology of the MTI energetics in Sect. 4.1.

From Eqs. (B.1) and (B.2), it is readily shown that the sum of the evolution equations for the available $E_A$ and background potential energy $E_B$ gives the right evolution equation for the total potential energy $E_P$, under the assumption that the reference state is at HSE ($\rho_* \nabla\Phi + \nabla p_* = 0$):

$$ \frac{\partial E_P}{\partial t} + \nabla \cdot (E_P \mathbf{v} + p \mathbf{v} + \mathbf{q}) = \rho \mathbf{v} \cdot \nabla \Phi + \mathbf{v} \cdot \nabla p + \boldsymbol{\Sigma} : \nabla \mathbf{v} + \eta \mu_0 j^2. \quad (B.3) $$

The conversion rate $\varepsilon_A$, defined as

$$ \varepsilon_A = \mathbf{q} \cdot \nabla \left( \frac{\partial T}{T} \right), \quad (B.4) $$

is absent of the latter evolution equation. This quantity is then a conversion rate between APE and BPE, due to heat flux. On the contrary, the last two terms on the right-hand side of Eq. (B.1) are only indirectly related to such a conversion process: they occur in the conversion rate of the sum of the kinetic energy and the APE into BPE. In the ocean and in the atmosphere, the conversion rate between APE and BPE $\varepsilon_A$ is usually referred to as a dissipation rate (going from APE into BPE) because, most of the time, it is negative in physical conditions representative of these systems, and even negative-definite in case of an incompressible Boussinesq stratified single-component fluid (Winters et al. 1995). However, Middleton & Taylor (2020), Middleton et al. (2021), Tailleux (2024) recently showed that this rate can turn positive when double-diffusive instabilities are developing. In this work, we clearly demonstrated that this conversion rate can also be positive (from BPE to APE) in the case of the MTI, and that it even controls its saturation. This feature should not be surprising given the obvious similarities that the MTI shares with other double-diffusive instabilities. For example, their maximum growth rates do not depend on the level of diffusivity itself, and they all modify the Schwarzschild's criterion for stability, on account on their diffusive nature (although different diffusive processes are at stake in each case). We therefore speculate that the energetics of double-diffusive instabilities (Garaud 2018), such as the fingering instability for instance (Traxler et al. 2011; Fraser et al. 2024), may be similarly controlled by the positive injection rate of BPE into APE. Finally, we stress that the possible indirect conversion of BPE into KE (through APE) is only made possible here thanks to non-ideal effects; while Lorenz (1955) first defined the APE so as to distinguish the part of the TPE that is available for reversible conversion into kinetic energy from the part that is not (the BPE). In thermodynamics, available energy is sometimes referred to as exergy, defined as the maximum useful work that can be produced by a physical system as it is brought to equilibrium with its environment through

an ideal process (Cengel et al. 2011). When the physical system under consideration is a perfect stratified fluid, the APE seems to play exactly this role (Kucharski 1997). As highlighted throughout this paper however, the relaxation of the MTI is very far from being adiabatic. Even more surprising, in this case, kinetic energy (i.e. useful work) can only be produced thanks to non-ideal effects such as parallel thermal conduction. Therefore, new ways of defining the exergy budget of a thermally conducting dilute stratified fluid, such as the ICM, should be explored to circumvent this theoretical limitation. For instance, partitioning the TPE through an isothermal, rather than adiabatic, rearrangement of masses might represent a promising avenue in this direction.

## Appendix C: Influence of the boundary conditions applied on the stratified thermodynamic fields

In all the simulations presented in Table 2, we set indirect inhomogeneous Dirichlet BC on the temperature field by imposing, all along the runs, both the density and the pressure in the ghost cells to be those of the initial ICM atmosphere at HSE. The main limitation of such BCs is the possible development of small boundary layers on the density and temperature fields, as can be seen at the inner and outer boundaries on the background density profile, $\langle\rho\rangle_\Omega$, in the 2D run S0 (green curve on the leftmost panel in Fig. 1), or on the background heat flux, $Q_{r,0}$, in the 3D run F0 (yellow curve on the top panel in Fig. 6).

As an attempt to dispose of such boundary layers, we implemented a first-order local HSE BC, as suggested by Zingale et al. (2002, Eq. 41):

$$ p_{i+1} - p_i = - (\Phi_{i+1} - \Phi_i) \frac{\rho_{i+1} + \rho_i}{2}, \quad (C.1) $$

which must be solved for the two unknowns $p_{i+1}$ and $\rho_{i+1}$ in the case of the outer (rightmost) ghost cells, or for both $p_i$ and $\rho_i$ at the inner boundary. Imposing inhomogeneous Dirichlet BC on the temperature field in the ghost cells is then enough to determine both quantities, thanks to the equation of state. With these new BCs on the stratified thermodynamic fields, coupled to homogeneous Dirichlet BC on the radial velocity, we reran the 2D simulation 2dDf1e-2 for 200 dynamical times $t_{dyn}$. In Fig. C.1, we plot the background density, temperature, and entropy profiles at various times during the simulation; while in Fig. C.2, we display the time evolution of the diverse energetic diagnostics introduced in Sect. 3.4. Both figures should be compared to their counterparts, Figs. 1 and 2, respectively. Several features of this new run are noteworthy. First, these new BCs do not prevent the formation of boundary layers on the density, temperature, and entropy profiles. Second, the thermodynamic profiles seem to keep evolving all along the simulations and to never saturate, or at least to do so on a timescale much larger than the usual simulation, such as the 2dDf1e-2 or F0 runs. As a result of the constantly evolving temperature profile, MTI-driven turbulence never reaches a proper quasi-steady state during the simulation. As seen on the middle panel in Fig. C.1, the average temperature through the domain is globally increasing with time. Indeed, a local HSE BC tends to suppress the thermal wind, whose intensity is then much lower in this run than in its counterpart 2dDf1e-2. Consequently, the thermal wind is no longer able to advect enough internal energy out of the domain: the fluid is constantly heated, without being sufficiently cooled down, since the conductive heat flux is larger at the inner than at the outer boundary. This result suggests that the thermal wind acts a self-regulation thermal mechanism, in a fluid heated from
<figure><img src="image_8.png" /><figcaption>**Fig. C.1.** Time evolution of the stratified ICM atmosphere during the simulation with the local HSE BC of Zingale et al. (2002). Shown are the density $\rho_0$ (in g/cm$^3$, left); temperature $T_0$ (in keV, centre); and entropy $S_0$ (right) profiles. In dashed lines: Initial profiles. In plain lines: Same profiles but displayed at various times. This figure should be compared to Fig. 1, which was obtained with the BCs described in Sect. 3.3.</figcaption></figure>

below. Despite the highlighted boundary layers and the unsteady
saturation, the regime of MTI-induced turbulence, and its sub-
sequent phenomenology, remain qualitatively unchanged in vol-
ume, with very similar zeroth-order levels of saturation than in
the initial run 2dDf1e-2, as seen from the comparison between
the opaque curves (current simulation), and their transparent
counterparts (2dDf1e-2 run), in Fig. C.2. Therefore, our con-
clusions on the nature of magneto-thermal turbulence are quite
independent on the BCs used in this work.

<figure><img src="image_9.png" /><figcaption>**Fig. C.2.** Comparison of the saturation regimes in the simulation with the local HSE BC of Zingale et al. (2002) (opaque lines) and in the 2dDf1e-2 run with the BCs described in Sect. 3.3 (transparent lines). Top: Volume-averaged time evolution of the available potential energy density $\langle E_A \rangle_\nu$, and of the different contributions to the kinetic $\langle E_K \rangle_\nu$ volume energy density. Bottom: Volume-averaged time evolution of the APE injection rate $\langle \epsilon_i \rangle_\nu$, APE thermal dissipation rate $-\langle \epsilon_k \rangle_\nu$, and the (opposite amount of) reversible buoyancy work $-\langle \mathcal{A} \rangle_\nu$ during the same run. Dashed lines are indicative of negative quantities whose sign has been switched for the purpose of visualisation.</figcaption></figure>