The Sunyaev-Zel'dovich effect of simulated jet-inflated bubbles in clusters

TL;DR

This paper investigates the Sunyaev-Zel'dovich effect of jet-inflated bubbles in galaxy clusters using 3D magnetohydrodynamical simulations, revealing limitations of common models and proposing methods to better infer bubble contents.

Contribution

It introduces a detailed simulation-based analysis of SZ signals from AGN-inflated bubbles and develops an approximate method to determine jet inclination and contents.

Findings

Ellipsoidal bubble models only match inner edges, missing shock cocoon emission.

Kinetic SZ effect becomes significant at small inclinations due to entrainment.

A new method combines jet parameters and SZ effects to infer bubble contents.

Abstract

Feedback by active galactic nuclei (AGNs) is essential for regulating the fast radiative cooling of low-entropy gas at the centers of galaxy clusters and for reducing star formation rates of central ellipticals. The details of self-regulation depend critically on the unknown contents of AGN-inflated bubbles. Observations of the Sunyaev-Zeldovich (SZ) signal of AGN bubbles provide us with the ability to directly measure the lobe electron pressure given a bubble morphology. Here we compute the SZ signal of jet-inflated bubbles in three-dimensional magnetohydrodynamical simulations of the galaxy cluster MS0735.6+7421 with the Arepo code, and compare our synthetic SZ results to inferences obtained with popular modelling approaches. We find that cutting out ellipsoidal bubbles from a double-beta pressure profile only matches the inner bubble edges in the simulations and fails to account for the emission of the shock-enhanced pressure cocoon outside the bubbles. This additional contribution significantly worsens the accuracy of the cut-out method for jets with small inclinations with respect to the line of sight. Also, the kinetic SZ effect of the bubbles, a previously neglected contribution, becomes relevant at these smaller inclinations due to entrainment and mixing of the intracluster medium with low-density jet material. Fortunately, the different signs of the kinetic SZ signal in opposite lobes allow modelling this effect. We present an approximate method to determine the jet inclination, which combines jet power and lifetime estimates, the stand-off distance between jet head and bow shock, and the kinetic SZ effect, thereby helping to correctly infer the bubble contents.

Equations

$$\begin{array}{rl}{n}_{\mathrm{e}}& =0.05{\left[1+{\left(\frac{r}{100\text{kpc}}\right)}^2\right]}^{-4.9}{\text{cm}}^{-3}\\ & +0.01{\left[1+{\left(\frac{r}{400\text{kpc}}\right)}^2\right]}^{-1.6}{\text{cm}}^{-3}.\end{array}$$

$$\begin{array}{rl}\delta i(x)=& g(x)y_{\mathrm{gas}}[1+\delta (x,T_e)]-h(x)w_{\mathrm{gas}}\\ & +[j(x)-i(x)]{\tau }_{\mathrm{rel}}.\end{array}$$

$$I(x)=i_{0}i(x)=i_0\frac{x^3}{\exp (x)-1},$$

$$f_{\mathrm{CRe}}(p,\alpha ,p_1,p_2)=\frac{(\alpha -1)p^{-\alpha }}{p_1^{1-\alpha }-p_2^{1-\alpha }},$$

Full text

1. Introduction

The radiative cooling time of the intra-cluster medium (ICM) in the center of cool core clusters is less than $1~{}\mathrm{Gyr}$. The central AGN provides a powerful heating source that offsets cooling and reduces star formation (McNamara & Nulsen, 2012). AGN jets power lobes, which detach and buoyantly rise in the cluster atmosphere. However, the detailed lobe content and thus, the heating mechanism remains uncertain. While X-ray observations can only provide lower limits to the temperature of lobes (Worrall, 2009), the SZ signal (Sunyaev & Zeldovich, 1972) is directly sensitive to the thermal and non-thermal heat contents. Hence, SZ observations of bubbles have been suggested (Pfrommer et al., 2005) and simulated (Prokhorov et al., 2012; Yang et al., 2018) to understand the physics of the heating mechanism. Recently, Abdulla et al. (2018) observed the cavities of the cluster MS 0735.6+7421 (hereafter MS0735) and found that they have very little SZ-contributing material. This suggests a lobe pressure support of diffuse thermal plasma with temperature in excess of hundreds of keV, or non-thermal relativistic particle populations.

Assuming energy equipartition between relativistic electrons and magnetic fields, radio synchrotron observations also necessitate an additional pressure component (Bîrzan et al., 2008), suggesting relativistic protons as a likely candidate that matches jet morphologies (Croston et al., 2018). As these cosmic ray (CR) protons escape into the ICM, they resonantly excite Alfvén waves. Damping of those waves provides a promising heating scenario (Guo & Oh, 2008; Enßlin et al., 2011; Fujita & Ohira, 2012; Pfrommer, 2013; Jacob & Pfrommer, 2017a, b). Jet-driven bubble simulations in a galaxy cluster (Ruszkowski et al., 2017; Ehlert et al., 2018) demonstrate that streaming CRs from the bubbles provide a sufficient heating rate to halt the cooling catastrophe. Alternative AGN heating mechanisms include mixing of hot-bubble gas with the ICM (Soker, 2016), dissipation of sound waves (Fabian et al., 2017), weak shocks (Li et al., 2017), turbulence (Zhuravleva et al., 2014; Bambic et al., 2018b) and/or gravity waves (Bambic et al., 2018a).

Early simulations of AGN bubbles reproduce the main features of observed X-ray cavities (Churazov et al., 2001; Brüggen et al., 2009). Simulations including magnetic fields (Robinson et al., 2004; Ruszkowski et al., 2007) and/or viscosity (Reynolds et al., 2005; Sijacki & Springel, 2006) stabilize the bubble against developing fluid instabilities. Turbulence and substructure show a significant impact on bubble dynamics (Heinz et al., 2006; O’Neill & Jones, 2010; Mendygral et al., 2012). The addition of CRs in bubbles leads to more elongated bubbles that reside closer to the cluster center, which is favoured by observations (Sijacki et al., 2008; Guo & Mathews, 2011). Recent simulations gained higher resilience against numerical mixing due to sophisticated refinement criteria (Weinberger et al., 2017; Bourne & Sijacki, 2017). This allows for more realistic (higher) density contrasts between simulated bubbles and ICM on long time scales.

For the first time, we use MHD jet simulations to study how the thermal and kinetic SZ signal of dynamical jet-blown bubble simulations compare to the simplified modelling of bubbles as ellipsoids. To facilitate observational comparison, we simulate CR-filled bubbles in a turbulent and magnetized cluster. We pick the observationally favoured, largest observed AGN outbreak in MS0735 to exemplify our analysis (McNamara et al., 2005; Colafrancesco, 2005).

We describe our simulation methods and setup in Section 2. In Section 3, we summarize the characteristics of jet evolution and the details of our SZ modelling. We show the expected signal of simulated bubbles with different fillings in Section 4 and compare the SZ signal from a simulated bubble to a modelled ellipsoidal bubble. In Section 5, we discus our findings regarding the kinetic SZ effect from our simulated bubbles and the influence of jet inclination on the SZ signal. We propose a method of combining simulations and observations to constrain jet inclination enabling more stringent limits on bubble content. We conclude in Section 6.

2. Simulations

To study the SZ signal from an AGN bubble, we simulate a jet, which self-consistently inflates an MS0735-like bubble in a turbulently magnetized cluster atmosphere. We use the same simulation techniques as for the fiducial run in Ehlert et al. (2018) with parameters adopted to the outburst and ICM in MS0735.

The dark matter profile is modelled after MS0735 as a static Navarro-Frenk-White (NFW) profile with $M_{200,\mathrm{c}}=1.5\times 10^{15}~{}\mathrm{M}_{\odot}$, $R_{200,\mathrm{c}}=2.43~{}\mathrm{Mpc}$ and concentration parameter $c_{\mathrm{NFW}}=3.8$ (Gitti et al., 2007). The electron number density follows a double-beta profile fit to MS0735 (Vantyghem et al., 2014) modified to obtain a gas fraction of 16% at $R_{200,\mathrm{c}}$:

[TABLE]

The Gaussian-distributed, turbulent magnetic field is generated in Fourier space and exhibits a Kolmogorov power spectrum on scales larger than the injection scale $k_{\mathrm{inj}}=37.5^{-1}~{}\mathrm{kpc}^{-1}$. On scales $k<k_{\mathrm{inj}}$ the spectrum follows a white noise distribution. The magnetic field is scaled in concentric spherical shells to obtain a predefined average magnetic-to-thermal pressure ratio of $X_{B,\mathrm{ICM}}=0.05$. Multiple nested meshes of magnetic fields with decreasing resolution from the central AGN, respectively, are required for the large range in spatial resolution within the box. Overlapping regions of neighbouring meshes are iteratively smoothed and cleaned off magnetic divergence. The initial Cartesian mesh is relaxed to obtain a honeycomb-like structure, which is more efficient for the unstructured, moving mesh code arepo. When evolved, the magnetic field drives turbulence through tension forces, which gradually decrease the magnetic field strength. Thus, we rescale the magnetic field to the desired $X_{B,\mathrm{ICM}}$ to obtain our initial conditions (see Ehlert et al., 2018, for further details).

We model the jets by injecting kinetic energy in two spherical regions with radius $r_{j}=1.65\ \mathrm{kpc}$ on opposite sides of the centrally placed supermassive black hole (BH) particle. Mass is removed from these injection regions and thermal energy is added from the surroundings to obtain low-density jets ($\rho_{\mathrm{jet}}=10^{-28}\ \mathrm{g}\ \mathrm{cm}^{-3}$, $\rho_{\mathrm{jet}}/\rho_{\mathrm{ICM}}\sim 10^{-4}$) in pressure equilibrium with the ICM (for more details, see Weinberger et al., 2017). Our fiducial run features a jet with power $P_{\mathrm{jet}}=2\times 10^{46}\ \mathrm{erg}~{}\mathrm{s}^{-1}$ and lifetime $\tau_{\mathrm{jet}}=150~{}\mathrm{Myr}$, amounting to an injected energy $E_{\mathrm{jet}}=P_{\mathrm{jet}}\tau_{\mathrm{jet}}\approx 10^{62}~{}\mathrm{erg}$ in MS0735 (Vantyghem et al., 2014). In addition, we inject a helical magnetic field in the jet region with magnetic-to-thermal pressure ratio $X_{B,\mathrm{jet}}=0.1$. Our lobes are defined via an advective scalar $X_{\mathrm{jet}}$, which is set to unity in the jet injection region. We define our lobes as the union of all cells with $X_{\mathrm{jet}}>10^{-3}$. During jet activity CRs are accelerated in the lobes by converting thermal energy to CR energy whenever the CR-to-thermal pressure $X_{\mathrm{cr}}$ ratio falls below a predefined value $X_{\mathrm{cr}}<X_{\mathrm{cr,acc}}=1$. We explicitly isolate our jet injection region magnetically to inhibit unphysical CR diffusion.

In addition to advection, CRs are expected to scatter on self-excited Alfvén waves in galaxy clusters (Kulsrud & Pearce, 1969; Enßlin et al., 2011). The low efficiency of Alfvén wave damping in the ICM limits the CRs to stream down their pressure gradient $\bm{\nabla}{P}_{\mathrm{cr}}$ along magnetic field lines at the Alfvén speed $\bm{v}_{\mathrm{A}}$ (Wiener et al., 2013). The damping of Alfvén waves effectively transfers CR to thermal energy, giving rise to the so-called Alfvén heating with a power \mathcal{H}_{\mathrm{cr}}=\left|\bm{v}_{\mathrm{A}}\mathchoice{\mskip 3.0mu\lower 0.86108pt\hbox{\scalebox{1.5}{\cdot}}\mskip 3.0mu}{\mskip 3.0mu\lower 0.86108pt\hbox{\scalebox{1.5}{\cdot}}\mskip 3.0mu}{\lower 1.29167pt\hbox{\scalebox{1.2}{\cdot}}}{\lower 1.29167pt\hbox{\scalebox{1.2}{\cdot}}}\bm{\nabla}{P}_{\mathrm{cr}}\right|.

The equations of MHD are discretized on a moving-mesh and evolved with second-order accuracy using the massively parallel arepo code (Springel, 2010; Pakmor et al., 2016b). Cosmic rays are treated as a second fluid including hadronic and Coulomb losses (Pfrommer et al., 2017). These losses are small in comparison to Alfvénic losses, which we include here. In combination with anisotropic diffusion (Pakmor et al., 2016a) (with a parallel diffusion coefficient $\kappa_{\parallel}=10^{29}\ \mathrm{cm}^{2}\ \mathrm{s}^{-1}$) this is used to emulate CR streaming. We employ mass-based refinement with target mass $m_{\mathrm{target}}=1.5\times 10^{6}\ M_{\odot}$. In addition, we impose a refinement criterion based on the density gradient, jet scalar and cell volume difference as in Weinberger et al. (2017) to maintain the large density contrast at the jet-ICM interface with target volume $V_{\mathrm{target}}$ ($V^{1/3}_{\mathrm{target}}=405\ \mathrm{pc}$).

3. Evolving SZ signal from simulated bubbles

Free electrons in the ICM Compton up-scatter cosmic microwave background (CMB) photons. We follow the procedure described in Pfrommer et al. (2005) to compute the resulting SZ signal of our ICM and bubbles. The relative change $\delta i(x)$ in the flux density as a function of dimensionless frequency $x=h\nu/(kT_{\mathrm{CMB}})$ is given by

[TABLE]

with the Planckian distribution of the CMB given by

[TABLE]

where $i_{0}=2(kT_{\mathrm{CMB}})^{3}/(hc)^{2}$ and $T_{\mathrm{CMB}}=2.725\ \mathrm{K}$. The first term in Equation (2) describes the thermal SZ effect that is proportional to the integrated thermal pressure of the ICM along the line of sight $y_{\mathrm{gas}}\propto\int dl\ n_{\mathrm{e,gas}}kT_{\mathrm{e}}$. Relativistic corrections $\delta(x,T_{\mathrm{e}})$ become relevant at high temperatures $kT_{e}\gtrsim 5\ \mathrm{keV}$ (Mroczkowski et al., 2018). Throughout the analysis, we include relativistic corrections from Itoh et al. (1998). The second term corresponds to the kinetic SZ effect due to gas motion relative to the Hubble flow $w_{\mathrm{gas}}\propto\int dl\ n_{\mathrm{e,gas}}v_{\mathrm{gas}}$, where $v_{\mathrm{gas}}<0$ if the gas is approaching the observer. The last term describes the relativistic SZ effect. For fully relativistic fillings the distribution function of the electrons determines the form of $j(x)$ and $\tau_{\mathrm{rel}}\propto\int dl\ n_{\mathrm{e,rel}}$. Only the terms $g(x)$, $h(x)$ and $i(x)$ depend on observing frequency.

Throughout the Letter, we focus on three different bubble fillings: 1.) a thermal electron distribution with $kT_{\mathrm{e}}=10\ \mathrm{keV}$, 2.) a single-temperature, relativistic Maxwellian with $kT_{\mathrm{e}}=500\ \mathrm{keV}$ and 3.) a single-temperature, relativistic electron population that follows a power-law in normalized momentum-space $p=\gamma_{\mathrm{e}}\beta_{\mathrm{e}}$ defined by

[TABLE]

where $\alpha=2$, $p_{1}=1$ and $p_{2}=10^{3}$. For all three models, we recompute the electron density of lobe cells while keeping the total (simulated) pressure $P_{\mathrm{tot}}$.

In Fig. 1, the initial jet inflates lobes which, after jet shut-down, rise buoyantly in the cluster atmosphere. The high jet power leads to jet velocities approaching the speed of light, where our non-relativistic treatment degrades in accuracy. While this influences the details of the kinematics and shock dissipation (Perucho et al., 2017), our results on bubble morphology and SZ signal are expected to be robust. The jet initially drives a shock wave into the ICM. Generally, the Mach number in the jet direction exceeds that perpendicular to the jet, thereby creating an ellipsoidal shock. The trailing contact discontinuity is clearly visible in the electron number density $n_{e}$ and temperature $T$ maps. The dimensions of the bubbles and the morphology of the contact discontinuity in the X-ray brightness $l_{X}$ are in good agreement with observations of MS0735 (Vantyghem et al., 2014).

Assuming a relativistic power-law distribution of electrons in the bubbles implies a larger contrast of the SZ signal $\delta i_{\nu_{0}}$. Note, a hypothetical thermal filling of the under-dense bubble with densities as shown in Fig. 1 implies relativistic temperatures. The high density contrast between bubble and ICM makes the bubble susceptible to the Rayleigh-Taylor instability. Dense gas streams into the lower part of the bubble, generates turbulence and mixes with the bubble gas. This process progressively dilutes the bubble material and transports it into the ICM until complete disruption.

4. Probing relativistic bubble fillings

The observed decrement in the SZ signal due to different bubble fillings can help to unravel its constituents. In the left panel of Fig. 2, we show the SZ signal across the jet axis for bubbles filled with thermal gas ($kT_{e}=10\ \mathrm{keV}$) or relativistic gas, which either follows a thermal distribution ($kT_{e}=500\ \mathrm{keV}$) or a power-law distribution (“relativistic power-law”). The larger SZ contrast for relativistic fillings becomes apparent.

Abdulla et al. (2018) compared cut-out ellipsoids from a fitted, smooth background cluster with their SZ observations to constrain the bubble filling. Rather than reproducing their analysis, we focus on the consequences for the SZ signal of modelling bubbles as ideal ellipsoids including variations in inclination.

On the right-hand panel in Fig. 2, we contrast the expected SZ signal from our simulated bubble to an ellipsoidal cut-out in our simulation, assuming the unperturbed initial profile. Throughout this Letter, we determine the dimensions of the cut-out ellipsoid and its position from unsharped masked X-ray maps of our simulated bubbles. We assume that the depth of the modelled bubble corresponds to its width. On the right-hand panel in Fig. 2, both bubbles are viewed perpendicular to the jet axis and contain a relativistic power-law filling. The profiles show good agreement at the inner rim of the bubble and start to deviate towards the bubble edges, where we see an enhanced SZ signal in the simulations out to the ellipsoidal shock. This feature corresponds to the increased pressure in the shocked cocoon powered by the outburst and is neglected in the simple cut-out model. This highlights the importance of including the bow shock in the model.

5. Degeneracies with jet inclination

In Fig. 3, we study the influence of inclination on the SZ signal. An angle of $\theta=90$ corresponds to a jet axis that is perpendicular to the line of sight. We compare the amplitude of the kinetic SZ effect $w$ (top row of Fig. 3) to the thermal SZ effect $y$ (middle row of Fig. 3). From left to right, we use later times for larger inclinations to keep the projected distance from the BH to the lower bubble edge approximately constant. We see a strong increase in signal strength of the kinetic SZ effect for larger inclinations. In the wake of the bubble, highly turbulent, dense gas pushes into the region previously occupied by the bubble. This bipolar structure develops significant velocities at high densities that contribute up to 10% of the thermal effect to the total SZ signal at $\theta=33^{\circ}$. Future high sensitivity observations of inclined jets should be able to measure this effect. The expected opposite sign of the signal in the two bubbles aids in modelling this effect.

The bottom panels of Fig. 3 show the expected SZ signal at $30\ \mathrm{GHz}$ for the differently tilted simulations and the cut-out model, for which we always assume an inclination of $\theta=90^{\circ}$. In both models, we adopt a power-law distribution of relativistic electrons in the bubbles. The model reproduces the SZ signal at the inner bubble edges well for large inclinations while the agreement becomes worse for smaller inclinations. When viewed face-on ($\theta=90^{\circ}$), the inclusion of the bow shock in the model is critical to correctly reproduce the SZ signal from the simulation. The bow shock loses its momentum and covers a larger volume as a function of time such that the signal from the shocked region flattens for later times (panels towards the right).

In addition, the SZ contrast due to the relativistically filled bubbles disappears almost entirely for smaller inclinations. Here, lines of sight pass through significantly larger portions of the shocked gas region. This is especially true for the inner part of the bubble where lines of sight intersect the dense central cool core. The effect is amplified for older bubbles through the advanced state of entrainment and mixing of ICM material with the relativistic bubble contents due to the Rayleigh-Taylor instability (see Section 3).

Contamination by the signal from the cluster core appears like a significant source of uncertainty for the analysis of observations as $\theta$ is unconstrained a priori. However, X-ray observations can inform simulations about jet power and bubble/shock ages (e.g., Diehl et al., 2008). Our simulation shows that the stand-off distance between bow shock and jet/bubble head is increasing with time. In particular, entraining the ambient ICM towards the end of the jet lifetime slows the jet down, see Fig. 4. Varying the viewing angle of the jet axis then helps to disentangle projection effects and potentially constrain inclination.

For large inclinations, the projected stand-off distance decreases for decreasing inclinations. However, at small inclinations the outer edge of the ellipsoidally shocked region enlarges the projected bow shock region beyond the projected distance of the upper tip of the bow shock. The ellipticity of the bow shock itself may aid in providing a coarse estimate of the inclination. Even though the details may depend on the concentration of the NFW potential, a spheroidal bow shock should generally favour small angles of inclination. Note that these features have a weak resolution dependence in our simulations, which is studied in Weinberger et al. (2017) and Ehlert et al. (2018). With high-sensitivity observations the kinetic SZ effect can be identified through its bimodality and used as an additional constraint for the suggested method.

6. Conclusion

We show that three-dimensional MHD simulations are instrumental to carefully model the SZ signal of jet-inflated bubbles and conclude:

• •

Relativistic bubble fillings imply a large SZ contrast, which is observable in high-resolution SZ observations.

• •

The SZ profiles of simulations and the cut-out model show good agreement at the inner rim of the bubble and start to deviate towards the bubble edge because the simplified model fails to account for the shock-enhanced pressure cocoon outside the bubbles.

• •

The match between simulations and model becomes worse when considering small inclinations between line-of-sight and jet axis ($\theta\lesssim 45^{\circ}$). This geometry probes a larger fraction of the shocked ICM, which leads to an increase of thermal SZ signal also towards the inner bubble rim region. Additionally, the light material of the bubble is progressively mixed with denser ICM due to the Rayleigh-Taylor instability, making the bubble indistinguishable from the surrounding ICM. Thus, the SZ signal from older bubbles is reduced in comparison to modelling them with cut-out ellipsoids.

• •

At small inclinations, the kinetic SZ effect reaches up to 10% of the thermal SZ effect. The wake of the bubble causes dense ICM to enter the bubble from below, causing a comparably large kinetic SZ signal. The opposing signs of the signal of inclined bipolar outflows are a smoking-gun signature for identifying this kinetic SZ signal.

• •

We propose to constrain the inclination with the stand-off distance between shock and jet/bubble, the elliptical appearance of the bow shock, and (if available) the amplitude of the kinetic SZ effect. To this end, a combination of high-resolution X-ray and SZ observations and full MHD simulations are crucial to break degeneracies due to projection effects.

This Letter opens up the possibility to understand biases associated with simplified SZ modeling of AGN bubbles and to finally constrain their contents so that we can observationally identify the physical processes underlying AGN feedback.

Acknowledgements

We acknowledge support by the European Research Council under ERC-CoG grant CRAGSMAN-646955.