© L. Di Mascolo et al. 2019

1. Introduction

Mergers play a crucial role in the formation of galaxy clusters, which are situated at intersections of the Cosmic Web. These spectacular events can have a profound impact on the intracluster medium (ICM) and the galaxies within these environments (see, e.g., Kravtsov & Borgani 2012). Mergers provide large-scale astrophysical laboratories for plasmas where the mean free path can be substantial (see, e.g., Markevitch & Vikhlinin 2007) and for measuring the self-interaction cross-section of dark matter (Markevitch et al. 2004; Randall et al. 2008; Wittman et al. 2018; Tulin & Yu 2018). Furthermore, the very existence of dark matter was conclusively demonstrated through the merging cluster 1E0657–56, or “Bullet Cluster”, which exhibits spatial offsets between its baryonic and total mass peaks in the X-ray and gravitational lensing maps (Clowe et al. 2006; Bradač et al. 2006; Paraficz et al. 2016).

Key to identifying merging clusters is the detection of shocks in the ICM. A “textbook example of a bow shock” is observed in the X-ray image of the Bullet Cluster (Markevitch et al. 2002). Using 500 ks of Chandra X-ray data, Markevitch (2006) reported a Mach number $\mathcal{M}=3.0\pm 0.4$ for the western, most-prominent shock in the Bullet Cluster, an estimate largely determined by the density jump conditions. We also note that Shimwell et al. (2015) revealed a second shock, on the eastern (opposite) side of the cluster, which we do not consider here.

Here we present deep, continuum Atacama Large Millimeter/submillimeter Array (ALMA) observations, sensitive to the thermal Sunyaev-Zeldovich (tSZ; Sunyaev & Zeldovich 1972) effect, of the main shock in the Bullet Cluster. These observations include data from both the 12 m array (hereafter “ALMA”) and 7 m Atacama Compact (Morita) Array (ACA). As the tSZ effect is linearly sensitive to the line-of-sight integral of the electron thermal pressure (see, e.g., Carlstrom et al. 2002; Mroczkowski et al. 2019a, for reviews), these observations complement the X-ray constraints on plasma density and, less accurately, electron temperature, yielding a ground-based, mm-wave view of the shock properties.

All the results presented in this work have been derived assuming a flat ΛCDM cosmology with Ω_m = 0.30, Ω_Λ = 0.70, and H₀ = 70.0 km s⁻¹ Mpc⁻¹. At the redshift of the Bullet Cluster (z = 0.296), 1″ corresponds to a physical scale of 4.41 kpc. Unless stated differently, all reported best-fitting parameters and their respective uncertainties are obtained from the 50th, 16th, and 84th percentiles of the marginalized posterior distributions, corresponding to the 68% credibility interval of the distribution¹.

2. Data and analysis overview

2.1. ALMA/ACA observations

As part of ALMA Cycle 2 operations ALMA and the ACA observed the Bullet Cluster for a total of 3.1 and 5.9 h integration time respectively in Band 3 (project code 2013.1.00760.S). These wideband observations span the frequency range 84−100 GHz in four 2 GHz-wide spectral windows, centered at approximately 85, 87, 97, and 99 GHz. Our strategy employed a single, deep observation centered approximately on the nose of the shock front, as inferred from the X-ray observations (Clowe et al. 2006; Markevitch 2006). The ALMA and ACA observations were respectively performed in 4 and 11 separate executions spanning 2014, obtaining root-mean-square noise levels of approximately 5 μJy and 45 μJy respectively (as measured in naturally-weighted imaging), and a synthesized beam with a main lobe of 4.01″ × 3.07″ FWHM (PA 81°).

We re-reduced the data using the ALMA pipeline (Shinnaga et al. 2015; Humphreys et al. 2016) in CASA 4.7 (McMullin et al. 2007), producing results consistent with the previous calibration using the script provided on data delivery. Our re-reduction provided a cross-check of the earlier reduction, and was necessary due to backwards-compatibility issues and bug fixes in subsequent CASA releases. The data were calibrated using the default calibration strategy of the ALMA observatory, which has nominal uncertainty ≤5%. However, since the flux calibrators, which included quasars, differed for each execution, we performed a manual cross-check of the values for the designated flux calibrators as well as the phase calibrators, finding they were consistent for the dates spanned by the observations.

An interferometer behaves as a spatial filter, sampling the sky only in Fourier modes corresponding to projected baselines in the array (see, e.g., Condon & Ransom 2016; Di Mascolo et al. 2019). This provides clean imaging free from atmospheric structure, but also leads to two major complexities: incomplete sampling in Fourier space even for the modes accessible to the array, and the lack of recovery of angular scales larger than those corresponding to the shortest projected distances between array elements (i.e., the “missing flux issue”). Based on the uv-space coverage of the ALMA and ACA data presented here, the largest recoverable scales are respectively ∼40″ and ∼55″. As detailed in the next section, we choose to forward-model the observed SZ signal using X-ray-motivated priors to address such issues. Additionally, to avoid known deconvolution biases intrinsic to the CLEAN algorithm (Högbom 1974; Thompson et al. 2015), we perform our analysis directly in visibility (uv) space. We extend the interferometric SZ analysis techniques presented in the appendix of Di Mascolo et al. (2019) to allow fitting pressure discontinuities due to shocks. Our approach builds upon the work of Basu et al. (2016), but incorporates several advances in the parameter-space sampling technique as well as more sophisticated and flexible models allowed by the deeper X-ray and SZ observations. In brief, we build an image-space model of the SZ signal by integrating numerically the three-dimensional pressure distribution model, and applying the proper SZ frequency scaling. The dependence of the SZ signal on the electron temperature is taken into account when modelling the SZ spectrum (Itoh & Nozawa 2004; Chluba et al. 2012, and end of Sect. 2.2.2). The pixel scale is chosen to fulfill the Nyquist sampling criterion for the smallest scales probed. The SZ model image is then Fourier transformed and sampled to the position of the sparse interferometric data. The resulting synthetic visibilities are then employed in combination with the observed ones to evaluate the likelihood at each step of the Bayesian inference procedure.

However, we choose not to model the raw post-calibration data, instead binning the data in each spectral window following the optimal averaging scheme described in Hobson et al. (1995). This is crucial for gaining a significant reduction in data volume and hence computational time.

2.2. Sunyaev-Zeldovich model

A summary of the model priors introduced in this section can be found in Table 1. We test for biases in the parameter reconstruction arising from our specific choice for the distribution of priors by performing a prior-only run, which is done by setting the likelihood to a constant value regardless of model fit (see, e.g., Di Mascolo et al. 2019). As expected, the result of this test simply returns the input distribution of priors.

2.2.1. Shock front

The common approach employed in the study of X-ray observations of shock fronts consists in describing them as spherical sectors within a specific region of the cluster image. This takes advantage of the image-space nature of the X-ray data to select a spatial region narrow enough to allow one to locally approximate the shock front as spherical. However, among the complexities of studying interferometric data is the difficulty of applying any spatial masking. This would entail convolution of the visibilities, inducing a non-trivial correlation between them. To avoid this, a complete two-dimensional model of the observed field is then required.

In order to allow more freedom in the description of the shock front than in the case of a spherical model, we describe the shock front as an axially-symmetric hyperbolic surface (see the dashed line in Fig. 1), with central axis coincident with the direction of the merger and lying in the plane of the sky. Since the interferometric data alone cannot constrain the line-of-sight distribution of pressure, we consider the curvature of the front to be symmetric with respect to the line-of-sight and the plane-of-sky direction. Although this is likely a reasonable assumption, any deviations from cylindrical symmetry may introduce non-negligible systematic errors into our results. In particular, the derived downstream pressure P_e, ds will be related to the true pressure $P_{\text{e},\text{ds}}^{\text{true}}$ as $P_{\mathrm{e},\mathrm{ds}}\approx P_{\mathrm{e},\mathrm{ds}}^{\mathrm{true}}({\ell }_{\mathrm{LoS}}^{\mathrm{true}}/{\ell }_{\mathrm{LoS}})$, where ℓ_LoS and ${\ell }_{\mathrm{LoS}}^{\mathrm{true}}$ are respectively the assumed and true line-of-sight extents of the shock front. We present specific estimates for the ratio ${\ell }_{\mathrm{LoS}}^{\mathrm{true}}/{\ell }_{\mathrm{LoS}}$ when discussing our results in the next sections. It is worth noting that, as already discussed by Wang et al. (2018), this is expected to scale approximately with the square root of the curvature radius of the shock surface. For line-of-sight and plane-of-sky curvatures radii r_LoS and r_PoS, we would then get ${\ell }_{\mathrm{LoS}}^{\mathrm{true}}/{\ell }_{\mathrm{LoS}}\approx \sqrt{r_{\mathrm{PoS}}/r_{\mathrm{LoS}}}$. We refer to Appendix A for a discussion.

A preliminary attempt to determine the shape of the shock front using SZ data alone shows that the parameters defining its geometry are heavily degenerate, and the small extent of the ALMA+ACA field of view does not allow for meaningful constraints. We would like to note this is a consequence of the sole parametrization of the shock geometry, since ALMA+ACA has proven to be able to identify edge positions with a beam-scale precision (Basu et al. 2016). We therefore derive a description of the morphology of the bow shock by finding the hyperbola that best describes the discontinuity observed in the Chandra surface brightness map. Analogous to Ueda et al. (2017), we find the best-matching shock geometry by minimizing the variance of the X-ray image within a defined region. The model is assumed to be simply given by a step function in which the discontinuity has a hyperbolic shape. The values of the function inside and outside the front itself are set equal to the mean photon counts in the respective regions of the X-ray image. In order to gain better leverage on the azimuthal geometry of the shock front, we further split the region in several angular sectors (see Fig. 1). The resulting maximum-a-posterior model for the hyperbolic surface is then employed for describing the profile of the shock front in all the following analyses. To account for a possible mismatch in the shock coordinates (e.g. due to astrometry errors) from the Chandra modelling with respect to ALMA+ACA, we allow for some additional freedom in the nose coordinates (RA, Dec)_shock and axis orientation θ. Specifically, we assign each one priors based on the respective marginalized posteriors derived in the X-ray-matching step described above. Any asymmetry in the recovered parameter uncertainties is modelled by means of split-normal distributions (Wallis 2014).

Fig. 1. Cut-out of the 0.5−2.0 keV Chandra X-ray surface brightness map of the Bullet Cluster. The solid wedges represent the sectors employed to derive the hyperbolic shape best-matching the shock front geometry. The maximum-a-posteriori model is shown as a dashed line, while the dotted contours indicate the corresponding 95% credible interval. Further, the dashed and solid yellow circles denote respectively the ACA and ALMA full-width-at-half-maximum fields of view. For reference, the upstream and downstream gas lie respectively west (right) and east (left) of the shock front.

Apart from its morphology, the main parameter defining the shock is the ratio of downstream to upstream pressure at the jump itself. In practice, for a fixed line-of-sight geometry, the available ALMA+ACA data are mainly sensitive to the absolute difference of the downstream and upstream electron pressures near the nose of the shock, i.e., ΔP_e = P_e, ds − P_e, us; the SZ signal associated with the large-scale distribution of the gas is effectively filtered out (see, e.g., Basu et al. 2016, and discussion in Sect. 3.2). Thus, the modelling of the ALMA signal remains only weakly sensitive to the assumed large-scale model. The immediate downside is the relative pressure jump at the shock, x_P = P_e,ds/P_e,us = 1 + ΔP_e/P_e,us, which serves as a proxy for ℳ, is poorly constrained by the interferometric SZ data alone. In fact, due to the lack of information on the pressure normalization, the marginalized posterior distribution of the Mach number inferred when performing an SZ-only analysis are found to entirely span the corresponding prior interval. To get a meaningful measure of the pressure jump from the ALMA+ACA data, we therefore employ an X-ray-informed analysis of the ALMA+ACA SZ observations as in Sect. 3.3 of Basu et al. (2016), and set the upstream electron pressure P_e, us to the value derived by modelling the Chandra data in a narrow sector centered on the shock nose (Markevitch 2006).

We do not include any model components describing the “bullet” itself (i.e., the contact discontinuity or “cold front”) or the subtle additional cold front between the bullet and the main shock reported by Markevitch (2006) and Markevitch & Vikhlinin (2007). The former lies outside the ALMA field of view, precluding any interesting constraints on the pressure difference (or lack thereof) across the cold front, while the latter is intrinsically faint, and is expected to be in thermal pressure equilibrium. For simplicity, we thus assume these features have a negligible effect on the measurements of the shock itself, and therefore ignore them. Further, we assume a single power-law profile for the downstream electron pressure (see below). A future analysis, joint with X-rays, will allow more model freedom for trying to build description of such features.

2.2.2. Bulk pressure distribution

To model the pressure distribution in the downstream region, we employ a power law radial profile with slope α, centered along the merger axis at a distance from the shock nose equal to the front curvature radius r_c,

$$\begin{array}{c}\hfill P_{\mathrm{e}}(r)=x_p{P}_{\mathrm{e},\mathrm{us}}{(r/r_{\mathrm{c}})}^{\alpha }.\end{array}$$(1)

On the other hand, we consider the pre-shock pressure distribution to be relaxed, thus to be described by a spherical generalized Navarro-Frenk-White (gNFW) profile (Nagai et al. 2007):

$$\begin{array}{c}\hfill P_{\mathrm{e}}(r)\propto M_{500}^{a_p}{P}_{500}{P}_0{(c_{500}r/r_{500})}^{-c}{[1+{(c_{500}r/r_{500})}^a]}^{(c-b)/a},\end{array}$$(2)

where r₅₀₀ and P₅₀₀ are functions of M₅₀₀, the total mass contained within an average overdensity 500× the critical density of the Universe at that redshift (Arnaud et al. 2010). We constrain M₅₀₀ so that the gNFW model always reproduces the X-ray value for the upstream pressure, P_e, us. For the main results reported here, the pressure normalization P₀, concentration parameter c₅₀₀, mass-dependence index a_P, and slopes (a, b, c) are fixed to the values reported in Arnaud et al. (2010) for the universal pressure profile. However, we show below that our results for ℳ are insensitive to the choice of gNFW parametrization and position of the gNFW model centroid. We therefore simply fix the gNFW centroid’s coordinates (RA, Dec)_gNFW to the position of main lensing κ-map peak inferred by Clowe et al. (2006), which we note does not coincide with the center of the post-shock profile.

Though pressure perturbations driven by the primary merger are confined to be within the shock front, it is possible for the passage of its associated dark matter component to affect, through infall, the ICM ahead of the shock (Springel & Farrar 2007). This may undermine our choice of the universal gNFW profile, reliable in the case of relaxed clusters, when describing the bulk pressure distribution. We tested against possible systematics introduced by this assumption. We found no significant deviations in the reconstructed parameters after changing either the slopes of the profiles or the position of the assumed centroids. The same applies to the structure of the gas on the downstream side as long as it is smooth, even though it may differ from expectations for a solid body moving through homogeneous fluid (see, e.g., Zhang et al. 2019). The results of the above tests are summarized in Fig. 2. The net result is that, as a consequence of the interferometric filtering, our shock model is largely sensitive to the pressure conditions right at the front, and not to the properties of the bulk pressure distribution (see Sect. 3.2). Thus, we consider wide uninformative priors on both the Mach number and post-shock slope, and marginalize over the latter.

Fig. 2. Bivariate posterior density function (left panel) for the inferred Mach number ℳ and slope of the downstream pressure distribution α for a set of upstream pressure profiles and respective centroids. To facilitate comparison, the best-fitting ℳ estimates are also plotted (right panel). We consider three different cases for the underlying gNFW profile by setting the slopes equal to the values reported in Arnaud et al. (2010) for the cool-core (cc; red), morphologically-disturbed (md; peach), or ensemble (up; blue) cluster samples. For each of them, the distribution centroid is then fixed to a number of different positions: far downstream and far upstream (right panel, left and mid-left points), respectively to arbitrary distances of 10′ east and 3′ west of the X-ray-derived shock nose coordinates; APEX-SZ centroid (mid-right point; Halverson et al. 2009); peak of the κ-map (right point; Clowe et al. 2006). In all the above cases, we assumed instantaneous shock heating of the electrons. The gray line in the right panel denotes the corresponding best-fitting ℳ reported in Sect. 3.1, while the darker and lighter bands the respective 68% and 95% credibility intervals.

In order to account for the high temperatures measured in the system, relativistic corrections to the SZ spectrum (Itoh & Nozawa 2004) are included in our modelling. In fact, variations in the measured SZ signal of the order of 5% up to 15% are expected for an electron gas with temperature ranging from 9 keV to 30 keV as measured from the X-ray data. As for P_e, us, we employ an X-ray-motivated prior on the upstream temperature T_e, us. In addition, we incorporate the 5% uncertainties on the ACA and ALMA flux calibration by introducing normalization hyperparameters κ_ACA and κ_ALMA (Di Mascolo et al. 2019). For all the modelling runs presented in the following sections, κ_ACA and κ_ALMA have been found to not deviate significantly from unity.

Given the plane-of-sky geometry of the merger involving the Bullet Cluster, any contribution from the kinetic SZ effect (Sunyaev & Zeldovich 1980) to the observed signal due to the motion of the single subclusters should be subdominant with respect to the thermal SZ effect. In fact, if we assume the velocity v ≈ 3000 inferred from the shock Mach number (Springel & Farrar 2007) to be measured with respect to the CMB rest frame and the merger direction to be oriented by around 8° with respect to the plane of sky (Markevitch et al. 2004), we find that the contribution of the kinetic SZ effect to the total SZ signal from the post-shock region would be of the order of 3% of the corresponding thermal component. This would induce a systematic error on the estimate of the Mach number ℳ lower than 2%. Given the small effect as well as the lack of robust constraints on the merger proper velocities and orientation, we then decide to not include the kinetic SZ effect in our model, keeping the merger axis aligned with plane of sky (see Sect. 2.2.1).

2.3. Implementation details

As noted earlier, the modelling algorithm and its specific implementation are detailed in Di Mascolo et al. (2019, see also references therein). However, we improved the posterior sampling algorithm by adopting the dynamic nested sampling by Higson et al. (2018). In particular, we employ the pure-Python implementation provided by dynesty (Speagle 2019).

3. Results

3.1. Instantaneous electron-ion temperature equilibration

We first consider a model for the shock front under the standard assumption of an instantaneous electron-ion temperature equilibration, i.e., T_e = T_i, at the shock front. This is consistent with the Chandra analysis by Markevitch (2006), who derived ℳ_X = 3.0 ± 0.4 from the density jump, and measured an electron temperature jump as expected if electrons reach the average post-shock temperature near-instantly. A more recent analysis of the same Chandra data by Markevitch (in prep.) provides an electron density jump x_n = 2.86 ± 0.16, assuming instantaneous electron heating (T_e, ds ≈ 30  keV) for converting the X-ray surface brightness to the density jump by accounting for the different emissivity in the post- and pre-shock regions (see, e.g., Ettori 2000 for a discussion about X-ray brightness modelling in the presence of temperature gradients). This corresponds to ℳ_X = 2.74 ± 0.25. The slight difference from the older estimate is due to a better-centered shock model and the inclusion of the azimuthal decline in x_n for angles away from the shock nose when projecting along the l.o.s. (as in Wang et al. 2018), based on the amplitude of the density jumps measured from the X-ray data in different sectors of the shock (Markevitch, in prep.).

The assumption of instantaneous heating implies that the nominal Rankine-Hugoniot condition can be used to relate ℳ to the measured amplitude of the electron pressure jump x_P relative to the upstream value P_e, us as

$$\begin{array}{c}\hfill \mathcal{M}={\left[\frac{(\gamma +1)x_p+(\gamma -1)}{2\gamma }\right]}^{1/2}.\end{array}$$(3)

Here, γ is the polytropic exponent, which we assume to be γ = 5/3, appropriate for non-relativistic fully-ionized gas. We further allow for the azimuthal variation of ℳ along the shock front, whose scaling with the azimuthal angle is derived using the same density jump decline discussed above. The omission of such azimuthal dependence would cause the Mach number to be averaged down with respect to its maximum value due to the effect of the wings with lower ℳ. For the results provided in this and the following sections, we estimate that the inclusion of the X-ray-based model for the azimuthal variation of the shock pressure jump increases the value of the inferred Mach number by only 5–7%. However, more severe effects should be expected for observations with larger field of views, which would include values from farther in the wings.

We obtain $\mathcal{M}=2.{08}_{-0.12}^{+0.12}$ (Fig. 3). While the model relies on the X-ray priors on the pre-shock pressure and temperature, the derived Mach number is inconsistent at a 2.4σ level with the X-ray estimate ℳ_X = 2.74 ± 0.25. Projection effects may play a non-negligible role in biasing the SZ-based measurement of the Mach number. However, a strong ellipticity of the shock front shape ${\ell }_{\mathrm{LoS}}^{\mathrm{true}}/{\ell }_{\mathrm{LoS}}\lesssim 0.6$ (see Sect. 2.2.1) would be required to bridge the gap between SZ and X-ray estimates. In reality, an even larger ellipticity would be necessary, given that the X-ray estimates would also be affected by geometry, albeit with a different dependence. Another potential source of bias is the X-ray-motivated prior on the upstream pressure, which comes from deprojected density and temperature estimates, used to compute the relative pressure jump. While the definition of a centroid for X-ray deprojection remains ambiguous, we found only extreme choices would alter our results significantly. A joint-likelihood X-ray+SZ analysis may be required to find a consistent geometry that fully reconciles such discrepancies.

3.2. Collisional electron-ion temperature equilibration

Here we consider the possibility that the electron and ion temperatures do not equilibrate instantaneously in cluster shocks (i.e., T_e ≠ T_i immediately inside the shock front; see, e.g., Fox & Loeb 1997; Markevitch 2006; Russell et al. 2012; Wang et al. 2018). Ions carry the majority of the gas bulk kinetic energy in collisionless shocks, and are heated dissipatively on scales comparable to their gyro-radii, while electrons might remain much colder (Vink et al. 2015), unless there is some process that equilibrates the ion and electron temperatures. The upper limit on the equilibration time scales is set by Coulomb collisions (Zeldovich & Raizer 1966), which for the downstream density and temperature in the Bullet Cluster is long (∼few 10⁸ yr), occurring over a distance comparable to the offset between the shock and the cold front.

Under the assumption of conservation of the enthalpy flux, electrons equilibrate with ions to the Rankine-Hugoniot downstream temperature at a rate driven by Coulomb collisions (Fox & Loeb 1997)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{T}_{\mathrm{e}}}{\mathrm{d}t}=\frac{1}{t_{\mathrm{eq}}}(1+\frac{n_{\mathrm{e}}}{n_{\mathrm{i}}})(x_t{T}_{\mathrm{e},\mathrm{us}}-T_{\mathrm{e}}),\end{array}$$(4)

where t_eq is the Coulomb collisional time-scale (Spitzer 1962), n_e and n_i are respectively the electron and ion densities, and x_T is the temperature ratio across the shock front. To build our SZ model, we convert the above equation in terms of the distance from the shock front by means of the downstream gas velocity u_ds = (ℳ/x_n) c_us, with c_us given by the upstream sound speed. Furthermore, we assume that electrons are first heated adiabatically (Vink et al. 2015), so that the electron temperature immediately inside the shock front equals $x_{\mathrm{n}}^{\gamma -1}{T}_{\mathrm{e},\mathrm{us}}$. As required by the conservation of charge neutrality across the shock front, the density jump x_n is also set to follow the Rankine-Hugoniot condition. In this case, the pressure jump x_P cannot be directly related to ℳ as in Eq. (3), and instead must be derived as the product of the density and temperature ratios at each three-dimensional model coordinate.

In the case of collisional equilibration, we find $\mathcal{M}=2.{53}_{-0.25}^{+0.33}$ (Fig. 3). This is consistent with ℳ_X = 2.57 ± 0.23, coming from the Chandra X-ray brightness fit if one uses the adiabatic-compression post-shock temperature (T_e, ds ≈ 20 keV) to convert to the density jump (Markevitch, in prep.). Unfortunately, due to the severe filtering of large spatial scales as well as the limited field of view, we are not able to put any significant constraint on the specific equipartition time-scale when treating t_eq as a free parameter. Instead, we find that assuming the electron-ion equilibration to be driven by Coulomb collisions is practically equivalent to setting t_eq = ∞.

Fig. 3. Marginalized posterior distributions for the shock Mach number ℳ derived under the assumptions of instantaneous (red) and collisional (blue) electron-ion equilibration. The dashed and dotted lines indicates the median of the posterior distributions and the 68% credible intervals, respectively.

For illustrative purposes, we present in Fig. 4 the dirty² images of the raw ALMA+ACA data employed in our analysis, the interferometric model corresponding to the collisional electron-ion equilibration scenario, and the respective model-subtracted data. As shown in the right panel, it is not possible to identify residuals that differ at a significant level from noise-like features.

Fig. 4. Dirty images of the raw (left), model (middle), and residual (right) ALMA+ACA interferometric data. They are generated by jointly gridding the ACA and ALMA data using a multi-frequency natural weighting scheme. For reference, this provides a synthesized beam of 4.07″ × 3.01″ FWHM (PA 81°; bottom-left corner of right panel). We cut off the fields at the 0.2 gain level of the joint ALMA+ACA antenna pattern. To better highlight the large-scale shock features, we avoid correcting for the primary beam attenuation and apply an additional 30 kλ taper. We note that the model subtraction is performed directly in Fourier space. The dashed line in the center and right panels indicate the inferred position of the shock front. X-ray contours are overlaid on the left panel from Fig. 1. We also note the positive signal at the shock front is not due to an increment of the SZ signal; rather it is an artifact of the high-pass filtering effects of ALMA+ACA (see also Fig. 5).

The analogous image for the instantaneous case is visually identical to Fig. 4, and therefore is not shown. This is evident in Fig. 5, where it is not possible to identify any significant difference between the filtered SZ models for the instantaneous- and collisional-equilibration scenarios. This confirms that the ALMA+ACA data are only sensitive to the properties of the pressure distributions near the shock edge, thus providing a direct estimate of the pressure difference across the shock front, i.e., ΔP_e rather than x_P. Moreover, our result is found to be practically independent of the specific assumption about the underlying gNFW profile, which is entirely filtered by the interferometric response (Fig. 5).

Fig. 5. Comparison of SZ signal profiles (upper panel) and the corresponding electron pressure profiles (lower panel) across the shock nose. The vertical gray line denotes the position of the shock front. We note that the upper panel contains both unfiltered input model fits for the SZ signal and the corresponding filtered (observed) profiles (see legend for details). As discussed in Sect. 3.2, the instantaneous (dash-dotted yellow) and collisional (solid red) shock models are indistinguishable after the spatial filtering of the interferometric ALMA+ACA observation, reflecting the fundamental limitation of ALMA+ACA to constrain any large-scale (≳1′) component of the SZ signal. Again, we note that while we measure a decrement due to the SZ effect, the filtered (observed) profiles can exhibit both positive and negative excursions (analogous to the Gibbs phenomenon; see, e.g., Bracewell 1978). For comparison, we also report the input unfiltered model (solid yellow) for the instantaneous equilibration, and both the raw and filtered underlying gNFW profile (dashed and dotted blue lines, as noted in the legend). We also present the X-ray expectation for both the filtered and raw SZ signal profiles (dashed and solid green lines) and corresponding pressure profiles that we would expect for the case of instantaneous equilibration, using the value ℳX = 2.74 ± 0.25 derived from fits to the X-ray data.

The fact ALMA+ACA is only sensitive in practice to the electron pressure difference ΔP_e across the shock is also reflected in the lack of any significant difference between the Bayesian evidences of the instant equilibration and adiabatic heating models ($\mathrm{\Delta }log\mathcal{Z}=1.{30}_{-1.48}^{+0.71}$). In fact, given that we cannot observe any large-scale feature in the SZ signal induced by the slow increase of the post-shock electron temperature in the case of collisional equilibration, the different heating scenarios practically differ only in the way we convert the pressure difference to an estimate of the shock ℳ.

We recall here a subtle cold front is observed in the Chandra X-ray image between the shock and the bullet boundary, ∼15″ east from the shock. Since the total pressure across a cold front is expected to be approximately continuous, this was not included in our modeling of the downstream pressure profile. In fact, no apparent signature of such feature can be distinguished in the ALMA+ACA SZ observation. However, because the flow of the post-shock gas would not cross the cold front, it is unlikely that any electron-ion temperature non-equilibrium would extend past it. If indeed there is an electron-proton temperature difference in the post-shock region, we should expect the electron temperature (and, hence, pressure) to reach its equilibrium value. This would in turn result in a discontinuity in the SZ signal in the direction of the front itself. Future, more sensitive ALMA observations may search for such a feature.

4. Conclusions

We further demonstrate the ability of using deep, high-resolution ALMA+ACA observations of the SZ effect to characterize shocks in merging clusters (see, for comparison, Basu et al. 2016). For this purpose, we studied the SZ effect across the shock in the Bullet Cluster, chosen as it is widely regarded as the “textbook example” of a cluster merger bow shock.

The application of our interferometric modelling technique – using X-ray-motivated priors – has allowed us to place constraints on the electron pressure discontinuity across the shock. Assuming a Rankine-Hugoniot shock adiabat, our pressure jump implies a Mach number $\mathcal{M}=2.{08}_{-0.12}^{+0.12}$, which is significantly lower than the one derived from Chandra data using the same geometric assumptions (ℳ = 2.74 ± 0.25). An interesting physical possibility to reconcile the two measurements is to allow that the electron and ion temperatures do not equilibrate instantly after the shock passage has heated the electrons adiabatically. For a given Mach number, this would lower the post-shock electron temperature and thus the observed electron pressure jump. Our Mach number would then become $\mathcal{M}=2.{53}_{-0.25}^{+0.33}$, in agreement with the X-ray estimate that assumes the adiabatic temperature jump for conversion between the X-ray brightness and density. We note that Chandra X-ray data constrain the gas density (from which the Mach number is derived) and electron temperature across the shock separately, and its post-shock temperature prefers instant equilibration over adiabatic heating of the electrons (at ∼2σ confidence; Markevitch 2006). However, while Chandra is free from ALMA+ACA interferometric limitations and can probe the upstream and downstream gas directly, the Bullet post-shock temperature is above the range where Chandra can measure electron temperatures reliably, and hence suffer significant systematic uncertainties.

To summarize, ALMA+ACA has proven to provide a clean measurement of the differential jump in pressure due to the shock, and, in combination with data that can access larger scales, can provide compelling constraints on shock properties such as the Mach number. In particular:

• Interferometric observations cleanly measure the projected pressure jump due to the shock. However, due to the inherent spatial filtering of ALMA and the ACA, which recover scales ∼0.5−1.1′ in Band 3 (compared to θ₅₀₀∼10′, corresponding to r₅₀₀ for the Bullet Cluster), X-ray priors on both the model geometry and upstream pressure are necessary in order to infer ℳ. On the other hand, the combination of SZ observations covering a broader range of spatial scales (i.e., from 0.1″ − 30′) could provide an SZ-only view of the shock properties.

• Once the geometry is fixed, the key quantities which drive the analysis and the interpretations of the results are the pressure difference across the shock front, the normalization of the pre-shock pressure, and the independent X-ray estimates of ℳ. We show their statistical uncertainties are small enough to allow us to differentiate between the instantaneous and adiabatic heating scenarios. Nevertheless, neither model is unambiguously preferred. Although the two scenarios result in SZ-based estimates for ℳ that deviate one from the other by ∼2σ, we find that the difference of the respective Bayesian log-evidence is not significant enough to completely rule out one versus another.

• We extensively tested our modelling choices – varying the geometry, pre- and post-shock pressure slopes, and underlying pressure distribution – and find our results to be robust for a broad range of possible assumptions motivated by the X-ray analyses. However, our model does not fully describe the complex morphology observed in the X-ray surface brightness. Together with the uncertainties on the three-dimensional morphology of the cluster, this may limit our ability to elucidate the nature of electron heating across the shock front.

Together, these illustrate the fundamental complementarity of X-ray and SZ effect observations in the study of the physics of galaxy clusters. It is then clear that a simultaneous, joint-likelihood analysis of the SZ and X-ray data on the Bullet Cluster, extending the approach of incorporating X-ray information in the form of priors (see our discussion in Sect. 2.2.1, and Sect. 3.3 of Basu et al. 2016), would benefit our understanding of the morphology of the galaxy cluster, as well as provide further insights into the physical mechanisms for shock heating of the intracluster medium. A forthcoming paper will present the results of a full joint-likelihood analysis of interferometric SZ and X-ray observations, as well as single-dish SZ measurements, building on the methodology discussed in Di Mascolo et al. (2019). Meanwhile, upcoming results from NuSTAR (Wik et al. in prep.) will better access the high photon energies corresponding to the high temperatures inferred from Chandra. Further, the number of observations of shocks with unambiguous geometry and sufficiently high Mach number that allow the detection of deviations from instantaneous electron heating is limited. Thus, along with improved modelling, observations of a larger sample of cluster shocks will be needed to improve our understanding. And finally, both improved spatial and spectral resolution, larger instantaneous field of view, and the ability to recover zero-spacing information will vastly improve future SZ-only studies. However, in order to provide sufficient overlap with the interferometric data in Fourier space, while also probing higher frequencies and spatial scales > 10′, a new wide-field (> 1°) single-dish facility, such as the Atacama Large Aperture Submillimeter/millimeter Telescope (AtLAST; see, e.g., Klaassen et al. 2019; Mroczkowski et al. 2019b) is required.

Acknowledgments

We thank the anonymous referee for the insightful comments. EC and RS acknowledges partial support from the Russian Science Foundation grant 19-12-00369. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00760.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. Basic research in radio astronomy at the US Naval Research Laboratory is funded by 6.1 Base funding, supporting T. Clarke. The scientific results reported in this article are obtained using the following softwares: AstroPy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2018); dynesty (Speagle 2019); galario (Tazzari et al. 2018).

References

Appendix A: Interferometric view of the SZ effect

The amplitude of the thermal SZ effect in the direction x of a galaxy cluster characterized by an electron pressure distribution P_e is proportional to Compton y,

$$\begin{array}{c}\hfill (\mathit{x})=\frac{{\sigma }_t}{m_{\mathrm{e}}{c}^2}{\displaystyle \int }{P}_{\mathrm{e}}(\mathit{x},\ell )\mathrm{d}\ell ,\end{array}$$(A.1)

where σ_T, m_e, and c are respectively the Thomson cross-section, the electron mass, and the speed of light, while ℓ is the coordinate along the line of sight.

As shown in Sect. 3.2 and Fig. 5, the available ALMA+ACA data (i) only probe a small region near the tip of the shock and (ii) any extended structures in this region are effectively filtered out. Therefore, a model that includes a jump of pressure across the shock front and is smooth otherwise should capture the most of the information contained in the ALMA+ACA data. To this end, we represent the electron pressure distribution P_e(x, ℓ) as a combination of two spatially smooth components, P_e, 0(x, ℓ) and P_e, 1(x, ℓ), where their sum is

$$\begin{array}{c}\hfill P_{\mathrm{e}}=(1-f_{\mathrm{ds}})P_{\mathrm{e},0}+f_{\mathrm{ds}}{P}_{\mathrm{e},1}.\end{array}$$(A.2)

Here f_ds ≡ f_ds(x, ℓ) is equal to 1 in the downstream region and 0 in the upstream region (i.e., is a Heaviside step function). We note both that here and below, for simplicity, we omit the explicit coordinate dependence of P_e or y. Rearranging the terms and integrating along the line of sight we obtain

$$\begin{array}{c}\hfill y\propto \int P_{\mathrm{e},0}\mathrm{d}\ell +\int (P_{\mathrm{e},1}-P_{\mathrm{e},0})f_{\mathrm{ds}}\mathrm{d}\ell .\end{array}$$(A.3)

Since the first term in the above expression corresponds to a smooth, large-scale pressure distribution its contribution to y is filtered out from the ALMA+ACA data. Moreover, the function (P_e, 1 − P_e, 0) in the second term is also spatially smooth and would be filtered too without the step function f_ds. Therefore, the signal $\stackrel{\sim }{y}$ measured by ALMA+ACA in the vicinity of the shock tip is effectively defined by the second term in Eq. (A.3), which is set by the pressure jump at the shock front and the length-scale of the downstream region. Thus,

$$\begin{array}{c}\hfill \stackrel{\sim }{\propto }\int (P_{\mathrm{e},1}-P_{\mathrm{e},0})f_{\mathrm{ds}}\mathrm{d}\ell \approx (P_{\mathrm{e},\mathrm{ds}}-P_{\mathrm{e},\mathrm{us}}){\ell }_{\mathrm{LoS}}=\mathrm{\Delta }{P}_{\mathrm{e}}{\ell }_{\mathrm{LoS}},\end{array}$$(A.4)

where ℓ_LoS is the line-of-sight extent of the probed post-shock region, and P_e, us and P_e, ds are the electron pressures measured just outside and inside the shock front, respectively (Sect. 2.2). Therefore, ALMA+ACA data effectively constrain a product of the electron pressure difference at the shock ΔP_e and the physical size of the region ℓ_LoS. The latter quantity can be easily determined if the merger is in the plane of the sky and the shock front possesses rotational symmetry. If the shape of the front can be approximated by a sphere with a radius R, then along the symmetry axis ${\ell }_{\mathrm{LoS}}\approx \sqrt{2rh}$, where h ≪ r is the distance from the tip of the shock. While the calculations in the paper were done without these simplifying assumptions, the Eq. (A.4) is useful to estimate the uncertainty introduced by the (unknown) geometry of the shock along the line of sight. In particular, if the curvatures in the sky plane r_PoS and along the line of sight r_LoS differ, the estimate of the pressure difference ΔP_e, which assumes r_LoS = r_Pos, will be biased by a factor $\sqrt{r_{\mathrm{LoS}}/r_{\mathrm{PoS}}}$ (the same argument is discussed in Wang et al. 2018).

All Tables

All Figures

Fig. 1. Cut-out of the 0.5−2.0 keV Chandra X-ray surface brightness map of the Bullet Cluster. The solid wedges represent the sectors employed to derive the hyperbolic shape best-matching the shock front geometry. The maximum-a-posteriori model is shown as a dashed line, while the dotted contours indicate the corresponding 95% credible interval. Further, the dashed and solid yellow circles denote respectively the ACA and ALMA full-width-at-half-maximum fields of view. For reference, the upstream and downstream gas lie respectively west (right) and east (left) of the shock front.

In the text

Fig. 2. Bivariate posterior density function (left panel) for the inferred Mach number ℳ and slope of the downstream pressure distribution α for a set of upstream pressure profiles and respective centroids. To facilitate comparison, the best-fitting ℳ estimates are also plotted (right panel). We consider three different cases for the underlying gNFW profile by setting the slopes equal to the values reported in Arnaud et al. (2010) for the cool-core (cc; red), morphologically-disturbed (md; peach), or ensemble (up; blue) cluster samples. For each of them, the distribution centroid is then fixed to a number of different positions: far downstream and far upstream (right panel, left and mid-left points), respectively to arbitrary distances of 10′ east and 3′ west of the X-ray-derived shock nose coordinates; APEX-SZ centroid (mid-right point; Halverson et al. 2009); peak of the κ-map (right point; Clowe et al. 2006). In all the above cases, we assumed instantaneous shock heating of the electrons. The gray line in the right panel denotes the corresponding best-fitting ℳ reported in Sect. 3.1, while the darker and lighter bands the respective 68% and 95% credibility intervals.

In the text

	Fig. 3. Marginalized posterior distributions for the shock Mach number ℳ derived under the assumptions of instantaneous (red) and collisional (blue) electron-ion equilibration. The dashed and dotted lines indicates the median of the posterior distributions and the 68% credible intervals, respectively.
In the text

Fig. 4. Dirty images of the raw (left), model (middle), and residual (right) ALMA+ACA interferometric data. They are generated by jointly gridding the ACA and ALMA data using a multi-frequency natural weighting scheme. For reference, this provides a synthesized beam of 4.07″ × 3.01″ FWHM (PA 81°; bottom-left corner of right panel). We cut off the fields at the 0.2 gain level of the joint ALMA+ACA antenna pattern. To better highlight the large-scale shock features, we avoid correcting for the primary beam attenuation and apply an additional 30 kλ taper. We note that the model subtraction is performed directly in Fourier space. The dashed line in the center and right panels indicate the inferred position of the shock front. X-ray contours are overlaid on the left panel from Fig. 1. We also note the positive signal at the shock front is not due to an increment of the SZ signal; rather it is an artifact of the high-pass filtering effects of ALMA+ACA (see also Fig. 5).

In the text

Fig. 5. Comparison of SZ signal profiles (upper panel) and the corresponding electron pressure profiles (lower panel) across the shock nose. The vertical gray line denotes the position of the shock front. We note that the upper panel contains both unfiltered input model fits for the SZ signal and the corresponding filtered (observed) profiles (see legend for details). As discussed in Sect. 3.2, the instantaneous (dash-dotted yellow) and collisional (solid red) shock models are indistinguishable after the spatial filtering of the interferometric ALMA+ACA observation, reflecting the fundamental limitation of ALMA+ACA to constrain any large-scale (≳1′) component of the SZ signal. Again, we note that while we measure a decrement due to the SZ effect, the filtered (observed) profiles can exhibit both positive and negative excursions (analogous to the Gibbs phenomenon; see, e.g., Bracewell 1978). For comparison, we also report the input unfiltered model (solid yellow) for the instantaneous equilibration, and both the raw and filtered underlying gNFW profile (dashed and dotted blue lines, as noted in the legend). We also present the X-ray expectation for both the filtered and raw SZ signal profiles (dashed and solid green lines) and corresponding pressure profiles that we would expect for the case of instantaneous equilibration, using the value ℳX = 2.74 ± 0.25 derived from fits to the X-ray data.

In the text