TY - JOUR AU1 - Zhang,, Congyao AU2 - Churazov,, Eugene AU3 - Schekochihin, Alexander, A AB - ABSTRACT Buoyant bubbles of relativistic plasma in cluster cores plausibly play a key role in conveying the energy from a supermassive black hole to the intracluster medium (ICM) – the process known as radio-mode active galactic nucleus feedback. Energy conservation guarantees that a bubble loses most of its energy to the ICM after crossing several pressure scale heights. However, actual processes responsible for transferring the energy to the ICM are still being debated. One attractive possibility is the excitation of internal waves, which are trapped in the cluster’s core and eventually dissipate. Here, we show that a sufficient condition for efficient excitation of these waves in stratified cluster atmospheres is flattening of the bubbles in the radial direction. In our numerical simulations, we model the bubbles phenomenologically as rigid bodies buoyantly rising in the stratified cluster atmosphere. We find that the terminal velocities of the flattened bubbles are small enough so that the Froude number Fr ≲ 1. The effects of stratification make the dominant contribution to the total drag force balancing the buoyancy force. Clear signs of internal waves are seen in the simulations. These waves propagate horizontally and downwards from the rising bubble, spreading their energy over large volumes of the ICM. If our findings are scaled to the conditions of the Perseus cluster, the expected terminal velocity is ∼ 100−200 km s−1 near the cluster cores, which is in broad agreement with direct measurements by the Hitomi satellite. hydrodynamics, waves, methods: numerical, galaxies: clusters: intracluster medium, X-rays: galaxies: clusters 1 INTRODUCTION Radio mode of active galactic nucleus (AGN) feedback has been widely accepted as the mechanism to prevent the gas in the cores of galaxy clusters from catastrophic overcooling (see e.g. Fabian 2012; McNamara & Nulsen 2012; Vikhlinin et al. 2014; Soker 2016 for reviews). Clear signs of the intracluster medium (ICM) interaction with radio lobes inflated by radio jets from supermassive black holes (SMBHs) have already been seen in ROSAT images of the Perseus and M87/Virgo clusters (Boehringer et al. 1993, 1995). The very same images provide an estimate (lower limit) of the power needed to inflate these bubbles based on the comparison of the inflation and buoyancy time scales and show that this power is comparable to the gas cooling losses (Churazov et al. 2000). This conclusion has been confirmed with much sharper Chandra images (e.g. McNamara et al. 2000; Fabian et al. 2000) and extended for large samples of clusters (e.g. Bîrzan et al. 2004; Hlavacek-Larrondo et al. 2012). Radio emission from the bubbles implies that at least part (and plausibly most) of the bubble pressure support is provided by cosmic rays (CRs) and magnetic fields. If CRs are confined within the bubbles, then the interaction with the ICM is purely mechanical, mediated by the motion of bubble boundaries. An expanding boundary can launch shocks and sound waves into the ICM, if it expands quickly. For a very rapid expansion, most of the energy goes into a shock-heated gas shell, while the fraction of the energy carried away by an outgoing sound wave is less than ∼12.5 per cent if the boundary expansion is spherically symmetric (Tang & Churazov 2017). In the opposite limit of a relatively slow expansion, most of the energy released inside the bubble is stored in the form of the bubble enthalpy (the sum of internal thermal energy and PV work done by the boundary). Here, we concentrate on the latter scenario, in which the bubble enthalpy accounts for most of the SMBH energy output. Bubbles of relativistic plasma are buoyant in the cluster atmosphere and will move to larger radii (Gull & Northover 1973).Energy conservation arguments imply that much of a bubble’s enthalpy will be transferred to the ICM once the bubble crosses several pressure scale heights (Churazov et al. 2001, 2002; Begelman 2001). While these arguments guarantee high coupling efficiency of the radio-mode AGN feedback, they do not depend on the properties of the ICM and, therefore, by themselves do not single out a particular process responsible for the energy transfer to the ICM. One can therefore pose a question as to the nature of the drag force that balances the buoyancy, once the bubble, which is assumed to be essentially massless, reaches its terminal velocity. For instance, it could be viscous or magnetic stresses, turbulence generated in the wake of the bubbles, advection of low-entropy gas to large radii, excitation of sound waves or internal waves. The latter possibility is attractive on several grounds. First, internal waves are trapped in the central region of a cluster, because the Brunt–Väisälä frequency (N) is a decreasing function of radius (Balbus & Soker 1990; Lufkin et al. 1995), implying that the energy will not leak outside the cluster core. Secondly, these waves can travel in the tangential direction (azimuthal) and spread the energy throughout the core. The question then arises: how efficient is the excitation of internal waves? This is the focus of this study. A comprehensive answer to this question requires detailed knowledge of the ICM micro-physics and the internal composition and structure of the bubbles, which is currently missing. Numerical methods face a serious challenge in modelling bubble dynamics and understanding the relevant heating processes. More specifically, ideal hydrodynamic models lead to rapid destruction of rising bubbles. Boundaries of buoyant bubbles are susceptible to the Rayleigh–Taylor (R-T), Kelvin–Helmholz (K-H) and Richtmyer–Meshkov (R-M) instabilities. Therefore, it is hard for the bubbles to survive for longer than ∼ 0.1−0.2 Gyr in hydrodynamic simulations, unless high viscosity and/or magnetic fields suppress the interface instabilities (e.g. Reynolds et al. 2005; Dong & Stone 2009; Bambic et al. 2018). Observations show that, however, some clusters (e.g. Perseus and M87/Virgo) have X-ray cavities with relatively regular shapes even far from the cluster centre (Forman et al. 2007; Fabian et al. 2011). Phenomenologically, this implies that a strong surface tension acts on the bubble surface and keeps the bubble stable, although the detailed physical description of this effective surface tension, presumably magnetic, is an open problem. To circumvent these problems, we make two major simplifying assumptions: (i) the ICM can be described in the framework of ideal hydrodynamics and (ii) the shape of the bubble does not change as the bubble rises in the cluster atmosphere. Thus, we model the bubbles phenomenologically as rigid bodies buoyantly rising in a stratified cluster atmosphere and discuss the nature of the drag force acting on these bubbles. This paper is organised as follows. In Section 3, we describe the models and simulation methods adopted in this work. The main results from our simulations are presented in Section 4, where we study systematically the generation of internal waves by bubbles of different shapes rising in a stratified atmosphere. In Section 5, we discuss the implications of our findings for the terminal velocity of the bubbles in galaxy clusters, velocity diagnostic with future high-energy-resolution missions, and limitations of our method. In Section 6, we summarize our conclusions. 2 PRELIMINARY CONSIDERATIONS In the scenario of a quasi-continuous radio-mode AGN feedback, the generation of X-ray cavities in the centres of galaxy clusters by the SMBH can be generally summarized as a two-stage process. First, a pair of bubbles are blown by bipolar jets, and subsequently expand until the expansion velocity becomes comparable to the velocity of their rise driven by the buoyancy force. At that moment, the expansion is subsonic and the bubble is close to pressure equilibrium with the surrounding ICM. Secondly, the relic bubbles detach from the cluster centre and buoyantly rise upwards. The bubbles finally reach their terminal velocity when drag balances buoyancy force, as shown in Fig. 1. Figure 1. Open in new tabDownload slide Sketch showing a bubble rising in a stratified medium. The bubble rises at the terminal velocity when the buoyancy force is balanced by the drag force. The grey, black, and orange lines show schematically sound waves, turbulence, and internal waves excited by the moving bubble, which can all contribute to the total drag. Figure 1. Open in new tabDownload slide Sketch showing a bubble rising in a stratified medium. The bubble rises at the terminal velocity when the buoyancy force is balanced by the drag force. The grey, black, and orange lines show schematically sound waves, turbulence, and internal waves excited by the moving bubble, which can all contribute to the total drag. Consider a stably stratified atmosphere in hydrostatic equilibrium, which is locally characterized by the Brunt–Väisälä frequency \begin{eqnarray*} N=\sqrt{\frac{g}{\gamma H_{\rm S}}}=\frac{c_{s}}{\gamma \sqrt{H_{\rm S}H_{\rm P}}}, \end{eqnarray*} (1) where g, γ, and |$\displaystyle c_s$| are the gravitational acceleration, gas adiabatic index, and adiabatic sound speed, respectively; HS = |dln S/dr|−1, where S is the gas entropy; HP = |dln P/dr|−1, where P is the gas pressure. Such an atmosphere supports sound and internal waves, the latter with frequencies ω below N, while the former have frequencies above the so-called acoustic cut-off frequency ωa = cs/2Hρ ≳ N, where Hρ = HP is the density scale height in an isothermal atmosphere. Now consider a buoyant bubble, whose mass can be neglected, rising with terminal velocity Uterm in such an atmosphere. Let us first neglect the stratification and the compressibility of the gas and assume that the drag is purely hydrodynamic, i.e. \begin{equation*} F_{\rm buoyancy}(\equiv g\rho V_{\rm b})=F_{\rm drag}\left(\equiv \frac{1}{2}C_{\rm d}\rho A_{\rm b}U^2_{\rm term}\right), \end{equation*} (2) where Vb and Ab are the volume and the horizontal cross-section area of the bubble, respectively; ρ is the ICM density and Cd is the hydrodynamic drag coefficient. From equation (2), the terminal velocity can be simply estimated to be \begin{equation*} U_{\rm term} = \sqrt{\frac{2gV_{\rm b}}{A_{\rm b}C_{\rm d}}}. \end{equation*} (3) Thus, for a sphere, \begin{equation*} U_{\rm term,\,sphere} = c_s\sqrt{\frac{4}{3\gamma C_{\rm d} } \frac{L}{H_{\rm P}}}\approx 1.3 c_s \sqrt{\frac{L}{H_{\rm P}}}, \end{equation*} (4) where for the last equality we used Cd ≃ 0.5, which is the value found for a sphere moving in an incompressible fluid with Reynolds number ∼103 (Roos & Willmarth 1971); L is the diameter of the bubble. We can now estimate the Froude number \begin{equation*} {\rm Fr}=\frac{U_{\rm term,\,sphere}}{NL}=\left( \frac{4\gamma }{3C_{\rm d}} \frac{H_{\rm S}}{L}\right)^{1/2}\approx 2.1 \sqrt{\frac{H_{\rm S}}{L}}. \end{equation*} (5) Thus, for plausible values of L being a fraction of the scale height, the Froude number is larger than 1, implying that internal-wave generation by a sphere is not going to be efficient. This is based on a number of simplifying assumptions, but the qualitative conclusion is that as long as simple hydrodynamic drag Fhydro dominates in Fdrag ( = Fhydro + Fwave + ...), the resulting Uterm is high and the wave drag Fwave is expected to be small. In the limit of Fr ≫ 1, the wave drag ∝ Fr−4 for a sphere in three dimensions (3D; see e.g. Gorodtsov 2007), implying that if one can reduce the terminal velocity, and consequently the Froude number, then the wave drag might become important. Let us keep the assumption that Fhydro is the dominant term, but consider a flattened bubble, whose vertical size is smaller than its horizontal size (see Fig. 1). For instance, consider a disc with diameter L and thickness h, so that its frontal area Ab is |$\pi$|L2/4 and the volume Vb = |$\pi$|L2h/4. Then the expression (3) for the terminal velocity becomes \begin{equation*} U_{\rm term,\,disc} = c_s\sqrt{\frac{2h}{\gamma C_{\rm d}H_{\rm P}}}\approx c_s\sqrt{\frac{h}{H_{\rm P}}}, \end{equation*} (6) where we use Cd ≈ 1.2 for a disc (Roos & Willmarth 1971). Now only the vertical size h of the bubble enters the expression for terminal velocity. The corresponding Froude number is \begin{equation*} {\rm Fr}=\frac{U_{\rm term,\,disc}}{NL}= \sqrt{\frac{2\gamma }{C_{\rm d}}} \sqrt{\frac{H_{\rm S}}{L}} \sqrt{\frac{h}{L}}\approx 1.7 \sqrt{\frac{H_{\rm S}}{L\epsilon _{\rm b}}}, \end{equation*} (7) where εb ≡ L/h is the ratio of the horizontal and vertical dimensions of the bubble. Note that, when evaluating the Froude number, we use the horizontal size L in the denominator, which is appropriate for εb ≳ 1. Comparing expressions (5) and (7), we see that in the latter expression there is an extra factor |$\epsilon _{\rm b}^{-1/2}$| that can make the Froude number smaller. Thus, for a given horizontal size L, the flattened bubble with εb > 1 will have a smaller terminal velocity and a smaller Froude number. The condition that Fr ≲ 1 corresponds to εb ≳ HS/L. For a sufficiently large bubble with large εb, the excitation of internal waves might become important. Of course, this would make the above estimates of the terminal velocity or the Froude number invalid, because they are based on the drag coefficient Cd corresponding to the unstratified case. It is expected that the drag will be higher due to the contribution of internal waves and other effects of stratification (see the experiment in Lofquist & Purtell 1984) and, therefore, the terminal velocity and the Froude number will be even lower accordingly. Thus, this exercise shows that flattened, buoyantly rising bubbles are potentially much more efficient at generating internal waves than spherical ones. We emphasize again that Uterm itself is a function of Fr, and, therefore, equation (7) is expected to overestimate the true Froude number for a buoyantly rising body. It is therefore important to verify the qualitative statements made in this section with numerical simulations, presented in Section 3. Is there any possible connection between the above discussion and the bubbles of relativistic plasma in clusters? A Chandra X-ray image of the Perseus cluster is shown in Fig. 2. The radio bubbles are seen as dark areas (X-ray dim) in this image. The central two bubbles are more or less spherical. A weak shock that encompasses the bubbles indicates that they are still expanding (see e.g. Fabian et al. 2006; Zhuravleva et al. 2016; Tang & Churazov 2017). However, older bubbles, which are presumably rising buoyantly, are clearly aspherical. For example, shown with a white line is the ellipse with the ratio of major-to-minor axis ∼3.6. While this ellipse is a poor match to the shape of the bubble observed at ∼30 kpc to the north–west (NW) from the nucleus, it nevertheless provides a rough estimate of the degree to which bubbles might be flattened. The flattened bubble shape could be seen as a natural consequence of surface tension acting on the bubble surface. The ram pressure gradient of the flow squeezes the bubble along the direction of its motion, but the surface tension prevents the bubble surface from shredding. Putting aside the question of the nature of this surface tension, we note that substituting into equation (7) the parameters of such a bubble, HS ∼ 40 kpc, L ∼ 25 kpc, εb ∼ 3.6, one gets Fr ≲ 1. It is, therefore, plausible that the NW bubble can excite internal waves. Figure 2. Open in new tabDownload slide Chandra 3.5–7.5 keV band image of the Perseus cluster. The bubbles appear as dark (X-ray dim) regions in this image. ‘Active’ bubbles (i.e. still expanding, its radius ∼7 kpc) are marked with black dashed circles. They are surrounded by quasi-spherical weak shocks (radius ∼14 kpc), shown by yellow circles. The outer bubble to the NW from the centre has the ‘horizontal’ and ‘vertical’ (radial) sizes L ∼ 25 kpc and h ∼ 7 kpc, respectively. Thus, its aspect ratio is εb = L/h ∼ 3.6. We argue that for such bubbles the effects of stratification dominate in the total drag force. Figure 2. Open in new tabDownload slide Chandra 3.5–7.5 keV band image of the Perseus cluster. The bubbles appear as dark (X-ray dim) regions in this image. ‘Active’ bubbles (i.e. still expanding, its radius ∼7 kpc) are marked with black dashed circles. They are surrounded by quasi-spherical weak shocks (radius ∼14 kpc), shown by yellow circles. The outer bubble to the NW from the centre has the ‘horizontal’ and ‘vertical’ (radial) sizes L ∼ 25 kpc and h ∼ 7 kpc, respectively. Thus, its aspect ratio is εb = L/h ∼ 3.6. We argue that for such bubbles the effects of stratification dominate in the total drag force. 3 MODEL AND SIMULATION METHOD In the previous section, we argued that the wave drag acting on flattened bubbles could be important, implying that a significant fraction of energy released during their buoyant rise would go into internal waves. In this section, we outline our programme of testing this argument using numerical simulations. These simulations are not intended to reproduce closely the behavior of real bubbles in real clusters, but rather to demonstrate the ability of flattened bubbles to excite internal waves. To do this, we need the modelled bubbles to maintain their integrity as they move. In other words, we need a sufficiently strong surface tension to act on the bubble surface. While magnetic fields have been considered as a possible candidate to provide this surface tension (Kaiser et al. 2005; Diehl et al. 2008; Bambic et al. 2018), we have chosen to model the bubble as a rigid object and to use hydrodynamic simulations to investigate its dynamic behavior in the ICM (see Fig. 1). This choicedramatically simplifies the simulations and allows us to isolate the issue of internal-wave excitation from the much more difficult problem of determining the configuration of the magnetic field inside and outside the bubble. Given the qualitative character of this study, we perform two-dimensional (2D) simulations, which corresponds to studying a flow past an elliptical cylinderin 3D. Rigid bodies moving in a stratified atmosphere have been widely explored in industrial applications both theoretically (e.g. Warren 1960; Voisin 1994; Torres et al. 2000; Hanazaki, Konishi & Okamura 2009) and experimentally (e.g. Mowbray & Rarity 1967). For uniform vertical motion of a body, the Froude number is considered to be the key parameter, characterizing the generation of the internal waves (Voisin 2007). This conclusion is broadly confirmed in our simulations. 3.1 Simulation method and settings We work in a 2D Cartesian coordinate system (x, y), and assume a fixed gravitational potential to model the cluster environment: \begin{equation*} \Phi (x,\ y) = 2V_{\rm c}^{2}\ln (y+R_{\rm c}), \end{equation*} (8) where Vc and Rc are the scaling parameters for the potential and the core radius in units of kpc, respectively. The plane-parallel ICM is initially assumed to be isothermal and in hydrostatic equilibrium. The corresponding gas density1 is \begin{equation*} \rho _{0}(x,\ y) = \rho _{\rm c}\exp \left [-\frac{\Phi (x,\ y)}{{c_{T}}^2}\right], \end{equation*} (9) where ρc and |${c_{T}}\equiv \sqrt{k_{\rm B}T_{\rm 0}/\mu m_{\rm p}}$| are the density scale and the isothermal sound speed of the gas, respectively; T0, kB, μ ( = 0.6), and mp are the initial gas temperature, Boltzmann constant, mean molecular weight per ion, and proton mass, respectively. The ICM is assumed to be an ideal monatomic gas with the adiabatic index γ = 5/3. The dynamic viscosity of the ICM is set to 0.9 g cm−1 s−1, which is smaller than 5 per cent of the Spitzer value (Braginskii 1958; Spitzer 1962). We model a 2D rigid bubble by an ellipse. Two parameters – the horizontal length (L) and the ratio of the horizontal and vertical dimensions (εb) of the bubble – define its size and shape. To obtain full-grown wakes and to investigate the effects of the bubble motion on the far-field regions, we set the computational domain to be x ∈ [ -150, 150 kpc] and y ∈ [0, 400 kpc], and fix the initial bubble position at (xb0, yb0) = (0, 40 kpc). Though in reality the interaction between the jets and the ICM could produce complex gas motions in the central region of the cluster while the bubble is formed, we assume that both the ICM and the bubble are initiallystatic. The force acting on the buoyant bubble is2 \begin{equation*} \mathbf {\boldsymbol{ F}}_{\rm bub} = \oint _{\rm bubble}{\left(P\mathbf {\boldsymbol{ n}}+\hat{\Pi }{\bf \cdot \boldsymbol{ n}}\right){\rm d}\ell }, \end{equation*} (10) which contains the pressure and viscous components (i.e. the integrals of the gas pressure P and of the shear stress |$\hat{\Pi }$| over the bubble surface). Here, |$\boldsymbol{ n}$| is the unit vector normal to the bubble surface. For simplicity, we only allow the bubble to move vertically and do not consider the component of the force acting on the bubble in the horizontal direction (it should be zero for a symmetric flow). The gravitational force on the rigid bubble is also neglected, because the bubble mass density is assumed to be much smaller than that of the surrounding ICM. The main parameters of our simulations are summarized in Table 1. These simulations can be split into three groups. In the first group, the bubble moves under the action of the buoyancy force in a stratified atmosphere. In the second group, the bubble moves with a constant vertical velocity U in the same stratified atmosphere. In the third group, the bubble moves with a constant vertical velocity U in a uniform atmosphere. The first group is the main content of this study, while the second and third groups are used to clarify the effect of stratification and to verify the dependence of the drag force on the bubble’s velocity. We use Vc = 300 km s−1, Rc = 60 kpc, ρc = 6.77 × 10−25 g cm−3, and kBT0= 1 keV for all our simulations except the runs FE4L12highN and GE4L12V80highN. The initial gas density and the corresponding Brunt–Väisälä frequency profiles are shown in Fig. 3. The pressure and entropy scale heights at y = 100 kpc are HP = 140 and HS = 210 kpc, respectively. The runs FE4L12highN and GE4L12V80highN are performed to investigate the effect of a different gravitational field, where we select Vc = 1000 km s−1, Rc = 60 kpc, ρc = 2.03 × 10−21 g cm−3, and kBT0 = 5 keV, which, for y > 50 kpc, provides a profile similar to that of Perseus cluster. Figure 3. Open in new tabDownload slide Initial gas density profile (in red solid line, left axis; see equation 9) used in all our simulations except the runs FE4L12highN and GE4L12V80highN, and the profile of the Brunt–Väisälä frequency N (in blue dashed line, right axis; see equation 1). Figure 3. Open in new tabDownload slide Initial gas density profile (in red solid line, left axis; see equation 9) used in all our simulations except the runs FE4L12highN and GE4L12V80highN, and the profile of the Brunt–Väisälä frequency N (in blue dashed line, right axis; see equation 1). Table 1. Parameters of simulations. ID εba L ( kpc)b U ( km s−1)c Atmosphered Typee Buoyancy-driven bubbles in stratified atmosphere FE1L6 1 6 – Stratified Float FE1L12 1 12 – Stratified Float FE4L6 4 6 – Stratified Float FE4L12 4 12 – Stratified Float FE4L24 4 24 – Stratified Float FE4L48 4 48 – Stratified Float FE4L12highNf 4 12 – Stratified Float FE9L12 9 12 – Stratified Float Bubbles moving with fixed velocity in stratified atmosphere GE1L12V120 1 12 120 Stratified Fix GE4L12V15 4 12 15 Stratified Fix GE4L12V30 4 12 30 Stratified Fix GE4L12V60 4 12 60 Stratified Fix GE4L12V80highNf 4 12 80 Stratified Fix GE9L12V15 9 12 15 Stratified Fix Bubbles moving with fixed velocity in uniform atmosphere NE1L12V120 1 12 120 Uniform Fix NE4L12V30 4 12 30 Uniform Fix NE9L12V15 9 12 15 Uniform Fix ID εba L ( kpc)b U ( km s−1)c Atmosphered Typee Buoyancy-driven bubbles in stratified atmosphere FE1L6 1 6 – Stratified Float FE1L12 1 12 – Stratified Float FE4L6 4 6 – Stratified Float FE4L12 4 12 – Stratified Float FE4L24 4 24 – Stratified Float FE4L48 4 48 – Stratified Float FE4L12highNf 4 12 – Stratified Float FE9L12 9 12 – Stratified Float Bubbles moving with fixed velocity in stratified atmosphere GE1L12V120 1 12 120 Stratified Fix GE4L12V15 4 12 15 Stratified Fix GE4L12V30 4 12 30 Stratified Fix GE4L12V60 4 12 60 Stratified Fix GE4L12V80highNf 4 12 80 Stratified Fix GE9L12V15 9 12 15 Stratified Fix Bubbles moving with fixed velocity in uniform atmosphere NE1L12V120 1 12 120 Uniform Fix NE4L12V30 4 12 30 Uniform Fix NE9L12V15 9 12 15 Uniform Fix aThe ratio of the horizontal and vertical dimensions of the bubble εb( ≡ L/h). bThe horizontal diameter of the bubble L. cThe velocity of the bubble. For the bubbles that buoyantly float upwards (the ‘float’ type), the velocity evolution is shown in Fig. 5. dThe environment where the bubble moves. eThe motion type of the bubbles: buoyantly rising bubbles (the ‘float’ type) and ones that move at fixed velocity (the ‘fix’ type). fThe runs FE4L12highN and GE4L12V80highN used a different gravitational field from other simulations: they had a higher gas temperature and Brunt–Väisälä frequency that are comparable to those of the Perseus cluster core. Open in new tab Table 1. Parameters of simulations. ID εba L ( kpc)b U ( km s−1)c Atmosphered Typee Buoyancy-driven bubbles in stratified atmosphere FE1L6 1 6 – Stratified Float FE1L12 1 12 – Stratified Float FE4L6 4 6 – Stratified Float FE4L12 4 12 – Stratified Float FE4L24 4 24 – Stratified Float FE4L48 4 48 – Stratified Float FE4L12highNf 4 12 – Stratified Float FE9L12 9 12 – Stratified Float Bubbles moving with fixed velocity in stratified atmosphere GE1L12V120 1 12 120 Stratified Fix GE4L12V15 4 12 15 Stratified Fix GE4L12V30 4 12 30 Stratified Fix GE4L12V60 4 12 60 Stratified Fix GE4L12V80highNf 4 12 80 Stratified Fix GE9L12V15 9 12 15 Stratified Fix Bubbles moving with fixed velocity in uniform atmosphere NE1L12V120 1 12 120 Uniform Fix NE4L12V30 4 12 30 Uniform Fix NE9L12V15 9 12 15 Uniform Fix ID εba L ( kpc)b U ( km s−1)c Atmosphered Typee Buoyancy-driven bubbles in stratified atmosphere FE1L6 1 6 – Stratified Float FE1L12 1 12 – Stratified Float FE4L6 4 6 – Stratified Float FE4L12 4 12 – Stratified Float FE4L24 4 24 – Stratified Float FE4L48 4 48 – Stratified Float FE4L12highNf 4 12 – Stratified Float FE9L12 9 12 – Stratified Float Bubbles moving with fixed velocity in stratified atmosphere GE1L12V120 1 12 120 Stratified Fix GE4L12V15 4 12 15 Stratified Fix GE4L12V30 4 12 30 Stratified Fix GE4L12V60 4 12 60 Stratified Fix GE4L12V80highNf 4 12 80 Stratified Fix GE9L12V15 9 12 15 Stratified Fix Bubbles moving with fixed velocity in uniform atmosphere NE1L12V120 1 12 120 Uniform Fix NE4L12V30 4 12 30 Uniform Fix NE9L12V15 9 12 15 Uniform Fix aThe ratio of the horizontal and vertical dimensions of the bubble εb( ≡ L/h). bThe horizontal diameter of the bubble L. cThe velocity of the bubble. For the bubbles that buoyantly float upwards (the ‘float’ type), the velocity evolution is shown in Fig. 5. dThe environment where the bubble moves. eThe motion type of the bubbles: buoyantly rising bubbles (the ‘float’ type) and ones that move at fixed velocity (the ‘fix’ type). fThe runs FE4L12highN and GE4L12V80highN used a different gravitational field from other simulations: they had a higher gas temperature and Brunt–Väisälä frequency that are comparable to those of the Perseus cluster core. Open in new tab We use the OpenFOAM3 code to carry out our simulations. This code was designed for a wide range of applications of computational fluid dynamics, using the finite-volume method. Its flexibility to deal with non-orthogonal and moving meshes and the implementation of the subgrid turbulence models (e.g. Large Eddy Simulation technique, LES) make it highly suitable for our purposes. The effective spatial resolution of the simulations is 0.5 kpc. We have confirmed the convergence of our simulation (i.e. the run FE4L12) by increasing the resolution by a factor of 2. The high-resolution run shows generally the same results as those from the default one. Further OpenFOAM settings adopted in this work are described in detail in Appendix A. 4 RESULTS In this section, we present the main results of our simulations. Section 4.1 describes a rigid bubble’s buoyant rise through the ICM and its impact on the ICM. Internal waves are generated downstream from the bubble. In Section 4.2, we make a comparison between bubbles moving in a stratified and an unstratified atmosphere. The stratification is found to play an important role in determining the drag coefficient. In Section 4.3, we study the dependence of the effective drag coefficient on the bubble’s parameters and on the gravitational field. In Section 4.4, we compare the internal-wave patterns formed in our simulations with those predicted by linear theory. 4.1 Buoyantly rising bubbles with different degrees of flattening Fig. 4 shows the time sequences of the gas entropy map S ( ≡ T/ρ2/3) from the runs FE1L12 (top row), FE4L12 (middle row), and FE9L12 (bottom row). In these simulations, the bubbles, propelled by the buoyancy force, start to rise at t = 0.4 Fig. 4 illustrates the differences in the entropy perturbations induced by bubbles with different degrees of flattening, characterized by the parameter εb. For the circular-bubble case (i.e. εb = 1, see the top row), a pair of vortices appears quickly on the downstream side of the bubble and survives until the end of the simulation. During the rise of the bubble, low-entropy gas is lifted up and forms a long and narrow channel trailing the moving bubble. As a result of the gravitational pull, the gas at the lower end of the channel decelerates after moving a distance ∼10h (see the top panels at t > 1 Gyr) and starts to fall back (it may further oscillate around equilibrium, although these oscillations are not visible in Fig. 4 due to the short duration of the run FE1L12). The key feature of these simulations is that the velocity of the bubble is high and all perturbations are confined to a narrow region in the wake of the bubble. Figure 4. Open in new tabDownload slide Time evolution of the entropy maps for bubbles with different degrees of flattening εb = L/h. Three simulations are shown FE1L12 (εb = 1, top row), FE4L12 (εb = 4, middle row), and FE9L12 (εb = 9, bottom row). The first one has a circular bubble, while in the other two the bubbles are significantly flattened. The interval between each two successive contour levels in all the panels is 3 keV cm−2. The time elapsed since the beginning of the simulation is indicated in each panel. This figure shows that the circular bubble (upper row) has the largest velocity and the perturbations induced by it in the ICM are confined to a narrow region trailing such a bubble. In contrast, flatter bubbles rise slower and cause perturbations over a larger region (see Section 4.1). Figure 4. Open in new tabDownload slide Time evolution of the entropy maps for bubbles with different degrees of flattening εb = L/h. Three simulations are shown FE1L12 (εb = 1, top row), FE4L12 (εb = 4, middle row), and FE9L12 (εb = 9, bottom row). The first one has a circular bubble, while in the other two the bubbles are significantly flattened. The interval between each two successive contour levels in all the panels is 3 keV cm−2. The time elapsed since the beginning of the simulation is indicated in each panel. This figure shows that the circular bubble (upper row) has the largest velocity and the perturbations induced by it in the ICM are confined to a narrow region trailing such a bubble. In contrast, flatter bubbles rise slower and cause perturbations over a larger region (see Section 4.1). Simulations of elliptical bubbles (see the middle and bottom rows in Fig. 4) show a very different picture. These bubbles rise much slower than circular ones (see Fig. 5for these rise velocities), while the vortices and the wake trailing the bubbles are much less prominent; the perturbations in the iso-entropy surfaces now extend in the horizontal direction over much larger distances. In the later stages of evolution (i.e. after the bubble has risen ∼50h), the downstream vortices behind the elliptical bubbles become unstable, separate from the bubble surface, and are shed downwards. We emphasize here that the vortex-shedding instability might be partly suppressed in our simulations, because the bubble is constrained to move vertically. Figure 5. Open in new tabDownload slide Time evolution of the bubble-rise velocity in runs FE1L12 (red line), FE4L12 (blue line), and FE9L12 (green line). The horizontal axis shows the vertical position of the bubble centre yb, which is a monotonic function of time (similarly in Figs 9 and 10). The terminal velocities for the bubbles with εb = 1,  4,  9 are ∼120,  30,  15 km s−1, respectively. As expected, flattened bubbles rise slower (see Section 4.1). Figure 5. Open in new tabDownload slide Time evolution of the bubble-rise velocity in runs FE1L12 (red line), FE4L12 (blue line), and FE9L12 (green line). The horizontal axis shows the vertical position of the bubble centre yb, which is a monotonic function of time (similarly in Figs 9 and 10). The terminal velocities for the bubbles with εb = 1,  4,  9 are ∼120,  30,  15 km s−1, respectively. As expected, flattened bubbles rise slower (see Section 4.1). The evolution of the bubble velocities is shown in Fig. 5. Independently of their shape, the bubbles experience a rapid acceleration in the early stage of the simulations. The bubbles first overshoot their respective terminal velocities, but then are decelerated by the pressure drag5 and settle approximately at the terminal velocity. Fig. 5 shows clearly that flattened bubbles rise with smaller terminal velocities. We note that the terminal velocity is not exactly constant in the simulations, and its evolution depends mildly on the shape of the potential well. In our simulations, bubbles generally require ∼2 Gyr to achieve the terminal velocity. This is related to the value of the relaxation factor β in the numerical scheme used for calculating the bubble acceleration (see Appendix A). The terminal velocities of the bubbles with different degrees of flattening are 120 km s−1 (εb = 1), 30 km s−1 (εb= 4), and 15 km s−1 (εb = 9). The smaller rise velocities of the strongly flattened bubbles in the simulations are not necessarily caused by stratification. Indeed, in these simulations the horizontal size of the bubble was kept constant, while εb increased from 1 to 9, leading to a decrease of the bubble volume (∝1/εb) and, consequently, of the buoyancy force. Since the frontal area of the bubble does not change, the normal hydrodynamic drag (at a given velocity) varies only as long as the drag coefficient varies. In 2D, the difference in the drag coefficient between a ‘sphere’ and a ‘disc’ is less than a factor of 2. Therefore, even in the absence of stratification effects, one can expect that when the driving force (the buoyancy force in our case) decreases by a factor of ∼9, the terminal velocity of a towed body will also decrease. We will elucidate the contribution of the stratification to the drag in the next section. Here, we only note that the presence of internal waves can already be seen in Fig. 4. The simulations FE4L12 and FE9L12 were run for a long enough time (>8 Gyr) to see the long-term signatures of interaction between the bubble and the ambient ICM. One can clearly see the signs of internal waves in the middle and bottom panels of Fig. 4, viz. the wave-like pattern revealed in the iso-entropy surfaces or the collapse of the eddies far from the bubble (see also Figs 7, 8, and 12 for a clearer view of the internal waves). Generation of internal waves by bubbles generally proceeds via two channels: (1) the bubble can directly excite the internal waves (body-driven generation); (2) the bubble first excites eddies within the wake (and uplifts the gas in the wake), which then generate the internal waves (wake-driven generation). Both of these mechanisms presumably appear in our simulations. In general, the relative amounts of body-driven waves and waves generated by the eddies in the wake depend on the Froude number. The body-driven waves are expected to dominate at low Froude numbers, when gravity prevents the formation of the eddies. For example, Brandt & Rottier (2015) experimentally studied the internal-wave field generated by a horizontally towed sphere and found that the former channel was dominant at Fr < 1 (see also Robey1997). The fraction of energy that goes to internal waves via the eddy excitation may depend on the viscosity and diffusivity of the fluid (both physical and numerical), because the excitation of waves by eddies competes with the viscous dissipation and mixing of the gas within each eddy. In other words, one can expect the Reynolds number to have a significant effect on the solution (see e.g. fig. 2 in Torres et al. 2000). To assess the level of dissipation in our simulations, we have calculated the ratio of the kinetic energy present in the simulation box to the total energy deposited by the bubble \begin{equation*} f_{\rm kin} = \frac{\int \int \rho (x,\ y)u^2(x,\ y)/2{\,\rm d}x\,{\rm d}y}{\int F_{\rm buoyancy}{\,\rm d}y}, \end{equation*} (11) where Fbuoyancy = ρ0gVb is the buoyancy force and u(x, y) is the gas velocity field. Note that the total energy associated with gas motions can be higher than just their kinetic energy due to the contribution of the potential energy. For instance, for internal waves, the contribution of potential energy doubles the energy fraction in the waves. The value of fkin is shown (colour-coded) in Fig. 6 as a function of the Froude and the Reynolds Re = UL/ν numbers at the moment when the bubble arrives at y = 250 kpc, where ν is the physical kinematic viscosity. This figure shows that fkin ranges ∼5−35 per cent for the simulation parameter sets that we have covered, being mainly dependent on the Reynolds number. For a bubble with a given shape, the simple scaling discussed in Section 2 predicts that |$U\propto \sqrt{L}$|⁠, i.e. the Reynolds number increases with L as L3/2, while the Froude number decreases as L−1/2. This generally explains those shown in Fig. 6. For smaller bubbles, more energy is dissipated relatively quickly; while for larger bubbles, most of the released energy survives for longer time (especially given that kinetic energy accounts only for a part of the energy). As we show in the next section, our simulations lead to the appearance of a ‘Christmas tree’ pattern formed by the internal waves in the wake of the bubble. According to the chart compiled by Torres et al. (2000), such patterns are characteristic for flows with Fr ∼ 1 and |$\rm Re\sim a\ few\ 100$| in agreement with Fig. 6. Figure 6. Open in new tabDownload slide The kinetic-energy fraction (defined by equation 11) as a function of the Froude and Reynolds numbers from different sets of simulations, at the moment when the bubble arrives at y = 250 kpc. The kinematic viscosity is set to be ν = 1.5 × 1026 cm2 s−1 (for FE4L12highN) and ν = 3 × 1026 cm2 s−1 (for all other runs); the Brunt–Väisälä frequency is taken at the height y= 150 kpc, viz. N = 8.8  rad Gyr−1 (for FE4L12highN) and N = 1.8 rad Gyr−1 (for all other runs). The ellipses illustrate the shape and the relative size of bubbles used in the simulations. This figure shows the parameter region |$(\rm {\rm Fr},\ Re)$| our simulations have covered and suggests that for smaller bubbles the energy dissipation is significant (see Section 4.1). Figure 6. Open in new tabDownload slide The kinetic-energy fraction (defined by equation 11) as a function of the Froude and Reynolds numbers from different sets of simulations, at the moment when the bubble arrives at y = 250 kpc. The kinematic viscosity is set to be ν = 1.5 × 1026 cm2 s−1 (for FE4L12highN) and ν = 3 × 1026 cm2 s−1 (for all other runs); the Brunt–Väisälä frequency is taken at the height y= 150 kpc, viz. N = 8.8  rad Gyr−1 (for FE4L12highN) and N = 1.8 rad Gyr−1 (for all other runs). The ellipses illustrate the shape and the relative size of bubbles used in the simulations. This figure shows the parameter region |$(\rm {\rm Fr},\ Re)$| our simulations have covered and suggests that for smaller bubbles the energy dissipation is significant (see Section 4.1). Torres et al. (2000) also simulated a sphere vertically moving in a uniformly stratified fluid. They found a caudal fluid column formed by the collapse of the vortices behind the body, with maximum velocity in the column up to 10 times that of the sphere (see also Hanazaki et al. 2009). This rear jet-like structure moving in the opposite direction to the sphere leads to a significant increase of the pressure drag acting on the body. In our simulations, we do see a long narrow channel trailing the circular bubble (together with a pair of vortices), which, however, has a velocity comparable to that of the bubble. It might not correspond to the jet-like structure found in Torres et al. (2000), whose formation seems to depend on the thin boundary layers near the body surface, which may not be well resolved in our simulations. 4.2 Comparison of bubbles moving in stratified and uniform atmospheres We have seen that buoyantly rising bubbles, in particular flattened ones, cause perturbations in the ICM outside the wake region. To confirm that these perturbations are indeed internal waves and to clarify the role of these waves in spreading energy over the volume, we make a direct comparison between simulations of bubbles moving with a fixed velocity in stratified and uniform atmospheres. For each value of εb, we use the terminal velocity obtained in Section 4.1 for buoyantly rising bubbles with the corresponding degrees of flattening (see Table 1 for the simulation parameters).6 4.2.1 Energy flow in the far field The specific kinetic energy u(x, y)2/2 of the gas in stratified and unstratified runs for an elliptical bubble with εb = 4 is shown in Fig. 7. The time-averaged potential and kinetic energies of linear waves are equal. Therefore, the specific kinetic energy is a good tracer of the energy flow associated with linear waves in the system. As seen in the top panels, in an unstratified atmosphere, most of the kinetic energy is confined to the vortices accompanying the bubble and the surrounding area. The bubble perturbs its local field and a narrow channel along its path. In the far-field region (say, |x| > 2L), the gas is barely perturbed. In contrast to the uniform-atmosphere run, in the stratified atmosphere, the bubble makes a significant impact on the far-field region (downstream of the bubble, see the bottom panels of Fig. 7). The excited internal waves form a ‘Christmas tree’ pattern. Figure 7. Open in new tabDownload slide Comparison of the specific kinetic energy of the gas in simulations with bubbles moving in unstratified (run NE4L12V30, top panels) and stratified (run GE4L12V30, bottom panels) atmospheres. In the stratified case, internal waves are excited, revealed by a characteristic ‘Christmas tree’ pattern (this pattern is also known as ‘flared skirt’ or ‘herring bone’). A significant amount of the undissipated energy is transported to the far field by the internal waves (see Section 4.2). Figure 7. Open in new tabDownload slide Comparison of the specific kinetic energy of the gas in simulations with bubbles moving in unstratified (run NE4L12V30, top panels) and stratified (run GE4L12V30, bottom panels) atmospheres. In the stratified case, internal waves are excited, revealed by a characteristic ‘Christmas tree’ pattern (this pattern is also known as ‘flared skirt’ or ‘herring bone’). A significant amount of the undissipated energy is transported to the far field by the internal waves (see Section 4.2). Fig. 8 shows the time evolution of 1D horizontal slices of the specific kinetic energy at y = 20, 60, 100 kpc for the run GE4L12V30. The periodic pattern of crests and troughs shows that the excited internal waves propagate in the horizontal direction away from the wake of the bubble. We can directly read the period Tiw and the wavelength (the reciprocal of the horizontal wavenumber) λiw|| of the internal waves from these time-space images (the resulting values are marked in the panels). For example, the pattern seen in the bottom panel (at y = 20 kpc) has the period Tiw ∼ 3−4 Gyr. The corresponding frequency ωiw( ≡ 2|$\pi$|/Tiw) is a factor of ∼0.5 smaller than the Brunt–Väisälä frequency at this height. The wavelength of the internal wave at y = 20 kpc is λiw|| ∼ 50 kpc. At y = 60 and 100 kpc, the period of the waves becomes longer, Tiw ∼  4−6 Gyr. Figure 8. Open in new tabDownload slide Time evolution of 1D slices of the specific kinetic energy of the gas at y = 20 kpc (bottom panel), 60 kpc (middle panel), and 100 kpc (top panel) for the run GE4L12V30. The periodic pattern reveals internal waves propagating outwards (up and down in this plot). The period of internal waves is approximately Tiw ∼  3−4 Gyr at y = 20 kpc and Tiw ∼  4−6 Gyr at both y = 60 and 100 kpc (marked in the panels; see Section 4.2). Figure 8. Open in new tabDownload slide Time evolution of 1D slices of the specific kinetic energy of the gas at y = 20 kpc (bottom panel), 60 kpc (middle panel), and 100 kpc (top panel) for the run GE4L12V30. The periodic pattern reveals internal waves propagating outwards (up and down in this plot). The period of internal waves is approximately Tiw ∼  3−4 Gyr at y = 20 kpc and Tiw ∼  4−6 Gyr at both y = 60 and 100 kpc (marked in the panels; see Section 4.2). Using the kinetic-energy maps, we can calculate the kinetic-energy fraction in the far field, \begin{equation*} f_{\rm far} = \frac{\int _{|x|\gt 2L}\int \rho (x,\ y)u^2(x,\ y){\rm d}x\,{\rm d}y}{\int \int \rho (x,\ y)u^2(x,\ y){\rm d}x\,{\rm d}y}, \end{equation*} (12) to quantify the efficiency of the energy flow away from the bubble’s wake. A homogeneous distribution of energy corresponds to ffar = (xmax − 2L)/xmax ≈ 0.84, where xmax = 150 kpc is the horizontal half size of the simulated box. The value ∼0.8 would imply that the kinetic energy is efficiently transported in the horizontal direction. The time evolution of ffar for the bubbles with different εb is shown in Fig. 9. For the unstratified cases (see the dashed lines), the fractions are lower than 10 per cent throughout the simulation, and show little dependence on the bubble shape. Thus, in the absence of stratification, the energy does not spread away from the wake in the horizontal direction. In contract, the internal waves formed in a stratified atmosphere provide an efficient way for transporting their energy in horizontal direction (see the solid lines in Fig. 9). The evolution of the energy fraction in the stratified runs shows significant dependence on the bubble flattening (viz. Froude number): thus, the kinetic-energy fractions for the bubbles with εb = 1, 4, 9 are approximately 6, 30, 60 per cent, respectively, at the moment when the bubble arrives at y = 300 kpc. Figure 9. Open in new tabDownload slide Time evolution of the far-field kinetic energy fraction ffar (defined by equation 12) for bubbles with different degrees of flattening, moving in stratified (solid lines, runs GE1L12V120, GE4L12V30, GE9L12V15) and unstratified (dashed lines, runs NE1L12V120, NE4L12V30, NE9L12V15) atmospheres. The result from run GE4L12V80highN is also plotted for comparison (blue dotted line). This figure shows that internal waves are efficient in transporting energy from the bubble into the far field (see Section 4.2). Figure 9. Open in new tabDownload slide Time evolution of the far-field kinetic energy fraction ffar (defined by equation 12) for bubbles with different degrees of flattening, moving in stratified (solid lines, runs GE1L12V120, GE4L12V30, GE9L12V15) and unstratified (dashed lines, runs NE1L12V120, NE4L12V30, NE9L12V15) atmospheres. The result from run GE4L12V80highN is also plotted for comparison (blue dotted line). This figure shows that internal waves are efficient in transporting energy from the bubble into the far field (see Section 4.2). Combining these results with the results presented in Fig. 6, we conclude that, according to our simulations, ∼5−30 per cent of the energy released by bubble is propagated into the far field. Given that the quantitative answer depends on the Reynolds number and, also, that our simulations are done in 2D, the conclusion is simply that a significant (order-unity) fraction of the energy can be transferred to internal waves. We plan to explore a wider parameter range |$\rm (Fr,\ Re)$| and the effects of 3D geometry in future studies. 4.2.2 Increase of effective drag coefficient in stratified medium We now proceed with the evaluation of the effective drag coefficient for bubbles moving in a stratified medium. The total force Fbub (see equation 10) acting on the bubble that rises at terminal velocity is zero. It can be decomposed into two components: the drag force (including the pressure and viscous drags) and the buoyancy force: \begin{equation*} F_{\rm bub}=F_{\rm drag}+F_{\rm buoyancy}= 0. \end{equation*} (13) The stratification-dependent drag coefficient is \begin{equation*} C_{\rm d}({\rm Fr}) = \frac{F_{\rm drag}}{\rho _{0}(x_{\rm b},\,y_{\rm b}) A_{\rm b}U^2/2}, \end{equation*} (14) where U is the velocity of the bubble and (xb,  yb) are the coordinates of the bubble’s centre. To evaluate Cd, we can, in view of equation (13), substitute Fdrag in equation (14) with the buoyancy force Fbuoyancy. Contributions to the drag coefficient can come from the excitation of the internal waves, advected low-entropy gas, viscous drag and turbulence generated in the wake of the bubble. We can separate the drag coefficient into two contributions (see Voisin 2007): one without the stratification, Cd(Fr = ∞), and the remaining part that contains the effects of stratification: \begin{equation*} C_{\rm d}({\rm Fr}) = C_{\rm d}({\rm Fr}=\infty ) + \Delta C_{\rm d}({\rm Fr}). \end{equation*} (15) In Fig. 10, we compare the evolution of the drag coefficient for bubbles moving with the same velocities in stratified and unstratified atmospheres. For the unstratified cases, the drag coefficient Cd(Fr = ∞) ranges between ∼0.6 and ∼1.2 (shown with dashed lines), which depends mildly on the bubble shape. The time-averaged drag coefficient of an infinitely long circular cylinder, as reported in the literature (obtained from numerical simulations or measured in experiments; see e.g. Cantwell & Coles 1983; Breuer 1998), is ∼1.2, which is larger than our estimate. This might be due to two reasons: (1) we do not treat the bubble boundary in such details as the studies focused on measuring the drag coefficient accurately, which requires much higher resolution near the bubble surface; (2) the flow passing a blunt body (not a Stokes flow) is usually unsteady: large-scale vortex shedding in our simulations takes place at later times t > 50h/U, so the system is still in the transition phase at the time when we evaluate the drag coefficient. Despite these issues, our simulations suffice for the purposes of this study, given that the deviations from the values reported in the literature are much smaller than the contribution to the drag caused by the stratification. The solid lines in Fig. 10 show the drag coefficient for a bubble moving at constant velocity (approximately equal to the terminal velocity of a buoyant bubble) in a stratified atmosphere. One can clearly see the prominent effect of stratification on the drag coefficient, which has a significant dependence on the bubble shape, viz. ΔCd(Fr) ∼ 0.5,  2.5,  5.0 for the bubbles with εb = 1,  4,  9, respectively. For comparison, we also show the drag coefficients evaluated for the buoyantly rising bubbles (shown with the dotted lines in Fig. 10). Overall, ΔCd(Fr) for buoyantly moving bubbles follows the corresponding values for the cases with constant velocity equal to the mean terminal velocity, but with obvious oscillations, caused by unsteady motion of the buoyantly rising bubbles. Figure 10. Open in new tabDownload slide The drag coefficient Cd for bubbles with different εb moving at constant velocity in stratified (solid lines, runs GE1L12V120, GE4L12V30, GE9L12V15) and unstratified (dashed lines, runs NE1L12V120, NE4L12V30, NE9L12V15) atmospheres. The drag coefficient for buoyantly rising bubbles is also plotted for comparison (dotted lines, runs FE1L12, FE4L12, FE9L12). This figure shows the dramatic effect of stratification on the drag coefficient for strongly flattened bubbles (see Section 4.2). Figure 10. Open in new tabDownload slide The drag coefficient Cd for bubbles with different εb moving at constant velocity in stratified (solid lines, runs GE1L12V120, GE4L12V30, GE9L12V15) and unstratified (dashed lines, runs NE1L12V120, NE4L12V30, NE9L12V15) atmospheres. The drag coefficient for buoyantly rising bubbles is also plotted for comparison (dotted lines, runs FE1L12, FE4L12, FE9L12). This figure shows the dramatic effect of stratification on the drag coefficient for strongly flattened bubbles (see Section 4.2). 4.3 Dependence of effective drag coefficient on bubble parameters and gravitational field Fig. 11shows the dependence of the effective drag coefficient on the bubble parameters (size, shape, and velocity) and on the gravitational field for a bubble moving in a stratified medium. The purpose of this plot is to verify that the values of the drag coefficient estimated in our simulations can be extrapolated to more realistic cases without doing additional simulations (see discussion of the Perseus and M87/Virgo clusters in Section 5). Two trends are clear from Fig. 11. First, fixing the bubble size and shape but varying the velocity shows that the drag coefficient goes down with the bubble velocity. In other words, ΔCd(Fr) in equation (15) is a decreasing function of the Froude number. Secondly, when the bubble rises buoyantly, the effective drag coefficient is a strong function of the bubble shape, but shows only weak dependence on the bubble size and on the degree of stratification, i.e. HS or g. Therefore, we can assume that the effective drag coefficient of a buoyantly rising bubble is only a function of the bubble shape εb. This conclusion is, of course, based on a limited set of simulations, but those are, arguably, sufficient for a qualitative consideration. As a direct inference, the terminal velocity is proportional to the square root of the bubble size, i.e. |$U_{\rm term}\propto \sqrt{h}$|⁠. Therefore, for two bubbles of the same shape εb, the larger one would rise at a higher terminal velocity. Figure 11. Open in new tabDownload slide The drag coefficient Cd from different sets of the simulations (see Table 1). The Brunt–Väisälä frequency is taken at the height y= 150 kpc, viz. N = 8.8 rad Gyr−1 (for FE4L12highN) and N = 1.8 rad Gyr−1 (for all other runs). The ellipses illustrate the shape and the relative size of bubbles used in the simulations. This figure shows how the drag coefficient depends on the bubble’s characteristics (i.e. its shape εb, size L, and velocity U) and on the level of stratification HS and g (see Section 4.3). Figure 11. Open in new tabDownload slide The drag coefficient Cd from different sets of the simulations (see Table 1). The Brunt–Väisälä frequency is taken at the height y= 150 kpc, viz. N = 8.8 rad Gyr−1 (for FE4L12highN) and N = 1.8 rad Gyr−1 (for all other runs). The ellipses illustrate the shape and the relative size of bubbles used in the simulations. This figure shows how the drag coefficient depends on the bubble’s characteristics (i.e. its shape εb, size L, and velocity U) and on the level of stratification HS and g (see Section 4.3). 4.4 Qualitative comparison of simulated pattern of internal waves with linear theory Voisin (1991, 1994) used linear theory to predict analytically the pattern of internal waves generated by a towed 3D body in a uniformly stratified incompressible fluid. In this section, we make a qualitative comparison between his predictions and our simulations. Fig. 12 shows the specific kinetic energy of the gas for simulations with bubbles moving at different velocities (U = 60,  30,  15 km s−1) but with fixed size (L = 12 kpc) and shape (εb = 4). The surfaces of constant phase of a body-driven internal wave estimated from the linear theory for incompressible fluid (see Appendix B) are overplotted as dashed lines. For calculating these surfaces, we used a dipole source to mimic the moving bubble, rather than its specific shape. This approximation is valid when the wavelength of the internal wave is much longer than the body size. As is seen from Fig. 12, the linear model agrees qualitatively with the pattern of internal waves found in our simulations. The latter pattern is more bent and concentrated near the bottom boundary in the simulations due to the reflective boundary condition and the space-dependent Brunt–Väisälä frequency N used in the simulation (see Fig. 3), whereas constant N = 2 was used in the analytic model. In the linear theory, the wavelength of excited internal waves is approximately proportional to the bubble velocity λiw ∝ U/N. Thus, the distribution of troughs and crests is expected to be more concentrated when the bubble velocity is smaller, consistent with the simulation results. In addition, the linear theory provides an explanation of the weak dependence of the gas velocity on the bubble velocity in the far field. In the linear theory, the gas velocity is |$u(x,\ y)\propto Nm_{\rm d}/Ur_{\rm h}^2$|⁠, where md is the dipole moment proportional to the volume and velocity U of the source (see equation B2), rh is the horizontal distance from the source. Therefore, u(x, y) is not a function of the bubble velocity U. Figure 12. Open in new tabDownload slide Surfaces of constant phase in the linear approximation superposed onto the map of the specific kinetic energy in our simulations. The legends are similar to those for Fig. 7, but for a bubble (ε = 4) moving in a stratified atmosphere at velocities U = 60,  30,  15 km s−1. The snapshots are selected from the runs GE4L12V60, GE4L12V30, and GE4L12V15 in such a way as to capture the bubbles when they have arrived at the same height. The dashed lines represent the surfaces of constant phase of the internal wave predicted by the analytical model. This figure shows that the analytical model of Voisin (1991, 1994) qualitatively reproduces the pattern of internal waves observed in our simulations (see Section 4.4). Figure 12. Open in new tabDownload slide Surfaces of constant phase in the linear approximation superposed onto the map of the specific kinetic energy in our simulations. The legends are similar to those for Fig. 7, but for a bubble (ε = 4) moving in a stratified atmosphere at velocities U = 60,  30,  15 km s−1. The snapshots are selected from the runs GE4L12V60, GE4L12V30, and GE4L12V15 in such a way as to capture the bubbles when they have arrived at the same height. The dashed lines represent the surfaces of constant phase of the internal wave predicted by the analytical model. This figure shows that the analytical model of Voisin (1991, 1994) qualitatively reproduces the pattern of internal waves observed in our simulations (see Section 4.4). Obviously, the non-linear features seen in the simulations (e.g. the wake) could not be captured by the linear model. Fig. 12 shows that stronger turbulence is generated in the wake of the bubble when it rises with higher velocity (i.e. larger Fr and |$\rm Re$|⁠). In this case, a larger fraction of the kinetic energy is associated with the wake and the uplifted gas. 5 DISCUSSION By considering a rigid bubble buoyantly rising in a stratified medium, we have found that flattened bubbles can efficiently excite internal waves, which possibly make an important contribution to the total drag that the ICM exerts on X-ray bubbles. These waves tap the potential energy of the rising bubble and transport this energy horizontally (i.e. azimuthally) away from the bubble’s wake. Given that in clusters internal waves are trapped in the radial direction by the decreasing Brunt–Väisälä frequency, the waves must eventually dissipate and heat the ICM. This conclusion depends critically on the ability of the bubbles to maintain their integrity. For instance, in simulations of Reynolds, Balbus & Schekochihin (2015), explosive energy release in localized volumes produced shock-heated bubbles of gas that were rapidly shredded by hydrodynamic instabilities and mixed with the surrounding ICM. No efficient excitation of the internal waves was seen (or expected) and heating took place locally via mixing of the hot gas in the bubbles with the ambient ICM (cf. Hillel & Soker 2014). Guided by the results of the simulations discussed in the previous section, we consider the implications of the buoyantly rising (flattened) bubbles on the heating of the ICM in the context of the radio-mode AGN feedback. 5.1 Constraints on bubble terminal velocity in Perseus cluster and M87 In our simulations, we have found that for bubbles with different degrees of flattening (viz. εb = 1−9), the effective drag coefficients are in the range ofCd ≃ 1−7 (see Sections 4.2 and 4.3). The weak dependence of the drag coefficient on the bubble’s size (as opposed to shape) and on the degree of stratification of the ICM allows us to estimate the bubble’s terminal velocity in galaxy clusters by combining equation (3) and the observed radial profiles of density and pressure. These estimated terminal velocities as functions of radius for the Perseus cluster and M87 are shown in Fig. 13. For these calculations, we used the radial profiles measured via the X-ray observations (see appendix F in Tang & Churazov 2017 and references therein). The shape and height of the bubble are fixed to εb = 4 and Vb/Ab = 5 kpc (motivated by the parameters of the NW bubble in the Perseus cluster shown in Fig. 2), for which Cd = 3. For the Perseus cluster, the estimated terminal velocity is in the range of ∼150−200 km s−1 within r < 400 kpc. This velocity corresponds to ∼0.2cs, a factor of 2−3 slower than that assumed in Churazov et al. (2000, 2001), whose estimates were made for a rising spherical bubble. Figure 13. Open in new tabDownload slide Estimated terminal velocity of a bubble as a function of the cluster radius for Perseus, for M87/Virgo and for the default model used in our simulations (from equation 3), assuming εb = 4 and Vb/Ab = 5 kpc, for which Cd = 3. In the Perseus, the expected terminal velocity of the NW bubble is ∼100−200 km s−1 (see Section 5.1). Figure 13. Open in new tabDownload slide Estimated terminal velocity of a bubble as a function of the cluster radius for Perseus, for M87/Virgo and for the default model used in our simulations (from equation 3), assuming εb = 4 and Vb/Ab = 5 kpc, for which Cd = 3. In the Perseus, the expected terminal velocity of the NW bubble is ∼100−200 km s−1 (see Section 5.1). The morphology of the wake formed behind the rising bubble could provide an independent way of constraining the rise velocity of the buoyant bubbles. In the simulations shown in Fig. 4, the vertical size l of the vortex behind the bubble is ≃  75,  30,  15 kpc for the bubbles with εb = 1,  4,  9, respectively. These simulations suggest that one can roughly estimate l from the condition that U/Nl ∼ 1. Therefore, the rise velocity of the buoyant bubbles in galaxy clusters could be estimated as Nl if we know the size of the vortex. This could be inferred, for instance, from the horseshoe-shaped H α filaments found in the Perseus cluster (Conselice et al. 2001; Fabian et al. 2003; Hatch et al. 2006) trailing the NW bubble, shown in Fig. 2 (see also fig. 3 in Fabian et al. 2003). The size of the region occupied by these filaments is l ≃  20−30 kpc, and therefore, the rise velocity is U ≃  200−300 km s−1. This estimate agrees with the results shown in Fig. 13. Moreover, the Hitomi Collaboration (2018) reported the velocity dispersion |$\sigma _{\rm v}=203^{+16}_{-15}{\rm \,km\,s^{-1}}$| around the region of the NW bubble, while the corresponding bulk velocity was relatively small: |$v_{\rm bulk}=60^{+20}_{-21}{\rm \,km\,s^{-1}}$|⁠. All these results broadly agree with the scenario that the NW bubble rises close to the plane of the sky with velocity ∼200 km s−1. On a more speculative note, it is possible to associate the absorption feature trailing the NW bubble in Perseus cluster with the formation of fluid rear jets in the bubble’s wake, caused by stratification. In simulations and some laboratory experiments (see e.g. Torres et al. 2000), the velocities in the direction opposite to the body motion were up to 10 times the body’s velocity. While formation of such jets may depend on the level of stratification, the shape of the bubble, and the properties of the ICM, it is tempting to consider the possibility that the so-called high-velocity system (Minkowski 1957) might originate from such a mechanism. 5.2 Velocity diagnostic with future high-energy-resolution missions Future missions like ATHENA (e.g. Nandra et al. 2013) and LYNX (e.g. Gaskin, Özel & Vikhlinin2016 ), which offer a combination of high angular and energy resolutions, will dramatically change our ability to probe the velocity field in the cool-core regions. Although a clean direct detection of internal waves generated by individual bubbles in the far field might be challenging due to projection effects and interference of patterns from several bubbles, the velocities in the near field might offer a cleaner test of the effects of stratification. In particular, the simplest test would be done by measuring the magnitude and the spatial extent of the bubble-induced velocities ahead of the bubble. If one neglects stratification, the velocity field ahead of a flattened body would resemble a potential flow affecting a region with size of the order of the transverse size of the body. Strong stratification works towards restoring the isopycnal (or isoentropic) surfaces, thus pushing these surfaces closer to the body. The result is a flow that has a smaller radial extent compared to the unstratified case, as clearly seen from the comparison of the top and bottom panels in Fig. 7. At the same time, the horizontal velocity increases to compensate for the reduced cross-section of the flow. Therefore, for a bubble rising approximately in the plane of the sky, mapping the velocity field ahead of the bubble could allow one to determine both the level of the stratification and the bubble’s rise velocity. Another interesting diagnostic could arise from studying the flow behind the bubble, which is sensitive to the combination of the degree of stratification and effective viscosity (see e.g. Torres et al. 2000). We defer quantitative predictions for these effects to future work. 5.3 Energy flow in cool cores If bubbles of relativistic plasma are capable of maintaining their integrity as they rise buoyantly, then the following simple scenario describes the overall energy flow in cool cores of relaxed clusters. The vertical transport of energy is provided by rising bubbles, which transfer energy to internal waves and turbulent wakes. These internal waves spread the energy over a larger volume of the ICM, principally confined to the cluster core. Since internal waves are trapped in the core of the cluster (Balbus & Soker 1990), they become non-linear and dissipate (e.g. Zhuravleva et al. 2014). Given that bubbles lose their energy after crossing a few pressure scale heights, the efficiency of coupling their energy to the ICM is very high, which is the characteristic feature of the rising-bubble scenario (Churazov et al. 2001, 2002; Begelman 2001). In the cool core of a cluster at any moment, the number of bubbles expected in the flattened-bubble condition is higher than that in the spherical-bubble condition due to a lower terminal velocity. Obviously, in the steady-state scenario, the total dissipation rate does not depend on the rise velocity of the bubbles, but is equal to the rate at which the AGN deposits energy into the bubbles. 5.4 Caveats Our model of a rigid-2D bubble has two obvious limitations. Firstly, because the volume of the bubble does not change with time, only the potential energy, i.e. PV, can be released to the fluid. The bubble’s internal energy (which is 3PV, assuming relativistic plasma inside the bubble) is conserved. However, the expansion of a ‘real’ bubble, which stays in pressure equilibrium with the ambient gas, would release this energy (see e.g. Churazov et al. 2001). This is a two-stage process: the expansion first converts the bubble’s thermal energy into the PdV work, which is then transferred to the fluid via the excitation of waves and wakes. An associated question is the nature and the possible role of effective surface tension needed to prevent the breakup of bubbles in the grid-based numerical schemes. Given the strongly flattened shapes of some of the observed bubbles, we conjecture that the surface energy does not dominate the energy budget. Another limitation of our approach is the 2D geometry adopted in the simulations. From the point of view of the order-of-magnitude estimates made in Section 2, the arguments would remain the same, but the 2D simulations may not capture correctly the properties of the eddies formed in the wake of the bubble or give quantitatively correct energy partition between near and far field. We plan to address this in a future study in 3D. 6 CONCLUSIONS We have shown that buoyant bubbles of relativistic plasma, ubiquitously found in the cores of galaxy clusters, can be efficient emitters of internal waves. A sufficient condition for this is a substantial flattening of the bubbles, so that their vertical (radial) size is a factor of a few times smaller than the azimuthal one. Their terminal velocity scales as the square root of their size, but the effective drag coefficient is larger than for a spherical bubble. For such flattened bubbles, the terminal velocities are small enough so that the Froude number Fr ≲ 1. The stratification effects make the dominant contribution to the total drag force balancing the buoyancy force, and clear signs of internal waves are seen in these simulations. The internal waves propagate horizontally and downwards from the rising bubble, spreading the energy over large volumes of the ICM, where it eventually dissipates. As in the majority of other buoyancy-controlled scenarios, the coupling efficiency of the AGN mechanical power with the ICM is large and the energy does not leak outside cluster cores. When scaled to the conditions of the Perseus cluster and M87, the expected terminal velocity is ∼100–200 km s−1 near the cluster cores (assuming the bubble vertical size of 5 kpc). For the Perseus cluster, the velocities are in broad agreement with direct measurements by the Hitomi satellite. The approach adopted in this study is qualitative. We used 2D simulations and modelled the bubbles as rigid bodies, rather than solving a self-consistent problem, which would require a much more sophisticated description of the stratified atmosphere and of the bubble itself. Nevertheless, we believe that the concept of flattened bubbles correctly captures the most important features of the proposed model, namely increased drag coefficient and, as a consequence, more efficient generation of the internal waves. We defer a study of the 3D geometry and probing broader parameter space in the |$\rm (Fr,\ Re)$| plane to future work. Footnotes 1 Hereafter, the variables with a zero in the subscript refer to the initial conditions. 2 Since we work in 2D geometry, the length of the infinitesimal bubble boundary element dℓ is used in lieu of the surface element of a 3D body. 3 Open Source Field Operation and Manipulation, www.openfoam.org. 4 A relaxation process for bubble acceleration is invoked to reduce the effects of the bubble’s impulsive start (see Appendix A). 5 The viscous drag is much smaller than the pressure drag in our simulations, though it is still included in the calculations. 6 To reduce the effect of an impulsive start, we make the bubbles start with zero velocity and accelerate them artificially to the preset velocity Uterm over ∼0.5 Gyr. 7 This theory was developed for 3D geometry. ACKNOWLEDGEMENTS EC acknowledges support by grant No. 14-22-00271 from the Russian Scientific Foundation. The work of AAS was supported in part by grants from UK STFC (ST/N000919/1) and EPSRC (EP/M022331/1). REFERENCES Balbus S. A. , Soker N. , 1990 , ApJ , 357 , 353 https://doi.org/10.1086/168926 Crossref Search ADS Bambic C. J. , Morsony B. J. , Reynolds C. S. , 2018 , ApJ , 857 , 84 https://doi.org/10.3847/1538-4357/aab558 Crossref Search ADS Begelman M. C. , 2001 , in Laing R. A. , Blundell K. M. , eds, ASP Conf. Ser. Vol. 250, Particles and Fields in Radio Galaxies Conference . Astron. Soc. Pac , San Francisco , p. 443 Google Preview WorldCat COPAC Bîrzan L. , Rafferty D. A. , McNamara B. R. , Wise M. W. , Nulsen P. E. J. , 2004 , ApJ , 607 , 800 https://doi.org/10.1086/383519 Crossref Search ADS Boehringer H. , Voges W. , Fabian A. C. , Edge A. C. , Neumann D. M. , 1993 , MNRAS , 264 , L25 https://doi.org/10.1093/mnras/264.1.L25 Crossref Search ADS Boehringer H. , Nulsen P. E. J. , Braun R. , Fabian A. C. , 1995 , MNRAS , 274 , L67 https://doi.org/10.1093/mnras/274.1.L67 Crossref Search ADS Braginskii S. I. , 1958 , Sov. J. Exp. Theor. Phys. , 6 , 358 Brandt A. , Rottier J. R. , 2015 , J. Fluid Mech. , 769 , 103 https://doi.org/10.1017/jfm.2015.96 Crossref Search ADS Breuer M. , 1998 , Int. Numer. Methods Fluids , 28 , 1281 https://doi.org/10.1002/(SICI)1097-0363(19981215)28:91281::AID-FLD759>3.0.CO;2-# Crossref Search ADS Cantwell B. , Coles D. , 1983 , J. Fluid Mech. , 136 , 321 https://doi.org/10.1017/S0022112083002189 Crossref Search ADS Churazov E. , Forman W. , Jones C. , Böhringer H. , 2000 , A&A , 356 , 788 Churazov E. , Brüggen M. , Kaiser C. R. , Böhringer H. , Forman W. , 2001 , ApJ , 554 , 261 https://doi.org/10.1086/321357 Crossref Search ADS Churazov E. , Sunyaev R. , Forman W. , Böhringer H. , 2002 , MNRAS , 332 , 729 https://doi.org/10.1046/j.1365-8711.2002.05332.x Crossref Search ADS Conselice C. J. , Gallagher J. S. III , Wyse R. F. G. , 2001 , AJ , 122 , 2281 https://doi.org/10.1086/323534 Crossref Search ADS Diehl S. , Li H. , Fryer C. L. , Rafferty D. , 2008 , ApJ , 687 , 173 https://doi.org/10.1086/591310 Crossref Search ADS Dong R. , Stone J. M. , 2009 , ApJ , 704 , 1309 https://doi.org/10.1088/0004-637X/704/2/1309 Crossref Search ADS Fabian A. C. , 2012 , ARA&A , 50 , 455 https://doi.org/10.1146/annurev-astro-081811-125521 Crossref Search ADS Fabian A. C. et al. , 2000 , MNRAS , 318 , L65 https://doi.org/10.1046/j.1365-8711.2000.03904.x Crossref Search ADS Fabian A. C. , Sanders J. S. , Crawford C. S. , Conselice C. J. , Gallagher J. S. , Wyse R. F. G. , 2003 , MNRAS , 344 , L48 https://doi.org/10.1046/j.1365-8711.2003.06856.x Crossref Search ADS Fabian A. C. , Sanders J. S. , Taylor G. B. , Allen S. W. , Crawford C. S. , Johnstone R. M. , Iwasawa K. , 2006 , MNRAS , 366 , 417 https://doi.org/10.1111/j.1365-2966.2005.09896.x Crossref Search ADS Fabian A. C. et al. , 2011 , MNRAS , 418 , 2154 https://doi.org/10.1111/j.1365-2966.2011.19402.x Crossref Search ADS Forman W. et al. , 2007 , ApJ , 665 , 1057 https://doi.org/10.1086/519480 Crossref Search ADS Gaskin J. , Özel F. , Vikhlinin A. , 2016 , Space Telescopes and Instrumentation: Optical, Infrared, and Millimeter Wave Proc. SPIE Conf. Ser. Vol. 9904 . SPIE , Bellingham , p. 99040N Google Preview WorldCat COPAC Gorodtsov V. A. , 2007 , Physics – Doklady , 52 , 165 https://doi.org/10.1134/S1028335807030093 Crossref Search ADS Gull S. F. , Northover K. J. E. , 1973 , Nature , 244 , 80 https://doi.org/10.1038/244080a0 Crossref Search ADS Hanazaki H. , Konishi K. , Okamura T. , 2009 , Phys. Fluids , 21 , 026602 https://doi.org/10.1063/1.3075953 Crossref Search ADS Hatch N. A. , Crawford C. S. , Johnstone R. M. , Fabian A. C. , 2006 , MNRAS , 367 , 433 https://doi.org/10.1111/j.1365-2966.2006.09970.x Crossref Search ADS Hillel S. , Soker N. , 2014 , MNRAS , 445 , 4161 https://doi.org/10.1093/mnras/stu2047 Crossref Search ADS Hitomi Collaboration et al. , 2018 , PASJ , 70 , 9 Hlavacek-Larrondo J. , Fabian A. C. , Edge A. C. , Ebeling H. , Sanders J. S. , Hogan M. T. , Taylor G. B. , 2012 , MNRAS , 421 , 1360 https://doi.org/10.1111/j.1365-2966.2011.20405.x Crossref Search ADS Kaiser C. R. , Pavlovski G. , Pope E. C. D. , Fangohr H. , 2005 , MNRAS , 359 , 493 https://doi.org/10.1111/j.1365-2966.2005.08902.x Crossref Search ADS Lofquist K. E. B. , Purtell L. P. , 1984 , J. Fluid Mech. , 148 , 271 https://doi.org/10.1017/S0022112084002342 Crossref Search ADS Lufkin E. A. , Balbus S. A. , Hawley J. F. , 1995 , ApJ , 446 , 529 https://doi.org/10.1086/175811 Crossref Search ADS McNamara B. R. , Nulsen P. E. J. , 2012 , New J. Phys. , 14 , 055023 https://doi.org/10.1088/1367-2630/14/5/055023 Crossref Search ADS McNamara B. R. et al. , 2000 , ApJ , 534 , L135 https://doi.org/10.1086/312662 Crossref Search ADS Minkowski R. , 1957 , Radio Astron. , 4 , 107 Mowbray D. E. , Rarity B. S. H. , 1967 , J. Fluid Mech. , 28 , 1 https://doi.org/10.1017/S0022112067001867 Crossref Search ADS Nandra K. et al. , 2013 , preprint (arXiv:1306.2307) Reynolds C. S. , McKernan B. , Fabian A. C. , Stone J. M. , Vernaleo J. C. , 2005 , MNRAS , 357 , 242 https://doi.org/10.1111/j.1365-2966.2005.08643.x Crossref Search ADS Reynolds C. S. , Balbus S. A. , Schekochihin A. A. , 2015 , ApJ , 815 , 41 https://doi.org/10.1088/0004-637X/815/1/41 Crossref Search ADS Robey H. F. , 1997 , Phys. Fluids , 9 , 3353 https://doi.org/10.1063/1.869448 Crossref Search ADS Roos F. W. , Willmarth W. W. , 1971 , AIAA J. , 9 , 285 https://doi.org/10.2514/3.6164 Crossref Search ADS Soker N. , 2016 , NewAR , 75 , 1 https://doi.org/10.1016/j.newar.2016.08.002 Crossref Search ADS Spitzer L. , 1962 , Physics of Fully Ionized Gases , 2nd edn. Interscience . New York Google Preview WorldCat COPAC Tang X. , Churazov E. , 2017 , MNRAS , 468 , 3516 https://doi.org/10.1093/mnras/stx590 Crossref Search ADS Torres C. R. , Hanazaki H. , Ochoa J. , Castillo J. , van Woert M. , 2000 , J. Fluid Mech. , 417 , 211 https://doi.org/10.1017/S0022112000001002 Crossref Search ADS Vikhlinin A. A. , Kravtsov A. V. , Markevich M. L. , Sunyaev R. A. , Churazov E. M. , 2014 , Phys. Uspekhi , 57 , 317 https://doi.org/10.3367/UFNe.0184.201404a.0339 Crossref Search ADS Voisin B. , 1991 , J. Fluid Mech. , 231 , 439 https://doi.org/10.1017/S0022112091003464 Crossref Search ADS Voisin B. , 1994 , J. Fluid Mech. , 261 , 333 https://doi.org/10.1017/S0022112094000364 Crossref Search ADS Voisin B. , 2007 , J. Fluid Mech. , 574 , 273 https://doi.org/10.1017/S0022112006004095 Crossref Search ADS Warren F. W. G. , 1960 , J. Fluid Mech. , 7 , 209 https://doi.org/10.1017/S0022112060001444 Crossref Search ADS Zhuravleva I. et al. , 2014 , Nature , 515 , 85 https://doi.org/10.1038/nature13830 Crossref Search ADS PubMed Zhuravleva I. et al. , 2016 , MNRAS , 458 , 2902 https://doi.org/10.1093/mnras/stw520 Crossref Search ADS APPENDIX A: THE OPENFOAM SETTINGS Here, we summarize the detailed OpenFOAM settings (including the solver, the mesh, and the subgrid turbulence model) used in our simulations. We performed our simulations by using the built-in rhoPimpleDyMFoam solver, which was developed for dealing with compressible flows. We modified this solver by adding a static gravitational field (see equation 8). The snappyHexMesh tool was applied to generate the mesh for the computational domain. The bubble was modelled as a moving wall with a no-slip boundary condition. Reflecting boundary condition was adopted at the edges of the domain in both of the x and y directions. At the bottom and top boundaries, the ‘fixedFluxPressure’ condition is used for the pressure to keep the stability of the atmosphere. The LES approach was employed in our simulations. This efficiently reduces the computational cost by subgrid modelling of the small-scale turbulence instead of directly simulating it. We select the one-equation eddy-viscosity model (kEqn) to handle the subgrid turbulence. The turbulent energy and the modified turbulent viscosity were uniformly set to be k = 1011 cm2 s−2 and |$\tilde{\nu }=3\times 10^{26}\,{\rm cm^{2}\,s^{-1}}$| in the initial conditions. We further tested different subgrid turbulence models (e.g. the k-omega-SST-DES model) and model parameters and found that they had little impact on our results. To reduce the effects of an impulsive start of the bubble and improve computational stability, we included a relaxation process for bubble acceleration. This limits the amount of acceleration changing from one time step to the next, i.e. |$a_i = (1-\beta )a_{i-1} + \beta \hat{a}_i$|⁠, where ai − 1 and ai are the bubble acceleration at steps i − 1 and i, respectively; |$\hat{a}_i$| is the expected acceleration obtained fromequation (10) at the step i; and β ( = 10−4) is the relaxation factor. The relaxation affects bubble motion until the bubble’s velocity approaches its terminal velocity; it also filters out fast variations in the bubble velocity. Since we mainly focus on the excitation of internal waves by bubbles moving at the terminal velocity and on time-averaged quantities (e.g. the effective drag coefficient), this artificial relaxation does not affect ourconclusions. APPENDIX B: LINEAR MODEL OF INTERNAL-WAVE GENERATION Here, we summarize the linear theory (Voisin 1991, 1994) applied in Section 4.4 to predict analytically the pattern of internal waves generated by bubbles.7 Assuming a uniformly stratified fluid (i.e. constant Brunt–Väisälä frequency N) and using the Boussinesq approximation, Voisin (1994) found a static solution for the Green’s function of the far-field velocity distribution umono in the case of a uniform vertical motion of a monopole source m0: in the rest frame of the source (see his equation 6.9), \begin{equation*} u_{\rm mono}\propto \frac{Nm_{0}}{Ur_{h}}\cos \Big [{\frac{N}{|U|}}f(r_{h},\ r_{z})\Big ], \end{equation*} (B1) where f(rh, rz) is a function of the horizontal (rh) and vertical (rz) distance from the source. In this study, we consider a dipole source to model the finite volume of the moving bubble: \begin{equation*} u_{\rm dipole}=\frac{m_{\rm d}}{m_{0}}\frac{{\rm d} u_{\rm mono}}{{\rm d}r_{z}}, \end{equation*} (B2) where md ∝ a3U is the dipole moment and a is the effective radius of the sphere modelled by the dipole. This approximation is valid in the far field where the wavelength of the internal wave is much longer than the bubble size L. In equation (B2), md/m0 determines the amplitude of the velocity distribution but has nothing to do with its pattern. The surfaces of the constant phase of udipole are shown in Fig. 12, where N is fixed to 2. © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) TI - Generation of internal waves by buoyant bubbles in galaxy clusters and heating of intracluster medium JF - Monthly Notices of the Royal Astronomical Society DO - 10.1093/mnras/sty1269 DA - 2018-08-21 UR - https://www.deepdyve.com/lp/oxford-university-press/generation-of-internal-waves-by-buoyant-bubbles-in-galaxy-clusters-and-6KQ17BRa4c SP - 4785 VL - 478 IS - 4 DP - DeepDyve ER -