© The Authors 2025

1. Introduction

Galaxy clusters are the largest gravitationally bound objects in the universe, with a typical length scale being the virial radius, R_vir ∼ Mpc, and with a total mass ranging from 10¹⁴ to 10¹⁵ solar masses (Girardi et al. 1998). The gravitational potential well maintaining the internal coherence of such massive structures originates from yet-to-be-detected dark matter (Zwicky 1933; Smith 1936), which is responsible for ∼85% of their total mass (see for example McNamara & Nulsen 2007, and references therein). Galaxies, observable in the optical wavelengths, only represent ∼3% of this mass. The remaining ∼12% mass is found as a tenuous and very hot plasma, usually referred to as the intracluster medium (ICM), that shines in X-rays (Fabian et al. 2003, 2006). Typical values of density and temperature in the ICM are ρ ∼ 10⁻²⁷ g/cm³ and k_BT ∼ 5 keV, respectively, where k_B is the Boltzmann constant. Radio observations of galaxy clusters also reveal that these astrophysical objects are magnetised up to B ∼ 10 μG (Ferrari et al. 2008).

From a plasma physics standpoint, such measurements allow for the computation of typical plasma parameters (Huba 1998). For instance, the Larmor radius of charged particles in the ICM reaches no more than the nanoparsec scale, while their collisional mean free path can reach the kiloparsec scale (Kunz et al. 2022). On the macroscopic level, the subsequent ‘dilute’ regime of the ICM plasma induces an anisotropic transport of heat, along magnetic field lines only (Braginskii 1965),

$$\begin{array}{c}\hfill {\mathit{q}}_{\mathbf{B}}=-\kappa \widehat{\mathit{b}}\widehat{\mathit{b}}·\mathbf{\nabla }T,\end{array}$$(1)

usually referred to as the Braginskii heat flux. In this equation, $\widehat{\mathit{b}}=\mathit{B}/B$ is the unit magnetic field vector andκ is the thermal conductivity, which should not be confused with the thermal diffusivity, χ = κ/(ρℛ), ℛ being the specific gas constant. The anisotropic heat flux triggers a new class of buoyancy instabilities, possibly relevant to the ICM: the magneto-thermal instability (MTI) where the background temperature is decreasing with radius (Balbus 2000) and the heat-flux buoyancy driven instability (HBI) in case of an increasing background temperature (Quataert 2008). As seen with X-ray observations, the ICM systematically features negative temperature gradients in outer cluster regions, from a typical cooling radius for relaxed clusters (Leccardi & Molendi 2008; Simionescu et al. 2011; Ghirardini et al. 2019), which are thus unstable to the MTI in principle. Perrone & Latter (2022a,b), Kempf & Rincon (2025, respectively, PL22a, PL22b, KR25, hereafter) showed that the non-linear saturation of the MTI can drive fluid turbulence at levels roughly consistent with those deduced from X-ray observations of galaxy clusters (Hitomi Collaboration 2016, 2018; XRISM Collaboration 2025a,b). The characterisation of such turbulence is critical for our understanding of the ICM physics: it could induce non-thermal pressure support (Pratt et al. 2019), turbulent heating (Zhuravleva et al. 2014), and magnetic field amplification (Donnert et al. 2018; Rincon 2019), to name but a few of the most outstanding problems related to galaxy clusters.

Both the MTI and the HBI are local instabilities that share a common maximum growth rate, M₀ ∼ Gyr⁻¹. The local dispersion relation for MTI elevator modes (i.e. those with zero vertical wavenumber) was originally found by Balbus (2000):

$$\begin{array}{c}\hfill D(\sigma ,k_x)={\sigma }^3+\stackrel{\sim }{\chi }{k}_x^2{\sigma }^2+N_0^2\sigma -\stackrel{\sim }{\chi }{k}_x^2{M}_0^2=0,\end{array}$$(2)

where σ and k_x are respectively the growth rate and the horizontal wavenumber of the local MTI mode, while $\stackrel{\sim }{\chi }$ and N₀ are respectively the reduced thermal diffusivity and the Brunt-Väisälä frequency, defined later in Sect. 2. Equation (2) is only valid for local elevator modes, whose typical length scales are much smaller than any ICM scale height. However, the MTI is known to preferentially develop at the conduction length, ${\ell }_{\chi }=\sqrt{\chi /M_0}$ (PL22a; PL22b), which can be a non-negligible fraction of the cluster size itself (Kempf et al. 2023; KR25). At these scales, appreciable variation of the background ICM thermodynamic profiles necessarily happens. Consequently, the validity of local models may be called into question when describing the linear stability and nonlinear saturation of the MTI. While KR25 previously addressed the question of the non-linear saturation thanks to global simulations of magneto-thermal turbulence, a global linear analysis of the MTI remains to be carried out (for a first approach to this problem though, see Appendix B in Berlok et al. 2021). Such a study would help to further understand why the phenomenology developed by PL22a; PL22b for MTI-driven turbulence in local ICM atmosphere successfully applies to global simulations too (KR25).

In the case of the HBI, previous global linear analyses unveiled interesting vertical structures of the corresponding eigenmodes (Latter & Kunz 2012; Berlok & Pessah 2016); but they did not establish any formal link between global modes and their local counterparts. On the other hand, Latter et al. (2015, LFF15 hereafter) carried out a global linear analysis of the magneto-rotational instability (MRI) in accretions discs, and elucidated the relationship between local and global modes of the MRI. As first intuited by Balbus (2001), the local stability properties of the MTI and the MRI are very much alike. Both instabilities modify, in the same way, the classical stability criteria of the gas. Therefore, we build a global linear theory of the MTI, drawing from that of LFF15 for the MRI, exploiting their common linear features. The main goal of this paper is to shed light on the connection between local and global MTI modes, and between the global linear MRI and MTI too.

Eventually, the non-linear study of such instabilities at saturation relies on large, complex, compressible, MHD codes for the simulation of astrophysical flows. However, only few detailed analytical benchmarks are found in the literature for global curvilinear geometries (LFF15), especially so for numerical operators as exotic as the anisotropic Braginskii heat flux (Latter & Kunz 2012; Berlok & Pessah 2016). A global linear analysis of the MTI in the ICM would therefore provide an additional test, in spherical geometry, to validate the computational fluid dynamics tools on which a lot of theoretical astrophysics relies nowadays.

The outline of the paper is as follows. In Sect. 2, we first describe the ICM model, along with the analytical and numerical methods, used to carry out a global linear analysis of the MTI. In Sect. 3, we present several new interesting analytical results (some of which are deferred to Appendix A). We derive the radial structures of global MTI eigenmodes and unveil their deeper relationships with local MTI modes, and with global MRI eigenmodes in accretion discs. In Sect. 4, these analytical results are further corroborated by numerical solutions obtained with both a Chebyshev pseudo-spectral solver and direct numerical simulations performed with the finite-volume code IDEFIX. Our main conclusions are summarised in Sect. 5.

2. Model and methods

In this section, we first introduce the fully compressible Braginskii-MHD (B-MHD) equations. We then present the equilibrium model of a stratified ICM, around which we carried out a global linear analysis of the MTI. Finally, we describe the numerical methods used both to solve the subsequent eigenproblem and to run direct numerical simulations.

2.1. Model of the intracluster medium

2.1.1. The compressible Braginskii magnetohydrodynamic equations

The B-MHD framework provides evolution equations for the density, momentum, and (internal) energy of a dilute, collisional plasma subject to anisotropic transport or heat and/or momentum. Our B-MHD equations read

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\mathbf{\nabla }·(\rho \mathit{v})& =0,\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \frac{\partial (\rho \mathit{v})}{\partial t}+\mathbf{\nabla }·(\rho \mathit{v}\mathit{v}+p\mathit{I}+\frac{B^2}{2{\mu }_0}(\mathit{I}-2\widehat{\mathit{b}}\widehat{\mathit{b}}))& =\rho \mathit{g},\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill \frac{\mathcal{R}}{\gamma -1}[\frac{\partial \left(\rho T\right)}{\partial t}+\mathbf{\nabla }·\left(\rho T\mathit{v}\right)]& =-\mathbf{\nabla }·{\mathit{q}}_{\mathbf{B}}-p\mathbf{\nabla }·\mathit{v},\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{B}}{\partial t}& =\mathbf{\nabla }\times (\mathit{v}\times \mathit{B}).\hfill \end{array}$$(6)

All the symbols here have their usual meanings, which are made explicit in KR25. The gravitational acceleration, denoted g, is specified in the next paragraph. The heat flux, q_B, in a dilute, collisional plasma is given by Eq. (1). In this analysis, we dismissed any non-ideal effects (such as isotropic viscosity or magnetic resistivity) except for the anisotropic transport of heat. Specifically, we neglected the Braginskii viscosity, which otherwise accounts for the effects of anisotropic pressure (Snyder et al. 1997). Besides, we ignored any plasma effects (such as kinetic micro-instabilities) that could affect the anisotropic heat flux on which the MTI relies (Berlok et al. 2021; Perrone et al. 2024a,b). The system of Eqs. (3), (4), (5), and (6) was closed thanks to the equation of state for a perfect gas with adiabatic index γ.

2.1.2. Initial ICM atmosphere at hydrostatic equilibrium

Since our goal was to carry out a global linear analysis of the MTI in the spherical equatorial plane, θ = π/2, of a stratified ICM, we therefore introduced a 2D spherical coordinate system, $({\widehat{\mathit{e}}}_r,{\widehat{\mathit{e}}}_{\phi })$. The position vector is $\mathit{r}=r{\widehat{\mathit{e}}}_r+\phi {\widehat{\mathit{e}}}_{\phi }$. The domain extends from an inner radius, r₁, to an outer radius, r₂, and from 0 to 2π in azimuth, φ. We took advantage of the global spherical geometry at hand to set up a background ICM atmosphere as realistic as possible. In galaxy clusters, the gravitational potential well of dark matter is well described by the Navarro-Frenk-White potential (NFW; Navarro et al. 1997). We used a NFW gravitational acceleration field, $\mathit{g}=-g_0{\widehat{\mathit{e}}}_r$, of intensity

$$\begin{array}{c}\hfill g_0(r)=\left(\frac{v_{\mathrm{vir}}^2}{R_{\mathrm{vir}}}\right)\frac{1}{r^2}\frac{log(1+cr)-cr/(1+cr)}{log(1+c)-c/(1+c)},\end{array}$$(7)

where r is naturally rescaled to the virial radius, R_vir, and c is a concentration parameter that we set to c = 5. We also introduced the virial velocity, v_vir = (GM_vir/R_vir)^1/2, where M_vir is the virial mass (not to be confused with the MTI frequency, M₀, in Eq. (2)).

We considered a simple, purely azimuthal, initial magnetic field, $\mathit{B}=B_0{\widehat{\mathit{e}}}_{\phi }$, with B₀ uniformly weak and independent of the radius, r, across the domain. This magnetic configuration is the most unstable to the MTI. Our choice was driven by pragmatism and tractability. In a 2D equatorial model of an ICM atmosphere, purely azimuthal fields are the only possible magnetic configurations unstable to the MTI, also ensuring the thermal equilibrium of any spherically symmetric temperature profile, since the background radial heat flux is initially shut off.

The temperature profile, T₀(r), is then a free parameter. Given the theoretical nature of the work undertaken here, we favoured a simple analytic function, rather than a complex fit to the observational data (Ghirardini et al. 2019; Kempf et al. 2023). We adopted a temperature profile at thermal equilibrium in the presence of a hypothetical background heat flux (with isotropic conductivity, or equivalently, anisotropic conductivity and purely radial magnetic field), a reasonable compromise between physical plausibility and analytical tractability. A thermal equilibrium in the presence of isotropic heat flux and uniform thermal conductivity is obtained by setting

$$\begin{array}{c}\hfill T_0(r)=T_1+(T_2-T_1)\frac{r_1-r}{r_1-r_2}\frac{r_2}{r},\end{array}$$(8)

subject to boundary conditions (BCs) T(r₁) = T₁ and T(r₂) = T₂ < T₁, at the inner and outer radii, r₁ and r₂, respectively. We set the values of these free parameters, based on observations of the Perseus cluster (Simionescu et al. 2011), to k_BT₁ = 6 keV and k_BT₂ = 2 keV, for r₁ = 0.3R_vir and r₂ = R_vir, respectively. For Perseus, we have R_vir = 1.1 Mpc. Other virial quantities are given in KR25. As shown in Fig. 1, with such choices, the resulting profile of MTI frequency is a well-behaved, monotonically decreasing function of the radial coordinate,

$$\begin{array}{c}\hfill M_0(r)=\sqrt{-g_0\frac{\mathrm{d}log{T}_0}{\mathrm{d}r}}·\end{array}$$(9)

Fig. 1. Local MTI frequency profile, M0, as a function of the radius, r.

The MTI frequency at the innermost radius, r₁, is M₁ = 2.44 Gyr⁻¹. Since the profile is decreasing with radius, it is also the maximum MTI frequency within the domain. For further reference, we also define M₂, the (smallest) MTI frequency at the outer boundary, r₂.

Finally, a thermal diffusivity profile, χ(r), needs to be specified too. Strictly speaking, this physical quantity should be the Spitzer diffusivity, χ_S ∝ T^5/2/ρ, in the ICM plasma (Spitzer 1962). However, we show in Sect. 3 that the precise dependency of the thermal diffusivity only plays a modest role, at least in the limit of fast thermal diffusion that we introduce later. Therefore, so as not to complicate the interpretation of both the analytical derivations and numerical results presented in the next sections, we used a uniform profile of thermal diffusivity, set to the Spitzer value it would have at the inner boundary,

$$\begin{array}{c}\hfill {\chi }_{\mathrm{s},1}=0.11R_{\mathrm{vir}}{v}_{\mathrm{vir}}{\left(\frac{k_{\mathrm{B}}{T}_1}{6\mathrm{keV}}\right)}^{\frac{5}{2}}{\left(\frac{n_{\mathrm{e},1}}{{10}^{-3}{\mathrm{cm}}^{-3}}\right)}^{-1},\end{array}$$(10)

where n_e, 1 is the electronic density number, close to 10⁻³ cm⁻³ at r₁ = 0.3R_vir in the Perseus cluster (Simionescu et al. 2011). As a result, the density profile, ρ₀(r), needs not be explicitly computed to solve the perturbation equations, introduced in the next subsection.

2.2. Global linear analysis

In order to carry out a global linear analysis of the MTI in the ICM, the compressible B-MHD Eqs. (3), (4), (5), and (6) were linearised to first order around the previous background HSE. All perturbed quantities were azimuthally expanded into monochromatic Fourier modes since this direction is naturally periodic in spherical geometry and the background state is spherically symmetric. As a starting point for this study, we discarded any dependency of global MTI modes on the polar angle, θ, and only considered radial and azimuthal variations. Their radial structure was then captured into a function of the radius only: δf(t, r, φ) = δf(r)exp(σt+imφ). This ansatz is appropriate to describe a 2D equatorial cut of MTI modes with large scale variation in the polar direction, but smaller scale (though still global) variations in the other directions. Moreover, guided by the local theory, we knew that the main features of the MTI derive in a 2D framework, so that the θ-variation can be safely neglected. This also leaves the problem analytically tractable. Accordingly, the amplitude of our initial magnetic field, B₀, is independent of the polar coordinate, θ, in the 2D equatorial cut. As discussed in Sect. 5, the cost for moving beyond this approximation would be to solve a 2D, rather than a 1D, eigenproblem because the Braginskii heat flux is a non-linear operator. A very similar approximation was recently used by Choudhury & Reynolds (2025), to linearly study cold fronts in a 2D equatorial ICM plane.

2.2.1. Perturbation equations

From the linearisation at first order of Eqs. (3), (4), (5), and (6), the perturbation equations are

$$\begin{array}{cc}\hfill \sigma \frac{\delta \rho }{{\rho }_0}& =-\frac{\mathrm{i}m}{r}\delta v_{\phi }-\frac{1}{r^2{\rho }_0}\frac{\mathrm{d}}{\mathrm{d}r}\left(r^2{\rho }_0\delta v_r\right),\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill \sigma \delta v_{\phi }& =-\mathcal{R}{T}_0\frac{\mathrm{i}m}{r}\frac{\delta p}{p_0},\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill \sigma \delta v_r& =-\mathcal{R}{T}_0\frac{\mathrm{d}}{\mathrm{d}r}\left(\frac{\delta p}{p_0}\right)+g_0\frac{\delta T}{T_0},\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill \sigma \frac{\delta T}{T_0}& =-(\gamma -1)[\frac{\mathrm{i}m}{r}\delta v_{\phi }+\frac{1}{r^2}\frac{\mathrm{d}}{\mathrm{d}r}\left(r^2\delta v_r\right)]-\frac{\mathrm{d}log{T}_0}{\mathrm{d}r}\delta v_r\hfill \end{array}$$(14)

$$\begin{array}{cc}& -\gamma \stackrel{\sim }{\chi }\frac{m^2}{r^2}(\frac{\delta T}{T_0}-\frac{\mathrm{d}log{T}_0}{\mathrm{d}r}\frac{\delta A}{B_0}),\hfill \\ \hfill \sigma \frac{\delta A}{B_0}& =-\delta v_r,\hfill \end{array}$$(15)

where we introduced the reduced thermal diffusivity, $\stackrel{\sim }{\chi }(r)=\chi (r)(\gamma -1)/\gamma$. The pressure perturbation is simply deduced from the perturbed equation of state,

$$\begin{array}{c}\hfill \frac{\delta p}{p_0}=\frac{\delta \rho }{{\rho }_0}+\frac{\delta T}{T_0}·\end{array}$$(16)

The perturbed potential vector, $\mathit{\delta }\mathit{A}=\delta A{\widehat{\mathit{e}}}_{\theta }$, relates to the perturbed magnetic field through δB = ∇ × δA. In the linearised Eqs. (12) and (13) for the velocity perturbations, we ignored the feedback from the magnetic field on the flow through the action of the Lorentz force. This simplification is justified when the plasma beta parameter is very large, that is in the limit β ≫ 1. This assumption holds to some extent for the plasma of the ICM. In this context, the only dynamical role of the magnetic field is to channel heat through the Braginskii heat flux, Eq. (1). The set of linearised equations above is solved numerically in Sect. 4, but is slightly simplified, in the next subsection, for the need of the analytical derivations carried out in Sect. 3.

2.2.2. Sound-proof approximation

Before the global linearised Eqs. (11), (12), (13), (14), and (15) can be reworked into a simple and useful ordinary differential equation (ODE), we needed to simplify Eq. (11). We adopted a type of anelastic model, in which we assumed that the flow is very subsonic and that the fractional thermodynamic perturbations remain small. Doing so allowed us to drop the time derivative of the density perturbation, σδρ/ρ₀, in Eq. (11),

$$\begin{array}{c}\hfill 0=\frac{\mathrm{i}m}{r}\delta v_{\phi }+\frac{1}{r^2{\rho }_0}\frac{\mathrm{d}}{\mathrm{d}r}\left(r^2{\rho }_0\delta v_r\right).\end{array}$$(17)

This approximation filters out sound waves that would otherwise interfere with the computation of the linear MTI modes, in which we were interested. This approach is justified as long as the Mach number of the flow is small enough¹, that is $\mathcal{M}=\delta v/c_{\mathrm{s}}\ll 1$, with the sound speed, $c_{\mathrm{s}}=\sqrt{\gamma \mathcal{R}{T}_0}$. This is guaranteed for arbitrarily small velocity perturbations. We stress that this procedure is not strictly equivalent to an anelastic approximation, at least not as usually found in solar physics (Miesch 2005), even though they look very similar. The main difference here is that we do not require the background temperature gradient to be the adiabatic gradient (Rieutord 2015), to within a small parameter. Indeed, the ICM is stable to thermal convection, precisely because this background gradient is usually much lower than the adiabatic gradient in galaxy clusters (Cavagnolo et al. 2009).

2.3. Numerical methods

We first describe the Chebyshev pseudo-spectral method designed to solve the full eigenproblem Eqs. (11), (12), (13), (14), and (15). We then present the Godunov-type code IDEFIX for astrophysical flows, used in Sect. 4 to cross-check the radial structure and growth rate of the analytically derived global MTI eigenmodes.

2.3.1. Chebyshev pseudo-spectral solver

We discretised Eqs. (11), (12), (13), (14), and (15) using a pseudo-spectral Chebyshev method (Boyd 2001) on a Gauss-Chebyshev grid. Only two BCs need to be specified because, as shown in Sect. 3, these equations can be rewritten as a second-order ODE. Dirichlet BCs were then imposed on the radial velocity field, through basis recombination of the Chebyshev polynomials. The eigenvalues, σ_n, of the resulting matrix were computed with SCIPY (Virtanen et al. 2020). The eigenvalue problem was then solved at fixed azimuthal order, m, for different n-indexed modes which are decreasingly ranked according to their growth rates. Our solver was adapted from previous quasi-global linear analyses of either the HBI or the MTI (Latter & Kunz 2012; Berlok et al. 2021), and was successfully tested against those. We did not use a single grid resolution for all computations. Instead, we required that the variation of the numerically computed eigenvalue of the eigenmode under consideration did not exceed a tolerance threshold of 10⁻⁵ when the grid resolution was doubled. The number of grid points needed for the numerical convergence of most modes presented here usually did not exceed 512. However, it was increased up to 2048 when we needed to solve for very large mode numbers, up to n ≈ 1000, in Sect. 4.1.

2.3.2. Direct numerical simulations with IDEFIX

IDEFIX is a new Godunov-type finite-volume code for astrophysical fluid dynamics (Lesur et al. 2023). It incorporates a Braginskii module, which includes in particular the anisotropic heat flux, Eq. (1). In this work, we used IDEFIX to solve the fully compressible B-MHD Eqs. (3), (4), (5), and (6). We note however that IDEFIX solves the conservation equation for the total energy (KR25, Eq. (9)) instead of the conservation equation for the internal energy, Eq. (5). The numerical methods used in all the simulations presented in Sect. 4 are the following: piece-wise linear reconstruction with a Van Leer slope limiter, HLLD Riemann solver, three-stage Runge-Kutta time integrator, super time-stepping Runge-Kutta-Legendre scheme (Meyer et al. 2014; Vaidya et al. 2017) for the Braginskii heat flux (no slope limiter though, Sharma & Hammett 2007), and upwind constrained transport contact algorithm (Londrillo & Del Zanna 2004). Even though the simulations were ran in spherical geometry, we only solved the equations in a 2D domain here, that is the equatorial plane, $({\widehat{\mathit{e}}}_r,{\widehat{\mathit{e}}}_{\phi })$, at θ = π/2, extending from r₁ = 0.3 to r₂ = 0.7 in radius and from 0 to 2π (or from 0 to π) in the azimuthal direction. When a smaller azimuthal extent is used, from φ = 0 to π, periodic BCs remain applicable provided that the azimuthal order, m, of the initial condition is even and greater than or equal to 2 because e^imφ is 2π/m-periodic, and, thus, at least π-periodic. Accordingly, the domain extent used in IDEFIX was always smaller in radii (and sometimes in azimuth too) than that used in the pseudo-spectral solver. The reason behind this choice was to ensure that IDEFIX solved the global MTI eigenmodes with enough radial (respectively azimuthal) resolution, when the azimuthal order, m, (resp. mode number, n), of the modes became large. The numerical resolution used was therefore 2048 × 4096. The initial parameters were χ = 0.11 (in line with Eq. (10)) and B₀ = 10⁻⁶, that is β ∼ 10¹².

To obtain a true static equilibrium upon which the MTI can grow, along the lines of Eqs. (11), (12), (13), (14), and (15), we chose an ICM atmosphere initially at HSE. However, imposing a HSE at machine precision in any Godunov-type codes, such as IDEFIX, is difficult (Zingale et al. 2002), because HSE is numerically obtained via the (only approximate) cancellation of two possibly large numbers: the pressure gradient and the gravitational pull. If not done properly, or without sufficient radial resolution, a spherically symmetric, m = 0, acoustic wave can develop. If the amplitude of this acoustic perturbation exceeds the initial amplitude of the seeded MTI eigenmodes, the latter can be disrupted, and unwanted modes with other azimuthal orders might appear. To avoid such pathological behaviours, we made use of a balance scheme, available in IDEFIX. The goal of this module is to enforce the initial condition, which needs to be an analytic equilibrium, as a true numerical equilibrium. It does so by initially computing, and then removing at each time step, the right-hand side numerical residual of the initial condition, so that the analytical equilibrium is recovered at machine precision. All simulations performed in Sect. 4 employed this balance scheme. The discussion about the BCs used in these simulations is deferred there, since they depend on the specific radial structure of the MTI eigenmodes, which we obtain later.

3. Analytical results

In this section, we present our analytical results on the global linear MTI eigenproblem. First, we elucidated the single ODE governing the MTI eigenmodes, and we applied it a specific limit of fast thermal diffusion to obtain a problem that is to some extent analytically tractable. Then, we qualitatively discuss the physical phenomenology of the global MTI modes induced by the simplified Schrödinger potential of the governing ODE. More quantitative results were obtained too, thanks to known theorems from the Sturm-Liouville theory; but they are deferred to Appendix A to avoid breaking the flow of the current section. Finally, we obtained convenient analytical approximations to the growth rates of the MTI eigenmodes, using WKBJ approximations in the limits of large azimuthal orders, m/n ≫ 1, or large mode numbers, n/m ≫ 1. The detailed derivations leading to these approximate growth rates are also deferred to Appendix B.

3.1. Single governing ordinary differential equation

First, we transformed Eq. (14), using Eqs. (11) and (16), so as to introduce the Brunt-Väisälä frequency, defined as

$$\begin{array}{c}\hfill N_0^2(r)=\frac{g_0}{\gamma }\frac{\mathrm{d}log\left(p_0{\rho }_0^{-\gamma }\right)}{\mathrm{d}r},\end{array}$$(18)

in:

$$\begin{array}{cc}\hfill \sigma (\frac{\delta T}{T_0}-\frac{\gamma -1}{\gamma }\frac{\delta p}{p_0})& =-\frac{N_0^2}{g_0}\delta v_r\hfill \end{array}$$(19)

$$\begin{array}{cc}& -\stackrel{\sim }{\chi }\frac{m^2}{r^2}(\frac{\delta T}{T_0}-\frac{\mathrm{d}log{T}_0}{\mathrm{d}r}\frac{\delta A}{B_0}).\hfill \end{array}$$(20)

The latter equation is exact, with respect to Eq. (14). No further approximation was used to derive it. The sound-proof Eq. (17), combined with Eqs. (12), (13), (15), and (20), can then be reworked into a single second-order ODE, for the variable r² ρ₀δv_r,

$$\begin{array}{cc}\hfill [\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}& +\frac{1}{H}\frac{\sigma /\gamma +\stackrel{\sim }{\chi }{m}^2/r^2}{\sigma +\stackrel{\sim }{\chi }{m}^2/r^2}\frac{\mathrm{d}}{\mathrm{d}r}\hfill \\ \hfill & -\frac{m^2/r^2}{{\sigma }^2(\sigma +\stackrel{\sim }{\chi }{m}^2/r^2)}D(\sigma ,m,r)]\left(r^2{\rho }_0\delta v_r\right)=0,\hfill \end{array}$$(21)

where we introduced the pressure scale height,

$$\begin{array}{c}\hfill H=-{\left(\frac{\mathrm{d}log{p}_0}{\mathrm{d}r}\right)}^{-1}>0.\end{array}$$(22)

The operator D is defined similarly to the local elevator MTI dispersion relation, Eq. (2), though adapted to the spherical geometry at hand,

$$\begin{array}{c}\hfill D(\sigma ,m/r)={\sigma }^3+\stackrel{\sim }{\chi }(r)\frac{m^2}{r^2}{\sigma }^2+N_0^2(r)\sigma -\stackrel{\sim }{\chi }(r)\frac{m^2}{r^2}{M}_0^2(r).\end{array}$$(23)

If m/r is formally replaced by a horizontal wavenumber k_x within the equation D(σ,m/r) = 0, we recover the original dispersion relation from Balbus (2000), Eq. (2), for local elevator modes (i.e. d/dr = 0). Moreover, if the derivative operator, d/dr, is locally approximated by a radial wavenumber, k_r, and in the limit k_rH ≫ 1 (which may be valid for modes with shortly varying radial structure), we then regain from Eq. (21) the dispersion relation for arbitrary local 2D perturbations (Parrish et al. 2012, Eq. (18), for vanishing viscosity).

3.2. Fast thermal diffusion limit and Sturm-Liouville problem

Unfortunately, the eigenproblem, Eq. (21), cannot be recast as a Sturm-Liouville ODE, which prevented us from applying results from this theory. Further theoretical progress was then only achievable at the cost of an additional hypothesis. Fast thermal diffusion (i.e. $\stackrel{\sim }{\chi }{m}^2/r^2\gg \sigma ,M_0,N_0$) is a common limit often used in studies of magnetised buoyancy instabilities, such as the MTI and the HBI (Balbus 2000; Quataert 2008). Under this assumption, Eq. (21) can formally be reworked into

$$\begin{array}{c}\hfill [\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}+\frac{1}{H}\frac{\mathrm{d}}{\mathrm{d}r}-\frac{m^2}{r^2}]\left(r^2{\rho }_0\delta v_r\right)=-\frac{m^2}{r^2}\frac{M_0^2(r)}{{\sigma }^2}\left(r^2{\rho }_0\delta v_r\right).\end{array}$$(24)

For our ICM model however, we have ℓ_χ ∼ H ≲ R_vir, therefore $\stackrel{\sim }{\chi }{m}^2/r^2\gg M_0$ also means ℓ_χ⁻² ∼ H⁻² ≪ m²/r². In other words, the limit of fast thermal diffusion corresponds to modes with length scales much smaller than the pressure scale height, H,

$$\begin{array}{c}\hfill \left|\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}\left(r^2{\rho }_0\delta v_r\right)\right|\sim \left|\frac{m^2}{r^2}\left(r^2{\rho }_0\delta v_r\right)\right|\gg \left|\frac{1}{H}\frac{\mathrm{d}}{\mathrm{d}r}\left(r^2{\rho }_0\delta v_r\right)\right|.\end{array}$$(25)

This fact is later used to approximate the growth rate, σ, with a WKBJ procedure, with m⁻¹ as the small control parameter. The first-order derivative term can then be dropped in Eq. (24), which can directly be rewritten as a Sturm-Liouville equation,

$$\begin{array}{c}\hfill [\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}-\frac{m^2}{r^2}]\left(r^2{\rho }_0\delta v_r\right)=-ϵ\frac{m^2}{r^2}\frac{M_0^2(r)}{M_1^2}\left(r^2{\rho }_0\delta v_r\right),\end{array}$$(26)

with homogeneous Dirichlet BCs, $\delta v_r(r_1)=\delta v_r(r_2)=0$. The weight function is −(m²/r²)(M₀²(r)/M₁²), with M₁ the MTI frequency at the innermost radius, and the eigenvalue is ϵ = M₁²/σ². The structure of this problem ensures a discrete set of real eigenvalues, ϵ, that we order increasingly: ϵ_n < ϵ_n + 1 for n ∈ ℕ. With this convention, the smallest Sturm-Liouville eigenvalue, ϵ₀, is reached for the highest growth rate, σ₀. Besides, the n^th eigenfunction associated to the eigenvalue ϵ_n possesses n zeroes in the corresponding solution domain (excluding those at the boundaries).

In Appendix A, we make use of the explicit Sturm-Liouville form of Eq. (26), in the limit of fast thermal diffusion, along with known results from the theory, to further constrain both the growth rate and the radial structure of the global MTI eigenmodes. Here, we simply state the obtained results, which are twofold. First, the growth rate, σ, of any global eigenmode is bounded by the maximum MTI frequency, M₁, in the domain, which is the growth rate of the fastest local elevator mode. Second, the global MTI eigenmodes are exactly the functions that maximise their power weighted by the maximum growth rate of local MTI modes, quantified by

$$\begin{array}{c}\hfill R\left[u\right]={\displaystyle {\int }_{r_1}^{r_2}}\frac{m^2}{r^2}\frac{M_0^2(r)}{M_1^2}{\left|u\right|}^2\mathrm{d}r,\end{array}$$(27)

under the constraint that the quadratic form,

$$\begin{array}{c}\hfill Q\left[u\right]={\displaystyle {\int }_{r_1}^{r_2}}(\frac{m^2}{r^2}{\left|u\right|}^2+{\left|\frac{\mathrm{d}u}{\mathrm{d}r}\right|}^2)\mathrm{d}r,\end{array}$$(28)

is kept constant. The latter Eq. (28) induces a natural subsequent norm, Q, for the MTI eigenfunctions, which is used to normalise any global MTI modes presented in what follows, so that $Q\left[r^2{\rho }_0\delta v_r\right]=1$. This mathematical optimality result is well in line with the physical phenomenology of the MTI eigenfunctions, dictated by the presence of a turning point at a critical finite radius, which we describe next.

3.3. Physics of the global MTI eigenmodes

3.3.1. Turning points – global and local eigenmodes

We first recall that M₀(r) is usually a decreasing function of the radius in the outskirts of galaxy clusters. This is also the case in our model of ICM atmosphere, described in Sect. 2, whose overall MTI frequency profile is displayed in Fig. 1. The maximum local growth rate, M₁, is thus located at the inner boundary, r₁, and reached for the local MTI elevator mode there. On the contrary, the (smallest) MTI frequency, M₂, is reached at the outermost radius, r₂. To better qualitatively describe the behaviour of the global MTI eigenmodes, we rewrote the single governing ODE, in the limit of fast thermal diffusion, under Schrödinger form,

$$\begin{array}{c}\hfill [\frac{{\mathrm{d}}^2}{\mathrm{d}{r}^2}-\frac{m^2}{r^2}(1-\frac{M_0^2(r)}{{\sigma }^2})]\left(r^2{\rho }_0\delta v_r\right)=0.\end{array}$$(29)

From the results obtained in Appendix A, we learnt that the growth rate, σ, lies somewhere between 0 and M₁. So eigenmodes with σ > M₂ present a turning point at some critical radius, r_c, defined as the only solution to M₀(r) = σ in r ∈ [r₁,r₂]. On the left (r < r_c) of the turning point, the modes oscillate in space with a varying wavenumber, $k_r(r)=\sqrt{E-V_{\mathrm{ftd}}}$, because their energy, E = (m²/r²)(M₀²(r)/σ²), is higher than the simplified potential, V_ftd = m²/r² (where the subscript ⋅_ftd stands for fast thermal diffusion). This radial wavenumber, k_r(r), is displayed on the top panel in Fig. 2. On the right (r > r_c) of the turning point though, the modes are evanescent. Physically, beyond the turning point, the local MTI frequency is no longer large enough to support the efficient growth of the global MTI mode, which has no other avenue but to decay there. For a given azimuthal order, m, the higher the mode number, n, the smaller the growth rate, σ, and hence the further away the turning point (and the wider the region of allowed growth, until it reaches the whole domain). As illustrated on the bottom panel in Fig. 2, which depicts a representative numerically computed MTI eigenmode, those modes are necessarily located at smaller radii, where the radial wavenumber is positive: k_r > 0. This phenomenological picture is very well in line with the radial structure of these global eigenmodes, deduced from the Sturm-Liouville theory (details are given in Appendix A): the MTI modes maximise their power weighted by the local MTI growth rate, Eq. (27). Global MTI modes, that are nowhere evanescent, could also develop with σ < M₂. Given their slower growth rate but greater radial extent, it is unclear whether such modes are dynamically important.

Fig. 2. Illustrative numerical global MTI eigenmode with m = 50 and n = 20. Top: varying radial wavenumber, $k_r=\sqrt{E-V_{\mathrm{ftd}}}$, as a function of the radius, r. Bottom: a high-order global MTI mode. The turning point, rc, of the mode is indicated by black dashed vertical lines. It can be found as the single solution to the equation kr(r) = 0, indicated by the grey dashed horizontal line on the top panel.

This result sheds light on the connection between local modes and global modes of the linear MTI: at a given radius, the fastest growing local mode, known as an elevator mode because it has no radial variation, is just a local zoom-in on the single global eigenmode that features a turning point at the same radius. Around this turning point, the global mode exhibits no further radial variation (similarly to the associated local elevator mode) since it becomes evanescent. Mathematically speaking, this is due to locally having k_r(r_c) = 0. On the contrary, global MTI eigenmodes that still exhibit radial structure at a given radius grow necessarily slower than the single mode which features a turning point at the same radius. Locally, these slower growing global modes correspond to radially varying local modes (in contrast to elevator modes). This phenomenological discussion is illustrated in Fig. 2. This reasoning was in fact inspired by the study of the global linear MRI (LFF15). In both cases, the fastest growing local modes exhibit no radial variation (Balbus 2001). For the MRI, they are called channel rather than elevator modes. Such an analogy could be drawn because both Schrödinger equations, Eq. (18) in LFF15 for the MRI, and Eq. (29) here for the MTI (in the limit of fast thermal diffusion), share a very similar structure. As a final note, we stress that a similar phenomenology could describe the localisation of global (possibly compositional) HBI modes in the ICM (Latter & Kunz 2012; Berlok & Pessah 2016), though a formal analysis of the kind led here would be necessary to confirm it. In these studies however, the Braginskii viscosity further helped to confine the HBI eigenfunctions to the lowest altitudes of the ICM atmosphere (Latter & Kunz 2012).

3.3.2. WKBJ solutions

To conclude this section on analytical results, we provide convenient WKBJ approximations to the growth rates of the solutions of the Schrödinger Eq. (29) in the limits of either large azimuthal orders, m/n → ∞, or large mode numbers, n/m → ∞.

Starting with the first limit, m/n → ∞, we found that such global MTI modes are the most unstable and feature a turning point at a radius, r_c, very close to the inner boundary, r₁. Given the decreasing profile of MTI frequency, M₀(r), the solution is wavy between r₁ and r_c, beyond which it becomes evanescent. A standard WKBJ estimate within the wavy region, r ∈ [r₁,r_c], combined with the quantification conditions issued from the inner and outer BCs, used in the limit m/n → ∞, leads to the approximated growth rate,

$$\begin{array}{c}\hfill {\sigma }_n(m)\approx M_1[1-{\left(\frac{3}{2\sqrt{2}}\frac{(n+\frac{3}{4})\pi }{m}{\left|\frac{\mathrm{d}log{M}_0}{\mathrm{d}logr}\right|}_{r=r_1}\right)}^{\frac{2}{3}}].\end{array}$$(30)

The details of this derivation are given in Appendix B.

Next, we also provide an approximated growth rate for global MTI modes, in the limit of large mode numbers, n/m → ∞. These modes are nowhere evanescent in our domain (i.e. σ < M₂). We show in Appendix B that the growth rates of these slowly growing MTI eigenmodes can be approximated as

$$\begin{array}{c}\hfill {\sigma }_n(m)\approx \frac{m}{n\pi }{\displaystyle {\int }_{r_1}^{r_2}}\frac{M_0}{r}\mathrm{d}r.\end{array}$$(31)

The integral on the right-hand side is easily computed numerically and its value depends only on the background structure, not on the particulars of any mode. In the next section, these analytical WKBJ approximations, Eqs. (30) and (31), are compared, in the relevant limits, m/n → ∞ and n/m → ∞, to the growth rates of numerically computed global MTI eigenmodes.

4. Numerical results

In this section, we describe our numerical work. First, we numerically computed the solutions of the global linear MTI eigenproblem, thanks to the Chebyshev pseudo-spectral solver presented in Sect. 2. Then, we performed direct numerical simulations of the MTI linear phase with the code IDEFIX, to cross-check the validity of the two methods.

4.1. Numerical solutions

We present here some numerical solutions to the system of Eqs. (11), (12), (13), (14), and (15) for the unknowns $\delta \rho /{\rho }_0,\delta v_{\phi },\delta v_r,\delta T/T_0,\delta A/B_0$, within a finite domain extending from r₁ = 0.3 to r₁ = 1. We stress that the equations solved here do not make use of the sound-proof approximation, or of the fast thermal diffusion limit, previously introduced in Sect. 3 to analytically probe the global MTI eigenmodes. A variety of numerically computed MTI eigenmodes are depicted in Fig. 3. All of them are normalised so that $Q[r^2{\rho }_0\delta v_r]=1$. According to the Sturm-Liouville theory discussed in Sect. 3, the mode number, n, at fixed m, is also the number of nodes displayed by the radial velocity component of the corresponding eigenmode. Similarly, the different eigenfunctions become evanescent beyond a certain radius, given by the mode turning point. The location of this turning point depends on the overall mode growth rate: most unstable modes are indeed more localised in the innermost regions of the domain, and they present a very simple radial structure when they are associated with the fundamental mode number, n = 0 (not shown). In contrast, more slowly growing modes extend further out in the atmosphere, and they exhibit richer radial structures when they take higher mode numbers, n > 0.

Fig. 3. Numerical solutions to the eigenproblem Eqs. (11), (12), (13), (14), and (15). From top to bottom: global MTI eigenfunctions for azimuthal order m = 30, 50, 100, respectively. From left to right: global MTI eigenfunctions for mode number n = 3, 12, respectively. From top to bottom in each subfigure: perturbations of density, azimuthal velocity, radial velocity, and temperature as a function of the radius. The potential vector perturbation is not shown because it simply relates to that of radial velocity through Eq. (15). Real parts of the modes are in plain lines, while imaginary parts are in dashed lines. The part whose corresponding line type is absent is zero. The modes oscillate in space up to their respective turning points, rc, and evanescent beyond.

Next, we focused on the dispersion relation of the global linear MTI. On the top panel in Fig. 4, we show the growth rates, σ_n, of the first six modes as a function of the azimuthal order, m. Appreciable growth of all the illustrated MTI modes occurs for m ≳ 30. Interestingly, the growth rates of all modes seem to converge to each other in the limit of large azimuthal orders, m/n → ∞, where they reach the asymptote σ → M₁. This behaviour is only possible because the present study does not include any other diffusive effect, such as isotropic viscosity or magnetic resistivity, nor magnetic tension, that would otherwise stabilise the small scales. In Fig. 4 (inset), we plot the WKBJ analytical approximation, Eq. (30), for modes n = 0, 2, and 4, on top on the numerically obtained growth rates. We see that they agree very well with each other, especially in the limit of large azimuthal orders, m/n → ∞, corresponding to modes with turning points, r_c, very close to the inner radius, r₁.

Fig. 4. Dispersion relations of the global MTI eigenmodes. Top: growth rates of the first six MTI modes as a function of the azimuthal order, m. The curve σ = M1 is indicated in dotted dark line. It seems to be a common asymptote to all curves, σn(m). Inset: same growth rates, σn, again (one over two for the sake of visibility though), together with their WKBJ approximations, Eq. (30), in dashed lines. Bottom: growth rates of the MTI eigenmodes with m = 100 and m = 50 as a function of the mode number, n, in plain lines. Their WKBJ approximations, Eq. (31), are shown too, in dashed lines. The curve σ = M2 is indicated in dotted grey line.

Additionally, we probed the behaviour of the global MTI dispersion relation as a function of the mode number, n, at fixed azimuthal order, m. On the bottom panel in Fig. 4, we plot the growth rates σ_n(m = 100) and σ_n(m = 50) for mode numbers up to n = 1000. Such modes exhibit slower growth rates (σ_n → 0 as n → ∞, as expected from the WKBJ approximation) but richer radial structure and further extent within the ICM atmosphere, since the critical radii of their turning points are far away from the inner boundary, r₁. In fact, such turning points for global MTI modes with growth rates, σ_n, lower than M₂ (below the horizontal grey line on the bottom panel in Fig. 4) are outside the numerical solution region, [r₁,r₂]. On the same panel, the growth rates, Eq. (31), approximated with a WKBJ procedure, are shown for both m = 100 and m = 50. The WKBJ approximation, Eq. (31), perfectly fits the numerically obtained growth rates, in the limit of large mode numbers, n/m → ∞. In particular, it captures the right asymptotic behaviour, σ_n ∼ n⁻¹.

4.2. Direct simulations

We present three different numerical simulations of the MTI linear stage, performed with the code IDEFIX, described in Sect. 2. We first discuss further the numerical implementation of the BCs, along with the diagnostics, used in these simulations. Then, we assess how well global MTI eigenmodes, and their growth rates, can be reproduced in these simulations, by comparing them to the numerical eigenmodes presented in Sect. 4.1.

4.2.1. Boundary conditions and diagnostics

In the light of the common radial structure shared by all global MTI modes seen in Fig. 3, we describe the BCs imposed in all our simulations. For all fields, we took periodic BCs in the azimuthal direction. More care was needed to choose the BCs in the radial direction. At the inner boundary for instance, the radial velocity, the temperature, and the density perturbations of all modes follow a homogeneous Dirichlet condition. In contrast, the fractional pressure perturbation is non-zero at the inner boundary, while its first derivative must be zero, according to Eq. (13). Therefore, we paid special attention in preserving these specific traits of the modes when designing our inner BCs: we either imposed symmetry (i.e. zero gradient) or anti-symmetry of the perturbations around their equilibrium values (which are known from the initial HSE condition). Regarding the outer boundary, we simply let the different fields take their equilibrium values. Such outer BCs are only suited for the fastest growing modes, which are evanescent beyond a given radius, but not for all MTI eigenmodes. This implementation ensures that the BCs disturbed neither the symmetry of the MTI eigenfunctions, nor the initial HSE on top of which they developed.

Moving on to the diagnostics, we systematically probed the radial structure of all the components (velocities, pressure, temperature, and density perturbations) of the linear modes developing in the simulation domain, through a Fourier analysis of their azimuthal dependency,

$$\begin{array}{c}\hfill \delta f_m(r)=\frac{1}{2\pi }{\displaystyle {\int }_0^{2\pi }}exp(\mathrm{i}m\phi )\delta f(r,\phi )\mathrm{d}\phi ,\end{array}$$(32)

where δf is either $\delta v_r,\delta v_{\phi },\delta T$ or δρ. The power contained in each component of the mode was then measured as

$$\begin{array}{c}\hfill {⟨\delta f_m⟩}_{\mathcal{V}}=\frac{3}{r_2^3-r_1^3}{\displaystyle {\int }_{r_1}^{r_2}}\left|\delta f_m(r)\right|r^2\mathrm{d}r.\end{array}$$(33)

The numerical growth rate of each mode was obtained by following the temporal evolution of these different powers, over a time window during which this evolution was linear (i.e. after efficient growth of the mode started and before it saturated).

In the next subsection, we present three different simulations, all of them initialised with very low amplitude perturbations, of the order of 10⁻¹⁰. The first run used a random white noise on the velocity field and extends from 0 to 2π in the azimuthal direction. The two other runs were seeded with a single global MTI eigenmode each and extended from 0 to π, so as to efficiently isolate the growth of a given mode with enough azimuthal resolution.

4.2.2. Numerical growth rates and MTI eigenfunctions

We plot a selection of numerical growth rates on top of the theoretical dispersion relations for the fastest growing MTI modes in Fig. 5. The chosen numerical growth rates, σ_num, corresponding to the black diamonds, are measured independently from each other in the first simulation, which was initialised with a random white noise. They are required to match the theoretical growth rates, σ_th, computed with the Chebyshev pseudo-spectral method, to within ζ ≲ 1% (the error was computed as ζ = |σ_num−σ_th|/σ_th). Additionally, for each of them, we visually checked that the radial structure of the associated mode is easily identifiable and that it matches well the structure of the theoretical eigenmode. At the end of this filtering process, 17 fastest growing MTI eigenmodes (with n = 0 and azimuthal orders between 15 ≤ m ≤ 62) are identified in the first simulation. Modes with m ≲ 14 can not grow efficiently because the support of their radial eigenfunctions extends beyond the outer boundary, r₂ = 0.7, used in the simulations. Similarly, we demonstrate below that modes with m ≳ 70 are not sufficiently well resolved to grow at the right rate, at least in the first simulation seeded with random white noise on the velocity field.

Fig. 5. Comparison between the analytical dispersion relation (in plain lines) and the numerically measured growth rates (markers). Black diamonds are the numerical growth rates of the most unstable modes with n = 0, found in the first simulation seeded with low amplitude random white noise on the velocity field. They match the theoretical dispersion relation to within ζ ≲ 1%. The blue star depicts the MTI eigenmode with m = 92 and n = 1, which could be successfully isolated in the same simulation. The purple circle stands for the numerical growth rate of the MTI eigenmode with m = 101 and n = 0. It agrees with the theoretical dispersion relation to within ζ ≈ 2.0% only. The green triangle represents the numerical growth rate of the exact same mode, but in a simulation seeded with a single mode. The error is very small, ζ ≈ 0.027%.

Isolating higher order MTI modes with n > 0 in the first simulation is trickier, because they grow slower and are thus always supplanted by n = 0 modes. Nevertheless, we managed to identify a persisting MTI eigenmode n = 1,  m = 92 in this simulation. The time evolution of this mode, δv_r,92(r), is roughly self-similar (not shown), despite growing by a factor ∼10⁶ through the course of the simulation. Its growth rate matches that from the theoretical dispersion relation to within ζ ≈ 0.9%, so it is added as a blue star in Fig. 5. However, comparing this numerical eigenmode with m = 92 and n = 1 against the corresponding theoretical mode reveals that its radial structure is not completely resolved at the innermost radii (not shown). Therefore, we also produced a simulation with a single seed of the same mode in order to provide an unequivocal comparison basis. For this simulation, we decreased the azimuthal extent of our solution domain from 2π to π, while keeping 4096 cells in this direction, to ensure that the mode is sufficiently well resolved and that it is not supplanted by the mode n = 0 (which can be triggered by discretisation errors) in the course of the simulation. In Fig. 6, the quantity δv_r,92(r) is plotted, at two different times during the simulation seeded with a single mode. It is normalised by its current maximum. This time, the mode evolution is completely self-similar across time, despite increasing by a factor ∼10⁷. Accordingly, the numerically measured growth rate agrees with that from the theoretical dispersion relation to within ζ ≈ 0.002%, a thousand times better than in the first simulation.

Fig. 6. Radial velocity of the global MTI eigenmode with m = 92 and n = 1 in the simulation seeded with a single mode. Radial eigenfunctions, δvr,92(r), are shown at two different times and normalised by their maximum values. Despite growing by a factor ∼107, the eigenmode still lies almost perfectly on top of the initial curve, and thus matches very precisely the theoretical eigenfunction.

Finally, we also added by hand the numerical MTI eigenmode with m = 101 and n = 0, identified in the very first simulation, as a purple circle in Fig. 5. This time, while the radial structure of the numerical mode, δv_r,101(r), matches quite well that of the theoretical MTI eigenmode (not shown), its numerical growth rate is slightly too low with respect to the dispersion relation. The discrepancy is of the order of ζ ≈ 2.4%. We therefore performed a dedicated simulation seeded with a single MTI mode (m = 101,n = 0) and with an azimuthal extent reduced again from 2π to π. The numerical growth rate measured in this simulation agrees with that from the theoretical dispersion relation to within ζ ≈ 0.027%, a hundred times better than for the randomly seeded simulation. This data point is shown as a green triangle in Fig. 5. We deduced that the azimuthal resolution in the very first run is not high enough to satisfactorily model the development of the fastest growing MTI eigenmodes with m ≳ 70, at least when other modes with equivalent azimuthal orders are developing concomitantly. Altogether, these results provide a valuable cross-validation of both the linear theory and the numerical codes.

5. Conclusion

In this paper, we have investigated the global linear stability properties of dilute stratified plasma subject to anisotropic transport of heat. Such a study might be relevant to the intracluster medium (ICM), which seemingly meets all the required conditions for the development of the magneto-thermal instability (MTI), and which is already known to sustain mild turbulence (Hitomi Collaboration 2016, 2018; XRISM Collaboration 2025a,b).

5.1. Summary

The main objective of this work was to unveil the deeper connections existing between local modes and their global counterparts, through a global linear analysis of the MTI in the ICM. A similar linear study was previously performed for the magneto-rotational instability (MRI) in the context of accretion discs by LFF15. We were able to build upon this earlier analysis, on account of the very similar linear behaviours shared by the MTI and the MRI (Balbus 2001).

Our main findings are summarised below:

• The global linear analyses of both the MRI and the MTI share many similarities. The most striking feature is the very similar overall eigenmode structure for these two instabilities. In both cases, the fastest growing local modes lack radial structure (they are called channel and elevator modes, in the case of the MRI and the MTI, respectively). These modes correspond, for both instabilities, to the evanescent parts of global eigenmodes; while the slowest growing local modes, which feature non-trivial radial structures, are in fact segments of the same global mode, but at smaller radii (LFF15). This feature may explain why local and global models of ICM atmosphere produce similar magneto-thermal turbulence (PL22a; PL22b; KR25), independently of the scale that they probe.

• From a numerical point of view, we successfully computed numerical solutions to the global MTI eigenproblem with our Chebyshev pseudo-spectral solver. These solutions were compared with direct numerical simulations of the MTI linear stage with the new generation finite-volume code IDEFIX, including a Braginskii module for anisotropic diffusion. Therefore, the results presented in this paper provide avaluable benchmark for any astrophysical fluid dynamics code implementing the Braginskii heat flux in spherical geometry.

This work is an additional step towards a better understanding of how local and global approaches to linear stability theory intertwine. It is our hope too that this paper will help to pave the way for a clearer view of the exact role played by global linear eigenmodes, both in terms of preferred amplitudes and locations, at turbulent saturation. More work is however needed to draw a complete and generic picture of the non-linear saturation of such instabilities. We briefly comment on this issue in the next subsection.

5.2. Caveats and future improvements

Although the model developed in this work could satisfactorily address the questions raised in the introduction, it could still be improved in many ways. First, an obvious refinement of the model would consist in including neglected non-ideal transport effects, such as magnetic resistivity or (possibly anisotropic) viscosity, or magnetic feedback on the flow through the Lorentz force. However, given the specific parameter regime of the ICM in terms of plasma beta, and regular and magnetic Prandtl numbers, we anticipate that these effects are not critical and would only modestly affect the results presented here.

Other clear caveats of the current analysis are related to our approximations of the modes being 2D, and of the magnetic configuration being purely toroidal, $\mathit{B}=B_0{\widehat{\mathit{e}}}_{\phi }$. The main advantage of these assumptions is to leave the model analytically tractable, while still capturing two of the main features of the ICM in galaxy clusters, namely stratification and anisotropic transport of heat, which are critical for the MTI. Going 3D and considering the full sphere, (r,θ,φ), by adding the so far neglected polar dependency and component, ${\widehat{\mathit{e}}}_{\theta }$, would be an easy way to implement more complex magnetic geometries while preserving the thermal equilibrium of spherically symmetric background profiles. Such a development on the full sphere would also enable us to probe the polar dependency, out of the equatorial plane, of the global MTI modes in the ICM. However, since the Braginskii heat flux is a non-linear operator coupling the magnetic field to the temperature field, a single mode analysis, on the basis of the (vector) spherical harmonics, as presented by Choudhury & Sharma (2016) in the case of the thermal instability, is precluded for the 3D linear MTI. Such a linear analysis on a basis of spherical harmonics would result in a 2D eigenvalue problem, whose solution is intrinsically more demanding than the 1D problem considered in this work. Such an extension of our model is outside the scope of the current study.

Finally, the global linear analysis of the MTI in the ICM carried out in this paper is completely oblivious to the specific energetic framework recently developed by KR25. In that work, KR25 showed that the notion of available potential energy (APE) is anchored at the very heart of the saturation mechanism of magneto-thermal turbulence. This phenomenology was inspired by fundamental studies in the context of stratified fluid dynamics (Middleton & Taylor 2020; Tailleux 2024), which suggested that having a net positive conversion rate of APE is a necessary condition for the linear development of double-diffusive instabilities, such as thermohaline or semi-convection, in the ocean. Therefore, it is only natural to wonder whether the same physics, namely having a net positive conversion rate of APE, could also be responsible for the triggering of linear magnetised buoyancy instabilities in dilute plasma, or not. This research direction represents an interesting avenue to better constrain the energetics of magnetised buoyancy instabilities in the ICM and to further understand how energetics and linear theory interlock. Specific work should be devoted to such questions in the future. Similarly, the role played by the most unstable MTI modes at saturation remains unclear. Future studies carefully examining whether linear theory can be of any help for understanding specific features of magneto-thermal turbulence at non-linear saturation, especially in terms of the turbulence preferred length scale and intensity (PL22a; PL22b; KR25), are therefore needed.

Acknowledgments

The authors are thankful to Thomas Berlok, Joshua Brown, Jean-Baptiste Durrive, Thomas Jannaud, Gordon Ogilvie, Prakriti Pal Choudhury, François Rincon, and Remi Tailleux, for many useful discussions which helped to improve the quality of this work. It is a pleasure to thank Geoffroy Lesur for sharing a version of the balance scheme that he had previously developed in IDEFIX. J.M.K acknowledges funding through the grant EUR TESS N°ANR-18-EURE-0018 in the framework of the Programme des Investissements d’Avenir, which funded a three-month stay in Cambridge during his PhD. J.M.K is grateful for the hospitality of the host institution, the DAMTP, where this work was initiated during summer 2024. This work was granted access to the HPC resources of CALMIP under the allocation 2023/2024-P16006, J.M.K thanks the associated support team for their work. The authors would also like to thank François Rincon for thoroughly reading a draft of the manuscript.

References

Appendix A: Sturm-Liouville theorems

In this appendix, we prove the two results stated in Sect. 3 about the growth rates and the radial structures of global MTI eigenmodes. This is possible thanks to a specific result of the Sturm-Liouville theory, known as the Rayleigh-Ritz rule. Before doing so, though, we raise a point of caution that intrinsically limits the scope of these results to the most unstable MTI eigenmodes, in the limit of fast thermal diffusion.

Formally speaking, unveiling an approximation to the pristine Eq. (21) under the form of a linear Sturm-Liouville equation, Eq. (26), does not automatically ensure that all known results of the Sturm-Liouville theory apply per se to the global linear theory of the MTI, precisely because Eq. (26) is only an approximation. In the relevant limit, one might hope however that conclusions relevant to the problem at hand can still be found in the results of the Sturm-Liouville theory. In the remaining part of this section, we do as if this was the case, even though further mathematical work would be needed to justify more rigorously the approach taken here.

A.1. Growth rate

First, we derived bounds for the Sturm-Liouville eigenvalue, ϵ. We multiplied Eq. (26) by ${\left(r^2{\rho }_0\delta v_r\right)}^{\ast }$ where ⋅^* denotes the complex conjugate, and integrated over the solution domain, [r₁,r₂]. The first term on the left-hand side was integrated by parts, and the homogeneous Dirichlet BCs were used to cancel the subsequently arising border terms, so as to obtain

$$\begin{array}{c}\hfill ϵ=\frac{{\displaystyle {\displaystyle {\int }_{r_1}^{r_2}}\frac{m^2}{r^2}{\left|r^2{\rho }_0\delta v_r\right|}^2\mathrm{d}r}}{{\displaystyle {\displaystyle {\int }_{r_1}^{r_2}}\frac{m^2}{r^2}\frac{M_0^2(r)}{M_1^2}{\left|r^2{\rho }_0\delta v_r\right|}^2\mathrm{d}r}}+\frac{{\displaystyle {\displaystyle {\int }_{r_1}^{r_2}}{\left|\frac{\mathrm{d}}{\mathrm{d}r}\left(r^2{\rho }_0\delta v_r\right)\right|}^2\mathrm{d}r}}{{\displaystyle {\displaystyle {\int }_{r_1}^{r_2}}\frac{m^2}{r^2}\frac{M_0^2(r)}{M_1^2}{\left|r^2{\rho }_0\delta v_r\right|}^2\mathrm{d}r}}.\end{array}$$(A.1)

From this expression, we saw that the Sturm-Liouville eigenvalue, ϵ, must be positive, provided that the domain is everywhere unstable to the MTI (i.e. M₀²(r) > 0 for r ∈ [r₁,r₂]). This requirement is, in fact, too strong: the eigenvalue, ϵ, is positive as long as the common denominator in both fractions is positive. In principle, this specific case remains accessible (for some eigenmodes at least), even though M₀²(r) < 0 in some regions of the domain. For such modes (which make the denominators positive), another bound can be found since the time scale M₁ is the maximum value of M₀(r) over the domain r ∈ [r₁,r₂]. In this case, the numerator of the first fraction on the right-hand side of Eq. (A.1) is automatically greater than the denominator of the same fraction, and ϵ must be greater than unity, since the second fraction is also positive. In other words, ϵ > 1 means that the growth rate of the global MTI eigenmodes cannot be higher than the growth rate of the most unstable local elevator mode in the solution domain (i.e. the maximum MTI frequency, M₁). This feature is satisfying from a physical point of view, and is further exploited in Sect. 3 to relate local and global MTI modes.

A.2. Radial structure

Next, we constrained the structure of the solutions to the Sturm-Liouville Eq. (26), which are the global MTI eigenmodes in the limit of fast thermal diffusion. A specific result of the theory, the Rayleigh-Ritz rule, provides a physical relationship linking the global structure of the different MTI eigenmodes to the MTI frequency profile, M₀(r). This criterion states that the Sturm-Liouville eigenvalue, ϵ, is the quotient between the two quadratic forms,

$$\begin{array}{cc}\hfill Q\left[u\right]& ={\displaystyle {\int }_{r_1}^{r_2}}(\frac{m^2}{r^2}{\left|u\right|}^2+{\left|\frac{\mathrm{d}u}{\mathrm{d}r}\right|}^2)\mathrm{d}r,\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill R\left[u\right]& ={\displaystyle {\int }_{r_1}^{r_2}}\frac{m^2}{r^2}\frac{M_0^2(r)}{M_1^2}{\left|u\right|}^2\mathrm{d}r,\hfill \end{array}$$(A.3)

which satisfactorily corroborates Eq. (A.1). Additionally, the rule states that the solutions of the Sturm-Liouville Eq. (26) are also solutions of the variational problem consisting in minimising the ratio Q[u]/R[u], amongst all functions, u, that satisfy the homogeneous Dirichlet BCs. The so obtained eigenmode, u₀, has then a minimal associated eigenvalue, ϵ₀, and the next smallest eigenvalue is obtained through the minimisation of the same ratio, under the additional constraint that

$$\begin{array}{c}\hfill {\displaystyle {\int }_{r_1}^{r_2}}\frac{m^2}{r^2}\frac{M_0^2(r)}{M_1^2}{u}_0^{\ast }u\mathrm{d}r=0,\end{array}$$(A.4)

and so on and so forth. We point out that minimising Q[u]/R[u] is equivalent to maximising R[u] under the constraint that Q[u] is kept constant. In other words, the global MTI eigenmodes, in the limit of fast thermal diffusion and at fixed azimuthal order, m, are exactly the n-indexed functions that maximise their power weighted by the local maximum MTI growth rate, represented by the quantity R[u_n]. These functions therefore maximise their global growth rate, σ_n, under the constraint of a constant weighted Sobolev norm (e.g. Q[u_n] = 1). The Sturm-Liouville theory indicates that the eigenfunctions, u_n, form an orthogonal basis under the weighted inner product defined by Eq. (A.4). Physically, this suggests that the MTI eigenmodes are tapping into the background profile of maximum local growth rate at different locations, in phase quadrature loosely speaking.

Appendix B: Analytical WKBJ approximations

In this appendix, we describe the WKBJ procedure used to deduce the analytical approximations, Eqs. (30) and (31), of the growth rates, σ_n(m). We also detail the derivations when necessary. We take m⁻¹ as our small WKBJ parameter.

We first focused on the fastest growing global MTI modes, which feature a turning point very close to the inner boundary, r₁. The critical radius, r_c, of the turning point is defined as the single solution to the equation M₀(r) = σ. Given the decreasing profile of MTI frequency, M₀(r), the solution is wavy between r₁ and r_c, beyond which it becomes evanescent. A standard WKBJ estimate is then used within the wavy region, r ∈ [r₁,r_c],

$$\begin{array}{c}\hfill r^2{\rho }_0\delta v_r\propto cos[{\displaystyle {\int }_r^{r_{\mathrm{c}}}}\frac{m}{s}{(\frac{M_0^2(s)}{{\sigma }^2}-1)}^{\frac{1}{2}}\mathrm{d}s+\frac{\pi }{4}].\end{array}$$(B.1)

The π/4 phase shift results from matching the WKBJ ansatz at the turning point, when the outer radius, r₂, is sent to infinity and the solution asked to vanish there (Riley et al. 2006). The quantification of the growth rate, σ_n, arises from the requirement that the eigenfunctions also cancel at the inner boundary, r₁,

$$\begin{array}{c}\hfill {\displaystyle {\int }_{r_1}^{r_{\mathrm{c}}({\sigma }_n)}}\frac{m}{r}{(\frac{M_0^2(r)}{{\sigma }_n^2}-1)}^{\frac{1}{2}}\mathrm{d}r=(n+\frac{3}{4})\pi ,\end{array}$$(B.2)

where the turning radius, r_c, depends on the eigenvalue, σ_n, itself. Generally, Eq. (B.2) is a non-linear equation of the unknown σ_n, and should therefore be solved numerically. However, we show that it is still possible to approximate the growth rate, σ_n, in the limit of large azimuthal orders, m/n → ∞. We consider the change of variable u = 1 − M₁/M₀. Introducing also λ = M₁/σ > 1, we get

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{1-\lambda }}\left(\frac{M_0^2}{r(u)M_1{M}_0^{\prime }}\right)\frac{\sqrt{{\lambda }^2-{(1-u)}^2}}{1-u}\mathrm{d}u=m^{-1}(n+\frac{3}{4})\pi .\end{array}$$(B.3)

We were interested in λ → 1. We wrote λ = 1 + ε, for 0 < ε ≪ 1. The integration range above becomes very small (between 0 and −ε) and thus we could expand the integrand around u = 0, that is r = r₁ or, equivalently, r = r_c. To leading order in ε, noting that in the integral u ∼ ε and that M₀′< 0, we obtained

$$\begin{array}{c}\hfill \sqrt{2}{\left|\frac{\mathrm{d}log{M}_0}{\mathrm{d}logr}\right|}_{r=r_1}^{-1}{\displaystyle {\int }_{-\epsilon }^0}\sqrt{u+\epsilon }\mathrm{d}u=m^{-1}(n+\frac{3}{4})\pi .\end{array}$$(B.4)

After performing the simple integration, we recovered the approximated WKBJ growth rate, Eq. (30), in the limit of large azimuthal orders, m/n → ∞.

Next, we focused on slow radially extended MTI eigenmodes, in the limit of large mode numbers, n/m → ∞. Such modes do not feature turning points in the solution domain, so a different WKBJ ansatz is used,

$$\begin{array}{c}\hfill r^2{\rho }_0\delta v_r\propto cos[{\displaystyle {\int }_{r_1}^r}\frac{m}{s}{(\frac{M_0^2(s)}{{\sigma }^2}-1)}^{\frac{1}{2}}\mathrm{d}s+\frac{\pi }{2}].\end{array}$$(B.5)

In this case, the π/2 phase shift results from requiring that the solution vanishes at the inner boundary, r = r₁, and the quantification of the growth rate arises from the outer BC,

$$\begin{array}{c}\hfill {\displaystyle {\int }_{r_1}^{r_2}}\frac{m}{r}{(\frac{M_0^2(r)}{{\sigma }_n^2}-1)}^{\frac{1}{2}}\mathrm{d}r=n\pi .\end{array}$$(B.6)

However, in the limit n/m → ∞, M₀/σ_n → ∞ everywhere in the domain too. In this case, the square root appearing within the integrand on the left-hand side of the latter equation is approximated as M₀/σ_n and Eq. (31) directly follows. We note that the WKBJ ansatz, Eq. (B.5), still requires m⁻¹ to be a small parameter. Our result is then only valid for n ≫ m ≫ 1. This is the reason why we used the azimuthal orders m = 50 and m = 100, together with mode numbers as large as n = 10³, on the bottom panel in Fig. 4.

Finally, for completeness, we provide a last analytical formula, which is valid for MTI modes in the limit n/m → ∞, but when the outer boundary of the domain is sent far away and the MTI frequency required to vanish there: M₀ → 0 when r → ∞. Modes in such an ICM atmosphere systematically feature a turning point. Even though their critical radii, r_c, are located very far away from the inner boundary, r₁, the same WKBJ ansatz as in the limit m/n → ∞, Eq. (B.1), and the same quantification condition, Eq. (B.2), can be used. In this case however, we were instead interested in the limit λ⁻¹ = σ/M₁ → 0. We then used the change of variable u = σ/M₀, so that Eq. (B.2) transforms to

$$\begin{array}{c}\hfill {\displaystyle {\int }_{{\lambda }^{-1}}^1}{\left|\frac{\mathrm{d}log{M}_0}{\mathrm{d}logr}\right|}^{-1}\frac{\sqrt{1-u^2}}{u^2}\mathrm{d}u=m^{-1}(n+\frac{3}{4})\pi .\end{array}$$(B.7)

The integral on the left-hand side is controlled by the nearly diverging u⁻² factor, as the lower integration bound is λ⁻¹ → 0. The integral is then dominated by the contribution of a small innermost region around u = λ⁻¹, that is near r = r₁. Assuming that, within this region, the factor |dlog M₀/dlog r| does not vary significantly and that the second factor is controlled by u⁻², we obtained the approximated WKBJ growth rate,

$$\begin{array}{c}\hfill {\sigma }_n(m)\approx M_1\frac{m}{(n+3/4)\pi }{\left|\frac{\mathrm{d}log{M}_0}{\mathrm{d}logr}\right|}_{r=r_1}^{-1}.\end{array}$$(B.8)

Both Eq. (31) and Eq. (B.8) share a similar scaling, ∝n⁻¹, although the physical prefactors in both equations differ.

All Figures

	Fig. 1. Local MTI frequency profile, M0, as a function of the radius, r.
In the text

Fig. 2. Illustrative numerical global MTI eigenmode with m = 50 and n = 20. Top: varying radial wavenumber, $k_r=\sqrt{E-V_{\mathrm{ftd}}}$, as a function of the radius, r. Bottom: a high-order global MTI mode. The turning point, rc, of the mode is indicated by black dashed vertical lines. It can be found as the single solution to the equation kr(r) = 0, indicated by the grey dashed horizontal line on the top panel.

In the text

Fig. 3. Numerical solutions to the eigenproblem Eqs. (11), (12), (13), (14), and (15). From top to bottom: global MTI eigenfunctions for azimuthal order m = 30, 50, 100, respectively. From left to right: global MTI eigenfunctions for mode number n = 3, 12, respectively. From top to bottom in each subfigure: perturbations of density, azimuthal velocity, radial velocity, and temperature as a function of the radius. The potential vector perturbation is not shown because it simply relates to that of radial velocity through Eq. (15). Real parts of the modes are in plain lines, while imaginary parts are in dashed lines. The part whose corresponding line type is absent is zero. The modes oscillate in space up to their respective turning points, rc, and evanescent beyond.

In the text

Fig. 4. Dispersion relations of the global MTI eigenmodes. Top: growth rates of the first six MTI modes as a function of the azimuthal order, m. The curve σ = M1 is indicated in dotted dark line. It seems to be a common asymptote to all curves, σn(m). Inset: same growth rates, σn, again (one over two for the sake of visibility though), together with their WKBJ approximations, Eq. (30), in dashed lines. Bottom: growth rates of the MTI eigenmodes with m = 100 and m = 50 as a function of the mode number, n, in plain lines. Their WKBJ approximations, Eq. (31), are shown too, in dashed lines. The curve σ = M2 is indicated in dotted grey line.

In the text

Fig. 5. Comparison between the analytical dispersion relation (in plain lines) and the numerically measured growth rates (markers). Black diamonds are the numerical growth rates of the most unstable modes with n = 0, found in the first simulation seeded with low amplitude random white noise on the velocity field. They match the theoretical dispersion relation to within ζ ≲ 1%. The blue star depicts the MTI eigenmode with m = 92 and n = 1, which could be successfully isolated in the same simulation. The purple circle stands for the numerical growth rate of the MTI eigenmode with m = 101 and n = 0. It agrees with the theoretical dispersion relation to within ζ ≈ 2.0% only. The green triangle represents the numerical growth rate of the exact same mode, but in a simulation seeded with a single mode. The error is very small, ζ ≈ 0.027%.

In the text

	Fig. 6. Radial velocity of the global MTI eigenmode with m = 92 and n = 1 in the simulation seeded with a single mode. Radial eigenfunctions, δvr,92(r), are shown at two different times and normalised by their maximum values. Despite growing by a factor ∼107, the eigenmode still lies almost perfectly on top of the initial curve, and thus matches very precisely the theoretical eigenfunction.
In the text