© ESO, 2017

1. Introduction

The assembly of massive clusters of galaxies represents the last step of the hierarchical formation of structures in the Universe. Galaxy clusters arise from smooth accretion of surrounding material, but also from the merging of subclusters and groups, at the intersection of filaments in the cosmic web (see Kravtsov & Borgani 2012, for a review). Merging clusters are not only useful as large scale astrophysical laboratories (e.g., Clowe et al. 2006), but are important from a cosmological point of view (e.g., Voit 2005). The physical processes involved in mergers significantly affect the cluster observables via the distribution of dark matter, the physics of the intracluster plasma, and the dynamics of galaxies (e.g., Ferrari et al. 2005; Maurogordato et al. 2011), i.e., all the components of the clusters. These observables are in turn essential to allow us to infer the masses of clusters when used as a cosmological probe (see Allen et al. 2011, for a review). In this context, measuring the collision velocity of merging systems provides insights of the cluster physics, such as shocks in the intracluster medium (ICM), the kinematics of the mergers, the hydrodynamical interactions between cluster cores and less dense shock heated regions, and the relation between the baryonic matter and the dark matter. In addition, such velocity measurements provide a complementary cosmological probe to constrain structure formation (Bhattacharya & Kosowsky 2007) and the study of the properties of the velocity field is a subject of active research (e.g., Hahn et al. 2015, using numerical simulations).

Measuring the peculiar velocities of galaxy clusters is a difficult task, in particular when looking at objects beyond the local Universe. Indeed, classical methods that are based on redshift surveys allow us to reconstruct the internal galaxy kinematics of clusters, but require an independent distance measurement to subtract the Hubble flow at objects’ positions in order to measure their true peculiar velocities (e.g., Tully & Fisher 1977). They are therefore limited to small redshifts as the errors scale with distance.

The Sunyaev-Zel’dovich effect (SZ, Sunyaev & Zeldovich 1972), on the other hand, provides a way to probe the ICM gas bound in clusters without suffering from cosmological dimming (Birkinshaw 1999; Carlstrom et al. 2002; Kitayama 2014). Its observation is therefore independent of redshift as long as the angular resolution of the observations is sufficient to resolve the clusters. We can distinguish between the thermal SZ effect (tSZ, Sunyaev & Zeldovich 1972), for which the cosmic microwave background (CMB) photons are spectrally distorted by the electronic thermal pressure of the ICM, and the kinetic SZ effect (kSZ, Sunyaev & Zeldovich 1980), arising from the CMB Doppler shift induced by the bulk motion of the cluster electrons. While the tSZ effect is related to the integrated pressure along the line-of-sight and can be used as a mass proxy in cosmological studies (see, e.g., Planck Collaboration XXIV 2016; Bleem et al. 2015; Hasselfield et al. 2013), the kSZ effect is sensitive to the integrated line-of-sight electron density and the gas velocity with respect to the CMB reference frame. The kSZ signal is subdominant to the tSZ unless the gas velocity reaches a few tenths of a percent of the speed of light (1000 km s^-1, see e.g., Birkinshaw 1999). In addition, the spectral dependence of the kSZ is the same as that of the CMB, so that they can only be separated from each other using the difference in spatial distribution of the respective signals. Direct observation of the kSZ signal within clusters therefore requires both high sensitivity and high angular resolution observations (Haehnelt & Tegmark 1996), which makes it particularly challenging to detect in individual clusters (see attempts on individual clusters by e.g., Holzapfel et al. 1997; Benson et al. 2003; Sayers et al. 2016; and the first clear detection on a statistical sample by Hand et al. 2012).

The cluster of galaxies MACS J0717.5+3745 at z = 0.55 is a striking example of a merging system and is an excellent target with which to measure the kSZ effect. A detailed review of the previous analyses of MACS J0717.5+3745 can be found in Mroczkowski et al. (2012) and Sayers et al. (2013), and we summarize some essential features here. The cluster is located at the end of a filamentary structure (e.g., Ebeling et al. 2004; Jauzac et al. 2012) and appears to be an extremely disturbed system. It is composed of four optically-identified main subclusters (Ma et al. 2009), which coincide with peaks in the observed strong lensing surface mass distribution (see e.g., Zitrin et al. 2009; Diego et al. 2015; Limousin et al. 2016; Kawamata et al. 2016). The X-ray emitting gas is very disturbed (e.g., Ebeling et al. 2007; Ma et al. 2009) and the spectroscopic temperature reaches up to ≳20 keV (Ma et al. 2009), indicating a scenario where two objects are undergoing a major merging event. Edge et al. (2003) discovered the presence of powerful radio emission near the high temperature region, indicating the presence of non-thermal component (see also, e.g., van Weeren et al. 2009). MACS J0717.5+3745 was also observed using 30 GHz interferometric SZ data by the Sunyaev-Zel’dovich Array (SZA), indicating that the cluster is a hot and massive system (LaRoque et al. 2003). Mroczkowski et al. (2012) reported high angular resolution SZ observations at 90 GHz using MUSTANG (The MUltiplexed SQUID/TES Array at Ninety GHz) at the Green Bank Telescope (GBT). The data revealed small scale structures in the pressure distribution, as expected from a merger scenario.

The velocity distribution of MACS J0717.5+3745 was first studied using optical spectroscopy (Ma et al. 2009), showing that the galaxy groups have exceptionally large relative line-of-sight velocities, providing a more detailed picture of the merger kinematics. Using Bolocam SZ data obtained at the Caltech Submillimeter Observatory (CSO) at 140 and 268 GHz, Mroczkowski et al. (2012) reported the first indication of the presence of kSZ signal in MACS J0717.5+3745. Sayers et al. (2013) then reported the first significant detection of the kSZ signal in a single cluster, MACS J0717.5+3745, after collecting additional data with Bolocam and refining the former analysis. As the Bolocam beams are larger than the target subclusters, the limited angular resolution of these data is a factor impacting the kSZ detection significance (see Sayers et al. 2013, for more details).

In this paper, we report deep, high angular resolution observations (<20 arcsec) of MACS J0717.5+3745 at 150 and 260 GHz, obtained at the Institut de Radio Astronomie Millimétrique (IRAM) 30 m telescope using the New IRAM KID Arrays (NIKA, Monfardini et al. 2011; Calvo et al. 2013; Adam et al. 2014; Catalano et al. 2014). The resolved nature of the observations allows us for a robust treatment of the contaminating astrophysical signal and the data are used to produce a map of the kSZ signal toward the cluster. Combining X-ray spectroscopic temperature measurements with NIKA, we also fit a parametric model to constrain the gas line-of-sight velocity distribution.

The paper is organized as follows. In Sect. 2, we present the SZ observations of MACS J0717.5+3745. The removal of the astrophysical contamination is presented and discussed in Sect. 3. In Sect. 4, we detail the X-ray data reduction. In Sect. 5, we present the mapping of the kSZ effect in MACS J0717.5+3745 and we use our data together with X-ray temperatures to constrain the gas line-of-sight velocity in Sect. 6. Throughout this paper we assume a flat ΛCDM cosmology according to the latest Planck results (Planck Collaboration XIII 2016) with H₀ = 67.8 km s^-1 Mpc^-1, Ω_M = 0.308, and Ω_Λ = 0.692.

2. Sunyaev-Zel’dovich observations and data reduction

2.1. The Sunyaev-Zel’dovich effect

The SZ effect is the spectral distortion of the CMB due to the inverse Compton scattering of CMB photons off energetic electrons in galaxy clusters. Only the thermal and the kinetic SZ effects have been observed so far (induced by a thermal electron population and the bulk motion of electrons with respect to the CMB, respectively), but in principle the SZ effect can also arise from the CMB interaction with more exotic electron populations, such as those originating from dark matter annihilation (Colafrancesco 2004) or ultra-relativistic electrons accelerated in radio lobes (Colafrancesco 2008). These other contributions are expected to be small and, in this paper, we assume that only the tSZ and the kSZ effects contribute to the SZ signal.

The tSZ effect follows the characteristic frequency dependence given by (Birkinshaw 1999) $\mathit{f}\mathrm{\left(}\mathit{x,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{x}}^{\mathrm{4}}{\mathrm{e}}^{\mathit{x}}}{{\mathrm{(}{\mathrm{e}}^{\mathit{x}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{\mathrm{2}}}\left(\mathit{x}\coth \left(\frac{\mathit{x}}{\mathrm{2}}\right)\mathrm{-}\mathrm{4}\right)\left(\mathrm{1}\mathrm{+}{\mathit{\delta }}_{\mathrm{tSZ}}\mathrm{\left(}\mathit{x,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\right)\mathit{,}$(1)where $\mathit{x}\mathrm{=}\frac{\mathit{h\nu }}{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{CMB}}}$ is the dimensionless frequency, with ν the observing frequency, h the Planck’s constant, k_B the Boltzmann’s constant, and the CMB temperature T_CMB = 2.725 K (Fixsen 2009). The term δ_tSZ(x,T_e) corresponds to relativistic corrections, which depend on the observing frequency and the electron temperature T_e, and which we compute using the work of Itoh & Nozawa (2003). These corrections are valid at a level of ≲5% up to 50 keV within the NIKA bands.

The spectral dependance of the kSZ effect follows that of the CMB in the non relativistic regime. It is given by (Birkinshaw 1999) $\mathit{g}\mathrm{\left(}\mathit{x,}{\mathit{v}}_{\mathit{z}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{x}}^{\mathrm{4}}{\mathrm{e}}^{\mathit{x}}}{{\mathrm{(}{\mathrm{e}}^{\mathit{x}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}}^{\mathrm{2}}}\left(\mathrm{1}\mathrm{+}{\mathit{\delta }}_{\mathrm{kSZ}}\mathrm{\left(}\mathit{x,}{\mathit{v}}_{\mathit{z}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\right)\mathit{,}$(2)where v_z is the gas line-of-sight velocity (positive for a cluster receding from the observer) and δ_kSZ are the relativistic corrections that depend on the observing frequency, T_e, and v_z. We use the analytical formula provided by Nozawa et al. (2006) to compute relativistic corrections, which we expect to be more accurate than ~10% up to 30 keV in the NIKA bands. The velocity dependence on these corrections is negligible because v_z is always much smaller than the speed of light.

The observed change of specific intensity expressed as a function of ν, ΔI_ν, with respect to the CMB, I₀, is then given by $\frac{\mathrm{\Delta }{\mathit{I}}_{\mathit{\nu }}}{{\mathit{I}}_{\mathrm{0}}}\mathrm{=}\mathit{f}\mathrm{\left(}\mathit{\nu ,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}{\mathit{y}}_{\mathrm{tSZ}}\mathrm{+}\mathit{g}\mathrm{\left(}\mathit{\nu ,}{\mathit{T}}_{\mathrm{e}}\mathit{,}{\mathit{v}}_{\mathit{z}}\mathrm{\right)}{\mathit{y}}_{\mathrm{kSZ}}\mathit{,}$(3)where the amplitude of the tSZ and the kSZ signal are given by ${\mathit{y}}_{\mathrm{tSZ}}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{T}}}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\mathrm{\int }{\mathit{P}}_{\mathrm{e}}\mathrm{d}\mathit{l}\mathrm{\equiv }\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{e}}}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\mathit{\tau }$(4)and ${\mathit{y}}_{\mathrm{kSZ}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{T}}\mathrm{\int }\frac{\mathrm{-}{\mathit{v}}_{\mathit{z}}}{\mathit{c}}{\mathit{n}}_{\mathrm{e}}\mathrm{d}\mathit{l}\mathrm{\equiv }\frac{\mathrm{-}{\mathit{v}}_{\mathit{z}}}{\mathit{c}}\mathit{\tau ,}$(5)respectively. The parameter σ_T is the Thomson cross section, m_e is the electron rest mass, c is the speed of light, and n_e is the electronic gas density. The kSZ and tSZ amplitudes can also be expressed using the optical depth, τ, as given in Eqs. (4) (assuming constant temperature along the line-of-sight) and (5) (assuming constant velocity along the line-of-sight). As the ICM is a dilute medium, we assume the ideal gas law to hold and express the electronic pressure as P_e = n_ek_BT_e. The functions f(ν,T_e) and g(ν,T_e,v_z) are integrated over the NIKA bands, and they are listed in Table 1 for a few electronic gas temperature values, accounting for the NIKA beam solid angle to express them in terms of Jy/beam. While the tSZ signal is always negative at 150 GHz and positive at 260 GHz, the kSZ signal can be either positive or negative, depending on the sign of the velocity along the line-of-sight, but its sign is the same at both NIKA frequencies. The kSZ maximum corresponds to the null of the tSZ signal and is located around 217 GHz if relativistic effects are negligible.

2.2. NIKA observations and data reduction

Fig. 1 Top: NIKA maps of MACS J0717.5+3745 at 150 GHz (left) and 260 GHz (right). The gray contours show the significance in units of standard deviation. They are multiples of 2σ at 150 GHz and 1σ at 260 GHz, starting at ±2σ. Both maps have been smoothed to have the same effective resolution of 22 arcsec FWHM, as represented by the white circle on the bottom left corner of the maps. The regions defined in Table 2 are shown as gray dot-dashed circles. Bottom: same as top-panel figures, in the case of half difference maps computed from two equivalent sub-samples. The steps between contours is 1σ for both frequencies.

MACS J0717.5+3745 was observed for 7.26 h in February 2014 and 5.85 h in January and February 2015. Most of the observations took place under good conditions, with a mean zenith opacity of 0.116 at 260 GHz and 0.085 at 150 GHz, and a stable atmosphere. The pointing center was (RA, Dec)_J2000 = (07:17:32.3, + 37:44:47) for the 2014 data set and (RA, Dec)_J2000 = (07:17:32.3, + 37:45:10) for 2015. The scanning strategy was the same as that employed for other NIKA clusters as detailed in Adam et al. (2015, 2016). Similarly, the detailed calibration procedure can be found in Adam et al. (2014, 2015). This results in a root mean square pointing error of ≲3 arcsec, the absolute calibration uncertainties are 7% at 150 GHz and 12% at 260 GHz, and the Gaussian beams FWHM were measured to be 18.2 and 12.0 arcsec at 150 and 260 GHz, respectively.

The data were reduced similarly and independently for the two bands, as described in Adam et al. (2015). The removal of the correlated noise across the arrays, due mainly to the atmospheric contribution and correlated electronic noise, leads to the filtering of the astrophysical signal at scales that are larger than the NIKA field of view (~2 arcmin). The resulting transfer function was computed using simulations. It is close to unity at small scales and vanishes smoothly at scales larger than the field of view. Apart from beam smoothing, the transfer function is the same at 150 and 260 GHz, which allows for direct combination of the two maps. The absolute zero level brightness of the maps is not measured by NIKA. In practice, the effective zero level is computed during the mapmaking procedure by setting the mean of the time ordered data to zero, taken in regions that are outside the cluster, flagged in an iterative manner. The flagged cluster region is based on the 150 GHz map only because of the higher signal-to-noise ratio, and corresponds to the 5σ contours of Fig. 1, top left panel. The beam efficiency correction as a function of elevation was optimized for extended emission (Greve et al. 1998), which can lead to up to ~10% loss in the flux of point sources for very high or very low elevations. However, most of the data where taken between elevations of 35 and 62 degrees, with a mean of 49 degrees, so we expect this bias to be less than a few percent and thus assume it is negligible.

The instrumental noise and the atmospheric noise residuals are estimated as detailed in Adam et al. (2016). In brief, we compute noise maps from the difference between two equivalent subsamples. Examples of such jackknife maps are provided in Fig. 1. We can observe that the noise is correlated over spatial scales that are larger than the beam size, in particular at 260 GHz. This can give rise to spurious signal at the scale of the subclusters, and the method we use allows us to account for such correlations as detailed below. The noise maps are normalized by the integration time spent per pixel in order to generate homogeneous noise in the observed sky patch. The homogeneous noise is then used to compute its power spectra (using the POKER software, Ponthieu et al. 2011), which is well described by a white component plus a pixel–pixel correlated contribution. We also check that the noise is Gaussian, and we do not observe deviations from gaussianity. This noise model allows us to simulate noise maps that account for the noise inhomogeneity and spatial correlations, and to compute the full noise covariance matrix from noise Monte Carlo realizations. The signal-to-noise ratio of the surface brightness maps are computed, at a given angular resolution (22 arcsec FWHM in the following), by dividing the signal by the root mean square of the noise realizations smoothed to the resolution of the signal map, and thus accounting for spatial correlations of the noise. The noise power spectrum is thus compared to the one of the signal, at the scales of the subclusters (after the filtering of the power on small scales), to infer the signal-to-noise ratio. We also estimate the cross-correlation between the 150 and 260 GHz noise using the cross spectra of the difference maps. We observe a small correlation due to the atmospheric noise residuals, which we expect because NIKA observes the sky simultaneously at 150 and 260 GHz. The corresponding correlated noise between the two bands corresponds to 0.22% of the variance at 260 GHz, and 1.25% of that at 150 GHz. However, this is negligible with respect to the cross correlation induced by the cosmic infrared background (CIB) contamination (see Sect. 3.1) and we do not account for it.

2.3. NIKA raw maps

The NIKA maps are shown in Fig. 1. They have been smoothed with a Gaussian filter to have the same effective resolution of 22 arcsec FWHM. Due to the integration time spent per pixel, the noise is relatively constant within a radius of 1.5 arcmin around the map center, and increases radially toward the edges. At the 22 arcsec effective resolution, the noise root mean square at the center of the field is 97 μJy/beam at 150 GHz and 226 μJy/beam at 260 GHz, which accounts for instrumental noise, residual atmospheric noise, and the astrophysical noise as discussed in Sect. 3.1. The root mean square noise on the flux of a point source located near the center is 0.14 mJy at 150 GHz and 0.58 mJy at 260 GHz. The positions of the four subcluster regions that we consider, A, B, C and D, are also shown in Fig. 1. Various coordinates have been provided in the literature, obtained for example by fitting dark matter profiles on strong lensing data (e.g., Limousin et al. 2012, 2016). Since we are interested in the gas distribution rather than the dark matter, we use the best-fit cored-component model coordinates of Limousin et al. (2016) as a starting point. These are adjusted slightly to match the peak distribution of the X-ray Chandra wavelet smoothed map, which is used as a prior in determining the gas distribution (see Sect. 4.2 and the top right panel of Fig. 3). These coordinates (see Table 2) are still consistent with the strong lensing analysis within the dispersion obtained from the different models by Limousin et al. (2016). Unless otherwise stated, they are used only for display purpose in this paper and the detailed choice of these centers has no impact on the results that we present.

The 150 GHz map shows the presence of a diffuse negative signal associated with the SZ effect. The signal is dominated by the tSZ signal at this frequency, but the kSZ signal is also significant (see Sect. 5). The signal is non-zero, but very small in the direction of subcluster A, in agreement with its expected low gas content (see Ma et al. 2009). The brightness in the direction of all other subclusters is significant with a peak between regions B and C, reaching −18.4σ. However, we do not observe a peaked signal in the direction of each group, but rather a diffuse signal that does not allow the clear identification of the subclusters using this map. We also report the detection of a foreground radio galaxy on the southeast edge of the cluster (see Sect. 3.2).

At 260 GHz, we observe a diffuse positive signal around region C reaching up to 4.3σ at the peak, but the surface brightness is negative in region B, with a peak at −3.7σ. The signal-to-noise ratio is insufficient to detect the signal in region A at 260 GHz. We also detect three compact sources at this frequency (see Sect. 3). Due to the relatively high astrophysical contamination observed in the 260 GHz map, its detailed interpretation requires the removal and the treatment of the point sources as will be discussed in Sect. 3. Nevertheless, we observe an excess of signal around region C with respect to the 150 GHz tSZ expectation (see Table 1 for conversions between tSZ and surface brightness at 150 and 260 GHz), and a negative signal in region B that counters the positive tSZ signal expected at this frequency. While the former could be attributed to another astrophysical signal, it is not the case for the latter, which we will interpret as a kSZ signature in Sect. 5.

3. Astrophysical contamination of the Sunyaev-Zel’dovich signal

In addition to the SZ signal from the cluster itself, we expect galactic, extragalactic, and cluster associated astrophysical signal to contaminate our data. This section describes how such contaminants are taken into account in the present paper.

3.1. Diffuse galactic emission and cosmological background

We apply the approach used in Adam et al. (2016) in order to check that the CMB and the galactic emission (we consider synchrotron, free-free and thermal dust) are negligible in the case of MACS J0717.5+3745 NIKA data. The CIB was also considered in Adam et al. (2016). Despite being the largest diffuse astrophysical noise contribution, it was shown to be negligible. The data we present here, however, are about seven times deeper, such that the CIB is no longer negligible and requires more careful investigation as we discuss below. To account for it, we produce CIB realizations, which we add to our Monte Carlo noise realizations as discussed in Sect. 2.2. The modeling and removal of the submillimeter sources that are individually detected by Herschel is treated separately as detailed in Sect. 3.3.

We consider the clustering of dusty star forming galaxies, the shot noise from dusty star forming galaxies and the shot noise from radio sources. The clustering term from radio sources is negligible (Hall et al. 2010). The clustering of dusty star forming galaxies is computed using the CIB power spectrum measured by Planck Collaboration XXX (2014) at 143 and 217 GHz. It is modeled by a one halo term and a two halo term and extrapolated to the NIKA frequencies. The shot noise originates from the discrete nature of the CIB signal and is due to unresolved sources below the detection threshold of NIKA. For dusty star forming galaxies, it is computed using the model from Béthermin et al. (2012), and we use the model of Tucci et al. (2011) for radio sources. We account for the beam smoothing and the large scale filtering implied by the NIKA processing. At the angular scales considered from the NIKA data, the shot noise dominates over the clustering by a factor of about five, and the shot noise from radio sources is negligible compared to that from dusty star forming galaxies at both NIKA frequencies. Using the results of Béthermin et al. (2013), we assume that the CIB (clustering and shot noise) is completely correlated between the two NIKA bands. At our effective angular resolution of 22 arcsec, we find that the noise root mean square of the CIB is 47 μJy/beam at 150 GHz and 130 μJy/beam at 260 GHz. This boosts the noise in our maps by 22.1% and 14.8% at 260 and 150 GHz, respectively, and induces a noise correlation between the two bands. The uncertainty on the CIB model that we use is expected to be about 20%. This translates to 7% and 5% uncertainty on the overall noise estimate at 260 and 150 GHz, respectively. Finally, we note that we expect CIB dimming due to the gravitational lensing by the cluster mass (Zemcov et al. 2013), but this effect was checked to be negligible in the case of our data.

Fig. 2 Point source model at 150 GHz (left) and 260 GHz (middle), and point source mask constructed from the 150 and 260 GHz point source models (right). The color scale of the models is the same as in Fig. 1. We also provide the labels of each source on the figure, in dashed circles for radio sources and dotted circles for submillimeter sources. The point sources subtracted maps can be seen in Fig. 7.

3.2. Radio sources

First, we consider the diffuse radio emission associated with the cluster. We extrapolate the surface brightness maps measured by NVSS (NRAO VLA Sky Survey, Condon et al. 1998) and FIRST (Faint Images of the Radio Sky at Twenty Centimeters, Becker et al. 1995) at an angular resolution of 45 arcsec and 5.4 arcsec, respectively, assuming a spectral index of −1.25 (van Weeren et al. 2009; Bonafede et al. 2009). We find that the radio emission is negligible in both cases within the NIKA bands. Moreover, this constitutes an upper limit on the radio contamination because we expect the radio spectra to be steeper at higher frequencies due to energy loss of the synchrotron-emitting electrons (Scheuer & Williams 1968). Therefore, we do not correct for the extrapolated signal and neglect its contribution.

In addition to the diffuse radio emission, a total of seven compact radio sources are identified in the NIKA field from the literature, at lower frequencies (Condon et al. 1998; LaRoque et al. 2003; Coble et al. 2007; Bonafede et al. 2009; van Weeren et al. 2009; Mroczkowski et al. 2012; Pandey-Pommier et al. 2013). Their coordinates are listed in Table 3 and they are labeled in Fig. 2. We also consider the more recent work of van Weeren et al. (2016) to add extra constraints on the flux densities of the sources. One of them, a foreground radio galaxy at z = 0.1546 (van Weeren et al. 2009), is significantly detected by NIKA at 150 and 260 GHz (source F). We assume the radio SED to follow a power law in order to model the contamination in the NIKA bands. Since we know that the field around MACS J0717.5+3745 is relatively crowded, and to avoid source confusion, we only use photometric data from which the observations are sensitive to scales smaller than 30 arcsec: FIRST, Very Large Array (VLA), Jansky Very Large Array (JVLA), Giant Metrewave Radio Telescope (GMRT), Owens Valley Radio Observatory/Berkeley Illinois Maryland Association (OVRO/BIMA) and MUSTANG/GBT. We also extract the flux density of the sources from the NIKA maps as detailed in Adam et al. (2016), except for sources HT and C4, which are located near the peak of the SZ signal. The available flux densities used to constrain the SED model, including the NIKA measurements, are listed in Table 3. In the case of source F, the spectrum is poorly described by our model over the full frequency coverage. Instead, the SED presents two distinct spectral indices at low (<2 GHz) and high (>10 GHz) frequencies. This could be due to significant free-free emission at high frequencies, in addition to the synchrotron at low frequencies, but might also indicate biases in the measurement of the flux density (obtained with different instruments) because of the source extension and location near the radio relic at low frequencies. Therefore, we assume that the power law model is valid at high frequency only and consider data at ν> 10 GHz. Despite indications of steepening, the other sources are correctly described by a power law, even if our NIKA constraints are in tension with respect to the best-fit model for C1 and C2 (see Table 3). In the case of sources C4, the limited data suggests to interpret the model prediction with caution. In addition, we emphasize that the NIKA photometry might be contaminated by submillimeter sources, as discussed in Sect. 3.3, and possibly also by SZ signal. The results flux density prediction in the NIKA bands are listed in Table 3 and are used as discussed in the following. The source model, i.e. the 150 and 260 GHz best-fit flux densities of Table 3, is given in Fig. 2, including submillimeter sources, as discussed in Sect. 3.3.

3.3. Submillimeter point sources

As mentioned in Sect. 2.3, submillimeter sources significantly contaminate our data, even at 150 GHz. We use the approach detailed in Adam et al. (2016) in order to subtract their contribution. To do so, in addition to the NIKA observation itself, we use the SPIRE (Spectral and Photometric Imaging REceiver, Griffin et al. 2010) catalog produced by Sayers et al. (2013) to clean their Bolocam images. Eleven sources are identified in the 4 arcmin × 4 arcmin region around the cluster and we report their flux densities in Table 4. The gray body model which we fit to the data, excluding the NIKA sources located near the SZ signal, allow us to better constrain the excepted flux densities in the NIKA bands (see Adam et al. 2016, for more details). The corresponding template to be subtracted from our data is provided in Fig. 2 together with the radio contribution. Most of the model predictions are fully compatible with the NIKA observations. Only SMG11 at 150 GHz differ by more than 3σ, and we note that we suspect the NIKA measurement of SMG11 to be strongly contaminated by SZ signal because of its location. The SPIRE data are also well described by the model and we do not observe significant excess in the residual between the model and the data for any of the sources.

In addition to the radio source, F, the excess seen in the NIKA 260 GHz map coincide with SMG09 (northwest of subcluster D) and SMG05 (south of subcluster D). The corresponding SED peaks at relatively low frequencies, indicating that these sources could be high redshift lensed galaxies. This is particularly true for SMG09, for which the SED peaks at about 600 GHz (500 μm). Therefore, we search for counterparts to the NIKA detected sources in strong lensing HST data set (Hubble Space Telecope, Frontier Field campaign, see Diego et al. 2015). By requiring that the projected distance is less than 5 arcsec between our identified sources and that of HST, we find two possible candidates for SMG09 at z = 4.5 (RA +07:17:31.269, Dec +37:44:41.10) at 2.0 arcsec from our coordinates, and z = 3.0 (RA +07:17:31.082, Dec +37:44:42.36) at 2.5 arcsec from our coordinates (ID25.1 and ID62.1, respectively, Diego et al. 2015). The source SMG09 is, therefore, very likely to be a high redshift submillimeter galaxy lensed by MACS J0717.5+3745, or even a blend of two galaxies. In the case of SMG05, no counterpart is found, even if we increase our threshold radius to 15 arcsec.

3.4. Diffuse submillimeter emission from the ICM

While galaxies are the major contributors to the submillimeter emission in galaxy clusters (Coppin et al. 2011), diffuse dust emission associated to the ICM could also lead to submillimeter emission (see, for example, Montier & Giard 2005; Planck Collaboration XLIII 2016, who study the overall contribution of dust in clusters). Such signal would likely be diffuse, at arcmin scales, and we expect it to be small compared to the galaxies, which give rise to emission on small scales. Using Herschel SPIRE maps toward MACS J0717.5+3745, we check that such signal is, at least, much lower than the one arising from point sources. Moreover, we do not observe any diffuse emission which would be correlated with the kSZ signal measured in Sect. 5. Therefore, the dust contribution associated with the ICM is neglected in the NIKA bands.

3.5. Construction of a mask for point sources

The spectral modeling of both radio and submillimeter point sources allows us for a good first estimate of the contamination expected in our data. However, the model uncertainties are relatively large and once subtracted to the respective channels, we observe small residuals in the case of radio sources, that are visible in particular for sources that are outside the region dominated by SZ signal. This is likely due to our modeling being too simplistic since it does not account for the steepening of the radio spectrum, as discussed in Sect. 3.2. Therefore, in addition to subtracting the point source model, we also construct a mask to avoid using pixels that we consider potentially biased when using the NIKA maps for photometry or for fitting purposes. This mask is constructed by applying a threshold (0.2 and 0.4 mJy/beam at 150 and 260 GHz, respectively) on the point source best-fit models at 150 and 260 GHz simultaneously. Such procedure is conservative and takes advantage of the high angular resolution of the NIKA data. The mask we use is given on the right panel of Fig. 2.

4. X-ray data analysis

In the work presented in this paper, we have complemented the SZ data with X-ray observations for two reasons. First, MACS J0717.5+3745 is known to be a very hot cluster (Ma et al. 2009) such that SZ relativistic corrections are expected to be fairly significant and should be taken into account by using an extra estimate of the gas temperature, available from X-ray spectroscopy. Second, X-ray data are also needed in order to disentangle the line-of-sight gas velocity from the density and the temperature (see Eq. (5)). In this section, we discuss the X-ray data reduction in terms of imaging and spectroscopy.

4.1. Data preparation

MACS J0717.5+3745 was observed by the XMM-Newton telescope using the European Photon Imaging Camera (EPIC, Turner et al. 2001; and Strüder et al. 2001) for 195 ks (obs-IDs 0672420101, 0672420201, and 0672420301) and by the Chandra Advanced CCD Imaging Spectrometer (ACIS, Garmire et al. 2003) for 250 ks in total (obs-IDs 1655, 4200, 16 235, 16 305). We processed XMM-Newton datasets using the Science Analysis System (SAS) version 15.0.0, applying the latest version of calibration files available in March 2016. We processed Chandra datasets using the Chandra Interactive Analysis of Observation (CIAO) version 4.6.5 and calibration database version 4.7. Both datasets suffer from contamination due to the high energy particle flux. To reduce this component, we removed from the XMM-Newton datasets all the events which keyword PATTERN is >4 and >13 for MOS1,2 and PN cameras, respectively. We applied Very Faint¹ (VF) mode filtering to obs-IDs 4200 and 16 305. Since our background subtraction technique for Chandra can be applied only to observations taken using the VF mode, we removed obs-IDs 1655 and 16235 from our analysis, thus reducing the effective observation time to 153 ks.

4.2. Spectroscopy analysis

Fig. 3 XMM-Newton and Chandra observations of MACS J0717.5+3745. Top left panel: exposure- and vignetting-corrected Chandra photon count image in the [0.5−2.5] keV band. The image has a resolution of 3.1′′. Top right image: wavelet de-noised map of the Chandra image shown in the top left panel. Iso-contours are shown with black solid lines to enhance substructures. The four regions in Table 2 are highlighted with black circles. Bottom left image: XMM-Newton wavelet de-noised temperature map, obtained using a threshold cut of 1σ. Contours from the wavelet map image are overlaid using white contours. Bottom right image: same as bottom left panel except that we show the temperature map computed using the Chandra dataset. All the maps shown are centered on the X-ray peak.

To remove observation intervals affected by flare episodes we followed the light-curve filtering procedures described in Pratt et al. (2007) and in the Chandra COOKBOOK² for XMM-Newton and Chandra datasets, respectively. We removed from the analysis all the intervals where the count rate exceeded 3σ from the mean value. We find no flare contamination in the Chandra dataset, so the full observation time of 153 ks was used. For the XMM-Newton dataset, the useful exposure times were 160 ks and 116 ks for MOS1,2 and PN cameras, respectively. We ran the wavdetect algorithm (Freeman et al. 2002) to identify point sources on exposure corrected images in the [0.3−2] keV and in the [0.5−2;2−8;0.5−8] keV bands for XMM-Newton and Chandra, respectively. Point source lists thus produced were inspected by eye, merged, and used as mask to remove point source contribution from the analysis. Point sources only detected by Chandra were masked from XMM-Newton dataset, using a circular region with radius 15′′ within 3′ from the aimpoint.

To perform spectral and imaging analysis we binned both datasets in sky coordinates and energy, creating an energy position photon cube for each dataset following the procedure described in Bourdin & Mazzotta (2008). To each event in the cube we assigned an effective area, which we use to correct for vignetting and exposure time, and the background noise. In the following, all the techniques described have been used on both Chandra and XMM-Newton datasets, unless otherwise stated.

The exposure- and vignetted-corrected Chandra image in the [0.5−2.5] keV band and its wavelet de-noised map are shown in the top left and right panel, respectively, of Fig. 3. The wavelet filtering is performed following the method described in Bourdin et al. (2013). Note that the resolution of the wavelet map is not fixed. It depends on the image statistics as the filtering is constrained to detect structures above 3σ at different scales. As already found by several teams (e.g., Zitrin et al. 2009; Ma et al. 2009), the cluster exhibits a very complex morphology, with at least 4 identifiable substructures. Because of the higher resolution of ACIS, we used the Chandra data as the baseline for the X-ray imaging analyses.

4.3. Background estimation

The X-ray background can be divided into a sky and an instrumental component. For the latter we used the analytical models produced by Bourdin et al. (2013) and Bartalucci et al. (2014) for XMM-Newton and Chandra, respectively. We normalized both models, extracting the events from a region free from the cluster emission in the [10−12] keV and in the [9.5−10.6] keV band for XMM-Newton and Chandra, respectively. The sky background is due to foreground Galactic emission and an additional extra-Galactic component. For the former we used two absorbed thermal models produced by the Astrophysical Plasma Emission Code (APEC, Smith et al. 2001) for which the temperatures are fixed to 0.248 keV and 0.099 keV (see Kuntz & Snowden 2000). The extra-Galactic component was modeled using an absorbed power law with γ fixed to 1.42, as proposed by Lumb et al. (2002). We determined the normalization of all the components by performing a joint fit of the spectrum extracted from a region free of source emission. We noticed the presence of a small excess in the [0.3−1] keV band in the XMM-Newton dataset, compatible with the solar wind charge exchange emission (Snowden et al. 2004). We found that one absorbed APEC model with a temperature of 0.86 keV is sufficient to reproduce this excess. We added this component in the joint fit. Once normalized, the background sky model is added to the instrumental model. Then, in all our spectral analysis procedures, we simply scale the resulting total background model by the ratio of the extraction areas.

All X-ray spectral analyses were performed following the scheme described in Bourdin & Mazzotta (2008). For a given region of interest, we add to the background model an absorbed APEC to account for the cluster emission, and fit it via χ² minimization to the exposure- and vignetting-corrected spectrum extracted from the photon cube. The absorption is fixed to the Galactic value along the line-of-sight n_H = 6.63 × 10^-20 cm^-2, as determined from the LAB survey Kalberla et al. (2005). All the models are convolved by the appropriate response matrix files, accounting for energy resolution loss due to instrumental effects. The source emission model and the sky background models are multiplied by the appropriate ancillary response file. From the fitting procedures performed in the [0.3−10] keV and [0.7−10] keV band for XMM-Newton and Chandra, respectively, we determined the APEC model temperature, abundance and normalization.

We produced temperature maps applying the wavelet filtering algorithm describe in Bourdin et al. (2004) and in Bourdin & Mazzotta (2008). Briefly, we sample the field of view in square grids where each element size increases by powers of 2. The minimum size of the meta-pixel is defined so that it contains at least 200 photons. In each meta-pixel we extract the spectrum and measure the temperature, fixing the abundance to 0.3 solar, obtaining a first set of temperature and associated fluctuation maps. We then convolve these maps with a B2-spline function evaluated at different scales to obtain the corresponding wavelet coefficients. The “de-noised” temperature map is constructed using a 1σ threshold for the wavelet coefficient maps. The de-noised temperature maps that we obtained from the XMM-Newton and Chandra datasets, are shown in the left and right bottom panel, respectively, of Fig. 3. Both temperature maps identify the same hot structure, with temperatures of the order of ~24 keV, and a cold region in the northwest sector associated with the most prominent structures. Note that the wavelet filtering is applied to maps where the element size is determined by the statistics. For this reason the XMM-Newton dataset allows a finer element sampling, due simply to the larger number of counts in this observation.

To check the validity of the values in the temperature maps, following Wang et al. (2016), we extract the spectra from the four circular regions shown in the top right image of Fig. 3, using circular radii of 15′′ and 30′′. We then measured the temperatures via our fitting procedure, which we report in Table 5 together with their 1σ errors. Because of Chandra’s shorter exposure time and poorer sensitivity, we measure the temperature only in a 30′′ circular bin aperture. The values reported for the two instruments are consistent with the values found in the respective temperature maps.

As expected, on average Chandra measures higher temperatures than XMM-Newton (see, e.g., Martino et al. 2014; and Schellenberger et al. 2015), but the values are consistent within the errors. Our results are also in agreement those of Sayers et al. (2013), even if the regions we consider here are not exactly centered on the same coordinates. In the case of subcluster C, however, the temperature that we find are about the same for both X-ray observatories. Our Chandra derived value is also lower for this subcluster than that reported in Sayers et al. (2013), but is in agreement within the errors. Because of the better resolution and significance, we use the XMM-Newton spectroscopic results as our baseline.

Fig. 4 Constraints on the kSZ and the tSZ amplitudes toward the regions associated with the four subclusters (from A to D). The black data points provide the average brightness measured by NIKA in the 22 arcsec radius regions with center coordinates given in Table 2 and shown in Fig. 1. The shaded areas give the 68% confidence level constraints on the kSZ (red), the tSZ (blue) and the sum of the two (gray). Possible point source residuals were masked when computing the average surface brightness.

5. A map of the kinetic Sunyaev-Zel’dovich signal

In the case of pure tSZ signal with small relativistic corrections, we expect the spatial distribution of the signal to be proportional at any frequency. However, this is clearly not the case for the NIKA data such that another contribution is necessary to explain our observations. Motivated by the results of Ma et al. (2009), Mroczkowski et al. (2012) and Sayers et al. (2013), we assume that the observed signal is due to the tSZ and the kSZ effects, and we use the NIKA data to disentangle the two. The main goal of this section is to produce, for the first time, a resolved map of the kSZ effect. This constitutes the main result of this paper.

5.1. Reconstruction of the kinetic Sunyaev-Zel’dovich signal

The two NIKA maps, cleaned of contamination (see Sect. 3), provide a measurement of the surface brightness, ΔI_ν, in the different regions of the cluster. The sensitivity to the tSZ and the kSZ contributions are known from the coefficients given in Table 1, and are corrected on each pixel of the sky for relativistic effects using the XMM-Newton X-ray spectroscopic temperature map in Fig. 3, under the assumption that the temperature is constant along the line-of-sight. We can therefore invert Eq. (3) in order to separate the contribution of the tSZ signal, given by ${\mathit{y}}_{\mathrm{tSZ}}\mathrm{=}\frac{\mathit{g}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{1}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\mathrm{\Delta }{\mathit{I}}_{{\mathit{\nu }}_{\mathrm{2}}}\mathrm{-}\mathit{g}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{2}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\mathrm{\Delta }{\mathit{I}}_{{\mathit{\nu }}_{\mathrm{1}}}}{{\mathit{I}}_{\mathrm{0}}\mathit{g}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{1}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\mathit{f}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{2}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathrm{0}}\mathit{g}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{2}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\mathit{f}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{1}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}}\mathit{,}$(6)and that of the kSZ signal, ${\mathit{y}}_{\mathrm{kSZ}}\mathrm{=}\frac{\mathit{f}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{1}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\mathrm{\Delta }{\mathit{I}}_{{\mathit{\nu }}_{\mathrm{2}}}\mathrm{-}\mathit{f}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{2}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\mathrm{\Delta }{\mathit{I}}_{{\mathit{\nu }}_{\mathrm{1}}}}{{\mathit{I}}_{\mathrm{0}}\mathit{f}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{1}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\mathit{g}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{2}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\mathrm{-}{\mathit{I}}_{\mathrm{0}}\mathit{f}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{2}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}\mathit{g}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{1}}\mathit{,}{\mathit{T}}_{\mathrm{e}}\mathrm{\right)}}\mathrm{·}$(7)We propagate the noise through Monte Carlo realizations as discussed in Sect. 2.2, including the CIB contribution and its induced correlation between the two NIKA frequencies. The statistical errors on both tSZ and kSZ are largely dominated by the noise at 260 GHz due to the intrinsic lower sensitivity of this band, with respect to the 150 GHz one.

5.2. Reconstruction of the kinetic signal toward subclusters

We first consider the mean SZ surface brightness measured in the individual regions provided in Table 2, i.e., within 22 arcsec radius disks. The point source mask of Fig. 2 was applied before computing the brightness, so that pixels that are potentially affected by mis-subtraction of the point sources are rejected. Following Sect. 5.1, we use these measurements to constrain the amplitude of both the tSZ and kSZ spectra. The results are presented in Fig. 4 for the four subclusters. From Eqs. (3) and (5), we can see that a positive kSZ amplitude corresponds to a negative gas line-of-sight velocity, which means that the cluster is moving toward the observer, with respect to the CMB reference frame. The signal-to-noise ratio is too low in region A and both tSZ and kSZ amplitudes are compatible with zero. However, at 150 GHz the SZ signal is significantly detected, and we assume this to be dominated by tSZ. In region B, we detect tSZ and kSZ signal at 2.6σ and −4.6σ, respectively. Moreover, we find that the kSZ induced brightness is slightly more negative than that of the tSZ, at the scales probed by NIKA. The kSZ signal is negative, which explains the decrement observed in this region at 260 GHz with respect to the 150 GHz expectation. In region C, we significantly detect tSZ signal and find indication for kSZ signal at 7.8σ and 2.1σ, respectively, but the situation is reversed with respect to region B and the kSZ signal is positive. The tSZ signal strength is very high, in agreement with the expectation from the merger scenario, responsible for gas heating in this region. This is also in agreement with our X-ray derived temperature map. Finally, tSZ signal is detected in region D at 4.3σ and the kSZ significance is 0.6σ.

5.3. Mapping of the kSZ and tSZ signal

The results of Sect. 5.2 provide quantitative illustration of our constraints, but they directly depend on the choice of the region coordinates and their aperture size. Therefore, we now take advantage of the resolved nature of the NIKA data and consider the full surface brightness maps (in addition to the XMM-Newton temperature map for relativistic corrections) and we use them to compute maps of the tSZ and kSZ amplitude, i.e., maps of y_tSZ and y_kSZ, that are independent of our baseline regions. Before combining the two data sets, we convolve the NIKA maps to the same angular resolution, i.e., 22 arcsec. The resulting maps are provided in Fig. 5. The residuals of point sources are not considered here, but the obtained maps are compared to the mask to check that the kSZ signal we detect is not affected by point sources. The coordinates of the point sources identified in Sects. 3.2 and 3.3 are also represented in Fig. 5. While these regions are potentially biased, they remain local and do not spatially correlate with the main significant kSZ structure that we observe. We observe a tSZ peak coincident with the region of subcluster C. The peak signal-to-noise ratio reaches 7.3σ, which is more than twice smaller than on the 150 GHz map due to the noise introduced by the 260 GHz data. The signal is more diffuse in the other regions but it covers the cluster extent on arcmin scale. The kSZ map shows two distinct peaks of opposite signs and similar amplitudes, reaching −5.1 and + 3.4σ, respectively. The negative peak is almost coincident with region B and the positive peak is located at the northeast edge of region C. No kSZ signal is significantly detected near subclusters A or D. Similarly, we do not observe any diffuse large scale kSZ signal that would result from the overall motion of the cluster and subclusters.

Fig. 5 Map of the thermal SZ effect, ytSZ (left), and the kinetic SZ effect, ykSZ (right). Gray contours are multiples of 1σ, starting at ±2σ. The map effective resolution, 22 arcsec, is shown as a white circle on the bottom left corner. Subcluster regions are represented in gray. The positions of the identified point sources (see Sect. 3), from which mis-subtraction residuals remain a potential local contaminant to the signal, are represented in light gray.

5.4. Systematic effects

In addition to the noise and astrophysical contaminants, which we account for as discussed in Sect. 3, our measurements are affected by three main systematic effects, which we discuss here.

The first is the absolute calibration uncertainties of the NIKA data. A change in the absolute calibration directly impacts the kSZ measurement since it changes the relative weights of the two NIKA bands. Assuming a typical gas temperature of 10 keV, Eqs. (6) and (7) can be expressed as y_tSZ ∝ −0.38 × ΔI_ν2 + 0.62 × ΔI_ν1 and y_kSZ ∝ 0.20 × ΔI_ν2 + 0.80 × ΔI_ν1, where the numerical coefficients indicates the fraction of the maps at 260 GHz and 150 GHz that contribute to the tSZ and kSZ maps. We can see that changing our 150 GHz calibration by 7% (the calibration uncertainty), leads to a change of 2.7% for the tSZ map, and 1.4% for the kSZ map. Similarly, changing our 260 GHz calibration by 12% leads to a change of 7.4% for the tSZ and 9.6% for the kSZ map. While the absolute calibration uncertainty does not significantly affect the significance of detection, it can slightly change the relative strength of the signal in the different regions. If we focus on region B and C, we notice that the 150 GHz surface brightness is negative in both regions, while it changes sign at 260 GHz, being negative in region B. On the other hand, the two linear coefficients that are used to compute the kSZ map from the primary NIKA maps (0.20 and 0.80 as written above) are both positive. Therefore, if the calibration changes in the same direction at both frequencies, the changes will accumulate in region B, and they will oppose each other in region C. Similarly, if the calibration changes are opposite in sign, they will oppose each other in region B, and they will accumulate in region C. Therefore, summing the calibration uncertainties, the kSZ signal in region B can increase (or decrease) by 11.0% while it increases (decrease) only by 8.2% in region C. Symmetrically, the same is true when inverting B and C. Note that we apply a simple sum of the calibration errors instead of a quadratic sum because they are likely to be strongly correlated. This effect will propagate linearly to the constraint on the gas line-of-sight velocity in Sect. 6, but it is small compared to the statistical uncertainties. Similar arguments apply to the tSZ map.

The second systematic effect is due to the uncertainty in the derived X-ray temperature, and the assumption that the temperature is constant along the line-of-sight. In order to estimate how changes in the temperature affect our results, we reproduce the kSZ map applying ±25% change in the normalization of the temperature map (i.e., about twice the typical difference between XMM-Newton and Chandra). This allows us to test the systematic effect responsible for the temperature difference between Chandra and XMM-Newton. We also add Gaussian noise of amplitude 3 keV at a resolution of 22 arcsec. Due to the fact that the relativistic corrections are close to proportional in the two NIKA bands, the changes in the temperature map lead to only small changes in the kSZ significance, i.e. less than 0.15σ for the signal peaks. The absolute amplitude of the signal, however, is affected by up to 0.4σ changes at the positive peak, which correspond to hot gas in region C, where relativistic effects are important.

The third systematic effect is due to the large angular scale filtering that affects the NIKA data. The filtering is the same in the two bands, and as the quantities y_tSZ and y_kSZ are a linear combination of the two NIKA maps, they are affected by the same filtering. Therefore, the tSZ and kSZ reconstructed signals are smoothly apodized at scales larger than ~2 arcmin. The zero level of the kSZ map is directly given by the respective zero levels of the 150 and 260 GHz map, being set to zero in the external regions of the map, as detailed in Sect. 2.2. While an arbitrary change in the brightness zero levels could change the kSZ signal zero level, the peak to peak difference between the positive and negative observed signal will remain unchanged. In fact, the kSZ signal we observe is much more compact than that of the tSZ effect, so we expect filtering and zero level effects to be less important for kSZ observations. The transfer function is nonetheless accounted for when fitting our data with a model as detailed in Sect. 6.

5.5. Comparison to maps at other wavelengths

Fig. 6 Comparison between the kSZ signal and other probes of the cluster components, shown on saturated linear scales. The astrometries are the same as those in Fig. 5, and the green contours reproduce the kSZ signal-to-noise ratio, given as a solid line if positive (+2 and +3σ) and dashed line if negative (−2, −3, −4 and −5σ). Our reference regions are provided as gray circles. The positions of B1 and B2 are also shown as gray diamonds on the left image. Left: Chandra X-ray photon count image, smoothed to 7 arcsec FWHM. Middle: HST optical combined image from the CLASH data (Postman et al. 2012). Right: lensing mass (κ) from Zitrin et al. (2011).

In Fig. 6, we compare the kSZ map obtained with NIKA to maps at other wavelengths, taken as tracers of the different cluster components. We use the Chandra X-ray photon counts as a proxy for the gas distribution (∝${\mathit{n}}_{\mathrm{e}}^{\mathrm{2}}$), HST optical imaging (CLASH data, Postman et al. 2012) to identify the galaxy spatial distribution, and the strong lensing mass reconstruction to probe the dark matter (Zitrin et al. 2011).

The X-ray image shows two main peaks in region B. These two peaks are also seen in SZ observations by MUSTANG at 90 GHz with a 13 arcsec effective angular resolution (Mroczkowski et al. 2012). The most significant peak is located on the east side of the region (hereafter B1), with the second on the west side (hereafter B2), the two being separated by about 20 arcsec. The group of galaxies associated with subcluster B is on average more coincident with B2, even if galaxies are also observed around B1. The dark matter map follows the distribution of the galaxies, and is therefore better aligned with B2. The negative kSZ signal coincides well with B1, but is significantly offset from B2. The morphology of clump B2 is that of an arc (or a tail) around the dense core, similar to what is observed in mergers on the plane of the sky (e.g., Clowe et al. 2006, for the bullet cluster). This could suggest, if this is the case, that B2 is an independent subcluster (with respect to B1) moving mostly perpendicular to the line-of-sight, and therefore, not responsible for the kSZ signal. The signal associated to B2 could also be due to the stripped gas of the main core B1, which would therefore move slower than B1 and is not significantly detected in kSZ. In both cases, our kSZ data, together with X-ray imaging, suggest that subcluster B is a complex object in itself. The main structure, B1, is rapidly moving along the line-of-sight, while the second one, B2, is moving either mostly perpendicular to the line-of-sight, or in the same direction than B1, but slower.

The positive kSZ signal coincides well with the northern half of the galaxies associated with subcluster C, while the other half of the galaxies are coincident with no kSZ detection. The X-ray peak is centered on the southern part of the galaxy group and the offset between the kSZ peak and the X-ray peak is about 20 arcsec. As this subcluster is highly disturbed, we can expect strong inhomogeneities in the gas velocity distribution, but any further interpretation would require higher signal-to-noise ratio kSZ imaging. Nonetheless, it is worth noting that the galaxy velocity dispersion within this subcluster is very high ($\mathrm{176}{\mathrm{1}}_{-607}^{\mathrm{+}\mathrm{234}}$ km s^-1 for 10 redshifts, see Ma et al. 2009), which could also suggest that the underlying distribution is multimodal.

No kSZ signal is detected around subcluster D but we observe an offset between the gas and the brightest galaxies and dark matter, the gas density peak being about 25 arcsec shifted to the south with respect to the other component. Subcluster A is not clearly identified in the X-ray data, but the galaxy and dark matter locations are consistent. Finally, the foreground radio galaxy, F, is visible in X-ray and in optical.

6. Constraints on the gas line-of-sight velocity distribution of MACS J0717.5+3745

The detection of the kSZ effect implies the detection of gas motion, and opens up the possibility of studying the velocity distribution of the gas contained within the cluster, i.e. gas substructure. We detect a positive line-of-sight gas velocity in the direction of subcluster B, i.e., the gas is moving away from us in the CMB reference frame, and negative gas line-of-sight velocity toward subcluster C, i.e. the gas is moving toward us. However, measuring the absolute line-of-sight velocity from our kSZ map requires us disentangling it from the integrated line-of-sight electron density (see Eqs. (3) and (5)), accounting for projection effects, and this cannot be done from our SZ observations alone. In this section, we model the cluster gas distribution in terms of electron density and gas line-of-sight velocity, and we fit our data using extra constraints available from the X-ray spectroscopic temperature. Unlike the results presented in Sect. 5, the results we detail in this section are therefore dependent on the model we use to describe the cluster. In particular, we stress that any mis-modeling of the gas density, which is intrinsically subject to strong assumptions for such a complex object, will be reflected in the optical depth, that is itself degenerate with the gas line-of-sight velocity.

6.1. Physical modeling of the gas distribution of MACS J0717.5+3745

As discussed in Sect. 1, MACS J0717.5+3745 is known to be multimodal with four main subclusters. Therefore, we model the cluster with a set of four spherically symmetric subclusters describing the gas distribution associated to each subcluster. We model each electron density profile with a β-model (Cavaliere & Fusco-Femiano 1978), ${\mathit{n}}_{\mathrm{e}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{n}}_{\mathit{e}\mathrm{0}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}{\left[\mathrm{1}\mathrm{+}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{c}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}}\right)}^{\mathrm{2}}\right]}^{\mathrm{-}\mathrm{3}{\mathit{\beta }}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathit{/}\mathrm{2}}\mathit{,}$(8)where i labels each subcluster, n_e0 is the central electron density, r_c the core radius and β is the outer slope of the profile. We also assume that the subclusters are isothermal, with a gas temperature ${\mathit{T}}_{\mathrm{x}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}$ estimated from XMM-Newton within 30 arcsec radius apertures, as reported in Table 5. The line-of-sight gas velocity of each subcluster is also assumed to be a constant, i.e. each subcluster has its own bulk velocity. The choice of the β-model is motivated by its simplicity and by the limited number of parameters needed to describe the gas density. Moreover, it is not guaranteed that a more complex radial profile would provide a better description of the data since it is the spherical symmetry assumption that is the major limitation of our model.

The NIKA surface brightness maps are the primary observables to which we compare our model, and they are modeled as $\begin{array}{ccc}\frac{\mathrm{\Delta }{\mathit{I}}_{\mathit{\nu }}^{\mathrm{model}}}{{\mathit{I}}_{\mathrm{0}}}& \mathrm{=}& {\mathit{\sigma }}_{\mathrm{T}}\sum _{\mathit{i}}{\mathit{f}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{x}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{\right)}\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{x}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\mathrm{\int }{\mathit{n}}_{\mathrm{e}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{d}\mathit{l}\\ & & \end{array}$(9)where ${\mathit{T}}_{\mathrm{x}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}$ and ${\mathit{v}}_{\mathit{z}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}$ are scalar quantities. The line-of-sight integrations are computed analytically for each sky pixel as ${\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{+}\mathrm{\infty }}{\mathit{n}}_{\mathrm{e}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{d}\mathit{l}\mathrm{=}\sqrt{\mathit{\pi }}{\mathit{n}}_{\mathit{e}\mathrm{0}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}{\mathit{r}}_{\mathrm{c}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\frac{\mathrm{\Gamma }\mathrm{(}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{\beta }}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{-}{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{)}}}{\mathrm{\Gamma }\mathrm{\left(}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{\beta }}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{\right)}}{\left[\mathrm{1}\mathrm{+}{\left(\frac{{\mathit{R}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}}{{\mathit{r}}_{\mathrm{c}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}}\right)}^{\mathrm{2}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}\frac{\mathrm{3}{\mathit{\beta }}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}}{\mathrm{2}}}\mathit{,}$(10)where R is the projected radius from each subcluster centers and the quantity Γ represents the Gamma function. Since the SZ surface brightness depends linearly on the electron density, the contribution of each subcluster is summed regardless of cross-cluster terms that would appear otherwise (e.g., for the X-ray surface brightness, which depends on ${\mathit{n}}_{\mathrm{e}}^{\mathrm{2}}$). This allows us to ignore the exact coordinate of each subcluster along the line-of-sight. The projected location of the subclusters centers, however, are parameters of our model.

6.2. Fitting algorithm

The fitting approach we use consists of processing test simulations to incorporate the same observational effects that affect our data (see Sect. 2.2), in order to predict the observable signal that we can directly compare to the SZ surface brightness maps. The SZ maps, obtained for each point of the parameter space that is explored are convolved to the beam and the transfer function of the NIKA processing. We also account for the overall zero level of the maps individually using two nuisance parameters, because they are not constrained by the NIKA observations³. The zero level is, however, not correlated to the constraint on the velocity, corresponding to the fact that the kSZ map does not show any signal on large scales. Indeed, while the value of the χ² would change with the zero level, the local peaks observed in the kSZ signal of Fig. 5, which are related to the constraints on the velocity we obtain, cannot be accounted for by any offset in the maps. The parameter space is sampled using Monte Carlo Markov Chains (MCMC). The Metropolis-Hasting algorithm (e.g., Chib & Greenberg 1995) is used to define the evolution of the chains, by using the Gaussian log likelihood computed as $\mathrm{ℒ}\mathrm{\propto }\sum _{\mathit{m,n}}{\mathrm{(}\mathrm{\Delta }{\mathit{I}}_{\mathrm{NIKA}}^{\mathrm{data}}\mathrm{-}\mathrm{\Delta }{{\mathit{I}}_{\mathrm{NIKA}}^{\mathrm{model}}}^{\mathrm{)}}}_{\mathit{m}}{\mathrm{(}{{\mathit{C}}_{\mathrm{NIKA}}^{-1}}^{\mathrm{)}}}_{\mathit{mn}}{\mathrm{(}\mathrm{\Delta }{\mathit{I}}_{\mathrm{NIKA}}^{\mathrm{data}}\mathrm{-}\mathrm{\Delta }{{\mathit{I}}_{\mathrm{NIKA}}^{\mathrm{model}}}^{\mathrm{)}}}_{\mathit{n}}\mathit{,}$(11)where m and n label sky pixels and subscript NIKA stands for the joint 150 GHz and the 260 GHz data. The quantity C_NIKA is the full noise covariance matrix of the NIKA surface brightness maps, including the noise correlation between the two bands (see Sect. 3.1). Pixels that are potentially contaminated by point source residuals are rejected using the mask of Fig. 2, so that they do not contribute to the likelihood, and therefore, do not affect the results.

Fig. 7 Comparison between the data and the best-fit model. From left to right, the figure shows the input data subtracted from point sources, the best-fit model F1, the residual F1, the best-fit model F2, and the residual F2. The top row provides the 260 GHz data and the bottom row the 150 GHz data. Contours are multiples of 2σ at 150 GHz and 1σ at 260 GHz, starting at ±2σ.

Our model includes a total of 26 parameters (4 × 3 for the density, 4 for the line-of-sight velocity, 2 × 4 for the coordinates, and 2 for the zero level of the maps). Even if the model describing each subcluster is relatively simple, the overall complexity of the cluster, and in particular the degeneracies between the characteristic radius, slope parameters, and subcluster centroids (which are not well determined), does not allow us to let all parameters free. In particular, the slopes of the density profile are not well constrained by our data. Therefore they are varied but limited to the range [1/3, 3]⁴, in order to marginalize over their uncertainties and thereby avoid diverging chains in the MCMC. We also apply a 10 arcsec standard deviation Gaussian prior on the coordinates of the subclusters, centered on the coordinates given in Table 2. The model is symmetric under the permutation of subclusters, and this ensures that we avoid exchanges between the different subclusters during the evolution of the chains, but we check a posteriori that the best-fit coordinates are consistent with our prior values.

In the case of subcluster A, the observed signal is almost consistent with noise at the 3σ level at 150 GHz and is not detected at 260 GHz. The velocity, which is fully degenerate with the gas density, is therefore constrained at the 3σ level. This prevents the MCMC sampling of the posterior likelihood from converging. Therefore, we do not fit for the velocity of subcluster A and set ${\mathit{v}}_{\mathit{z}}^{\mathrm{\left(}\mathrm{A}\mathrm{\right)}}\mathrm{=}\mathrm{0}$.

The convergence of the chains are verified using the Gelman–Rubin convergence criteria (Gelman & Rubin 1992). After convergence, we remove a burn-in phase and account for the correlation length of the chains. Finally, the chains histograms provide the probability density function in parameter space.

6.3. X-ray prior on the gas density

The fitting algorithm discussed in Sect. 6.2 does not use X-ray imaging information because it suffers from projection effects that are not accounted for in our model. Our model would require additional parameters to predict the expected X-ray surface brightness because, unlike for SZ images, X-ray images requires the knowledge of the relative coordinates of the subclusters along the line-of-sight. This is due to cross-terms between subclusters, which appear when squaring the density, and would be responsible for additional emission on large scales. However, NIKA SZ data alone are subject to strong degeneracies between the line-of-sight velocity and the optical depth, which can be broken with additional information from X-rays. This can significantly improve the constraint on the gas velocity, but at the cost of requiring extra assumptions on the line-of-sight geometry of the cluster.

In the following, we therefore consider two cases when fitting the model of Sect. 6.1 to our data.

• 1.

  As our baseline, hereafter F1, we do not consider X-rayimaging. The fit relies only on SZ imaging and X-ray estimates ofthe gas temperature (see Table 5). In this case, wedo not use all the available information to constrain the velocity,but we are not significantly affected by assumptions concerningthe line-of-sight distribution of the gas.

• 2.

  The second fit, hereafter F2, makes use of the extra information from X-ray imaging. In this case, we consider the deprojected XMM-Newton density profile centered on the X-ray peak (extracted as in Adam et al. 2016), i.e. within a few arcsec of region B. We fit the profile with a β-model and use the constraint obtained on n_e0 (the parameter that is the most degenerate with v_z) as a prior in the fit described in Sect. 6.2. Such an approach is only possible toward subcluster B because any density profile extraction in the other regions, i.e. away from the X-ray peak, would be strongly affected by deviations from spherical symmetry.

While fit F1 is conservative, fit F2 allows us to break the v_z – τ degeneracy, and to check the consistency between our SZ-based constraint and X-ray data.

Fig. 8 Constraints of the model on the ICM distribution in the plane line-of-sight velocity – central optical depth, ${\mathit{v}}_{\mathit{z}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}$ – ${\mathit{\tau }}_{\mathrm{0}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}$, for subclusters B–D. The purple, dark blue and light blue are constraints contours at 68, 95 and 99% confidence limit in the case of our baseline fit (F1, without X-ray imaging prior on the gas density). For subcluster B, the red contours are similar to the blue ones in the case of the fit F2 (with the X-ray imaging prior). The marginalized probability density distribution are also given by the histograms in the case of F1. The dashed line give the zero location of the velocity axis. As discussed in the text, given the limitation of the models to describe the cluster, we stress that the error contours do not necessarily reflect the true underlying uncertainties.

6.4. Constraints on the velocity

The results of the MCMC fit F1 and F2 are shown in Fig. 7, including the input data, the maximum likelihood model, and the residual between the two. Despite the complexity of the cluster, a simple multi-component β-model is able to reproduce the SZ surface brightness data fairly well in both cases.

The 150 GHz residual map is consistent with noise over most of the cluster extent for both F1 and F2, but we observe several features exceeding 2σ. In the southwest region, an excess larger than 4σ is observed for both F1 and F2, indicating significant deviation from spherical symmetry for subcluster D. We also observe a significant negative excess in the direction of the foreground galaxy F, indicating that our radio model predicts a larger flux than measured, perhaps because we have neglected the possible steepening of its spectrum at high frequencies. The corresponding pixels are masked in the fit (see Sect. 3.5).

The 260 GHz residual maps are on overall consistent with noise because of the lower signal-to-noise ratio at this frequency, but we notice that the shape of our model for subcluster B is too extended (or not peaked enough) compared to the 260 GHz data. It results in a positive residual around the kSZ peak, and a negative residual at the peak location, reaching −2.3σ for F1 and −3.1σ for F2. This could be due to the fact that subcluster B is itself made of two substructures, as discussed previously. Both B1 and B2 give rise to tSZ signal so that the signal at 150 GHz is relatively extended, while the kSZ, dominant at 260 GHz, is only seen toward B2, which would lead to a more compact signal than the one at 150 GHz. Alternatively, the constant temperature assumption could be invalid. If subcluster B is a cool-core system, as expected from Ma et al. (2009), the pressure profile would be more extended than the density profile. The best-fit density profile would therefore be too extended because it is partly driven by the tSZ signal, i.e. the pressure. This is what we observe, and it is certainly one of the limitations of our model. The fact that the residual is larger for model F2 in the direction of subcluster B could also indicate that all the gas, as measured from the X-ray density, does not contribute to the kSZ signal we observe. Indeed, a lower amount of moving gas is preferred to explain simultaneously our 150 and 260 GHz data, assuming the given XMM-Newton temperature (see also Fig. 8). Therefore, our data favor (with low significance though) the hypothesis in which subcluster B is made of two separate subclusters moving with different velocities along the line-of-sight.

The MCMC chains are used to compute the optical depth at the center of each subcluster, expressed as ${\mathit{\tau }}_{\mathrm{0}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{=}\sqrt{\mathit{\pi }}{\mathit{\sigma }}_{\mathrm{T}}{\mathit{n}}_{\mathit{e}\mathrm{0}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}{\mathit{r}}_{\mathrm{c}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\frac{\mathrm{\Gamma }\mathrm{(}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{\beta }}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{-}{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{)}}}{\mathrm{\Gamma }\mathrm{\left(}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{\beta }}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{\right)}}\mathrm{·}$(12)This is the quantity that best reflects the degeneracy between the ICM density distribution and the line-of-sight velocity (see Eq. (5)). Figure 8 provides the posterior probability density function in the line-of-sight velocity versus central optical depth plane for each subcluster B–D. The one dimensional marginalized distributions are also provided as histograms.

As expected, in Fig. 8 we observe a strong degeneracy between ${\mathit{\tau }}_{\mathrm{0}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}$ and ${\mathit{v}}_{\mathit{z}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}$ for all the subclusters in our baseline fit F1 (in blue), even though it is slightly weaker for subcluster C. This degeneracy is broken in the case of F2 when using priors (subcluster B only) on the X-ray density, and the constraint are significantly improved (in red). While the line-of-sight velocity is fully compatible with zero for subcluster D, subcluster B and C show velocity distributions that are significantly different from zero.

The constraints from F1 and F2 are consistent, and the residual with respect to the best-fit model is not significantly degraded in case of F2 (see Fig. 7). Nevertheless, we observe a small tension (less than 2σ) in the case of subcluster B between the fits F1 and F2, which could be due to subcluster B being made of two subcomponents with different line-of-sight velocity, as discussed in the case of Fig. 7.

In addition to the central optical depths, our model allows us to compute the optical depth map of the cluster. This map is not independent of the best-fit velocities because of the ${\mathit{\tau }}_{\mathrm{0}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{-}{\mathit{v}}_{\mathit{z}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}$ degeneracies, but it does not depend significantly on the way the velocity is modeled because the density model is mostly driven by the tSZ signal. We use this τ model together with our kSZ map to obtain a gas line-of-sight velocity map via Eq. (5), which we use to identify structures in the velocity and assess the validity of our assumption about the velocity being constant in each region. The optical depth map and the line-of-sight velocity map are shown in Fig. 9 for both fit F1 and F2. The kSZ map is subject to statistical uncertainties, while the optical depth model imposes systematic uncertainties not accounted for in the velocity map. Therefore, the error contours provide the signal-to-noise ratio but they do not reflect the true uncertainties of the map because it depends on the model (the signal-to-noise contours are thus identical for the two models, and to that of the kSZ map of Fig. 5). In particular, the constant temperature assumption can affect the shape and amplitude of the subclusters in the optical depth map. As the uncertainties of the velocity map increase with decreasing optical depth, we mask pixels for which the noise rms is boosted by a factor larger than three, with respect to the minimum noise rms of the map.

In the case of our baseline fit, F1 (Fig. 9, left panels), the optical depth presents a smooth structure which is maximum in region C. This is counter intuitive with respect to the X-ray imaging where the peak is clearly located toward subcluster B. The X-ray imaging being sensitive to ${\mathit{n}}_{\mathrm{e}}^{\mathrm{2}}$, this could indicate that subcluster B is highly peaked while subcluster C is more diffuse and extends further in space. Our best-fit model F1 being limited by degeneracies, we check this result by computing the best-fit model F2 (right panels). In this case, the optical depth map increases toward subcluster B by about 40%, and its value between subcluster B–D stays almost constant. Nonetheless, we do not observe a peak toward subcluster B and our model suggests that it is more compact than subclusters C and D.

On the velocity map, we observe a dipolar distribution associated to subclusters B and C. The map is mostly negative between regions B and D but it is not significant. In the case of the fit F2, the velocity toward subcluster B reduces by about 40%. The dipolar structure becomes therefore more symmetrical in terms of amplitude. We also note that unlike the optical depth model, the velocity map is convolved by the NIKA transfer function similarly to the kSZ map. We report our line-of-sight velocities from our two approaches, alongside velocities from other studies in Table 6.

Fig. 9 Top: map of the best-fit optical depth model. Bottom: projected gas line-of-sight velocity map. The resolution of the map is 22 arcsec FWHM as provided by the white circle in the bottom left corner. The left panels correspond to fit F1 and the right panels to fit F2. Regions where the rms of the noise is larger than 3 times the minimum value have been masked for display purpose. Contours show the detection level of the velocity (thus proportional to that of the kSZ map of Fig. 5) and are spaced by 1σ, starting at ±2σ. We stress that they do not reflect the true uncertainty on the velocity because they do not account for the statistical and systematic effect of the optical depth models. The field shown in here is smaller than the other ones displayed in this paper.

In summary, the best-fit residuals indicates that subcluster B is more complex than a single isothermal β-model. Our constraints indicate large velocities, but they are limited by degeneracies with the optical depth. They can be broken using X-ray imaging at the cost of extra assumptions. We finally extract a gas velocity map, showing a dipolar structure, using the kSZ map and our optical depth model.

6.5. Comparison to previous results and discussion

The first kSZ detection toward MACS J0717.5+3745 was made by Bolocam using 140 and 268 GHz observation (Sayers et al. 2013), following constraints on the velocity of galaxies by Ma et al. (2009), as discussed in Sect. 1. Table 6 summaries velocity measurements performed toward MACS J0717.5+3745.

The raw NIKA and Bolocam maps are difficult to directly compare because they probe different angular scales as a result of beam smoothing and large scale filtering. However, the best-fit model of Sayers et al. (2013) and our best-fit model are deconvolved from the respective transfer functions and can be quantitatively compared. While we consider a physical parametric description of the ICM, Sayers et al. (2013) use X-ray data to predict a tSZ template, to which they add a β-model component to include the kSZ signal. Each model has its advantages and disadvantages and both are limited in the case of a complex system like MACS J0717.5+3745. Our model is clearly an idealization of the cluster, but it allows the building of a self-consistent description of the ICM. The model of Sayers et al. (2013), on the other hand, does not assume any parametrization of the ICM for the tSZ signal, but it does use a β-model for the kSZ, and the two are not necessarily related to the same underlying gas density distribution. The tSZ template built by Sayers et al. (2013) is based on X-ray data and is therefore likely to better reproduce the projected geometry of the cluster than our model. However, it requires the conversion of ${}^{\mathrm{\int }}\mathit{n}_{\mathrm{e}}^{\mathrm{2}}\mathrm{d}\mathit{l}$ to ${}^{\mathrm{\int }}\mathit{n}_{\mathrm{e}}\mathrm{d}\mathit{l}$, which is done via an effective line-of-sight extent of the ICM, ℓ_eff. The latter is left as a free parameter, but is assumed to be a constant for the entire cluster extent. This is a strong assumption and our baseline model is not sensitive to such an effect. In fact, ℓ_eff is directly related to the clumping of the gas, $\mathrm{⟨}{\mathit{n}}_{\mathrm{e}}^{\mathrm{2}}\mathrm{⟩}$/⟨ n_e ⟩ ², which is expected to be of the order of 1.4–1.6 for disturbed clusters (e.g., Zhuravleva et al. 2013). At the scales probed by NIKA, this would lead to local inaccuracies in the tSZ template (∝$\sqrt{{\mathit{\ell }}_{\mathrm{eff}}}$) that could be mis-interpreted as kSZ signal. Sayers et al. (2013) assume that the temperature is constant along the line-of-sight, while we assume isothermal subclusters, both being based on X-ray spectroscopy. Finally, Sayers et al. (2013) use fixed coordinates for their β-model center and check a posteriori the effect of changing them within reasonable uncertainties. Due to the higher angular resolution of NIKA, we fully marginalize over the subcluster coordinates in our fit.

Sayers et al. (2013) also consider a nearly model independent approach, in which they use the direct SZ flux integration of their maps in the regions B and C (corrected for a model-based zero level), together with the isothermal assumption. Using this direct approach (referenced as “direct” in Table 6), they obtain consistent results with respect to their model dependent approach (referenced as “model” in Table 6), but with reduced significance.

Despite the fact that the cluster is modeled differently, our constraints on the gas line-of-sight velocity are compatible with the results from Sayers et al. (2013). In the case of subcluster C, however, the two results are in mild tension (≳2σ). While this could easily be attributed to mis-modeling, it could also be due to differences in the observations. Indeed, our data are sensitive to scales as low as ~20 arcsec, and the constraint in region C is driven by the positive peak seen at 260 GHz at small scales. Due to Bolocam’s larger beams, the corresponding scales are likely to be diluted over the cluster extent, leading to velocity averaged over a larger extension, and hence smaller. We also use our model to infer the integrated Compton parameter in the regions considered by Sayers et al. (2013). We compute $\mathit{Y}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}\mathrm{=}{}^{\mathrm{\int }}\mathit{y}_{\mathrm{tSZ}}\mathrm{d\Omega }$, where Ω is taken as a 60 arcsec diameter circular region centered on subclusters B and C, using the coordinates of Sayers et al. (2013). We obtain ⁽0.173 ± 0.024⁾ × 10^-3 arcmin² for region B ((0.221 ± 0.015) × 10^-3 arcmin² in the case of F2), and ⁽0.242 ± 0.024⁾ × 10^-3 arcmin² for region C. These values are in agreement within 1σ (2σ in the case of F2) of the values measured by Sayers et al. (2013).

Our constraints are also roughly compatible with the optical spectroscopic measurement of Ma et al. (2009), within the error bars. However, as discussed in Sect. 5.5, our data hint at substructure within subclusters B and C themselves. As a possible consequence, our measurement in the direction of subcluster C, which is dominated by the northern part of the region, indicates higher velocities than those measured in optical wavelengths. Additional optical spectroscopy would be necessary to better address the comparison between the galaxy velocities and the gas velocity.

The X-ray photon count image presents two main peaks, associated with subclusters B and C, which is elongated to subcluster D (Fig. 3). The presence of hot adiabatically compressed gas between B and C/D suggests a pre-merger scenario between the two main structures. While the morphology of the temperature, forming a bar (see Fig. 3), indicates that the merger is likely oriented on the plane of the sky, our kSZ data and the constraints we set on the gas velocity of each subcluster shows that there is a strong line-of-sight component to the merger axis. Our best-fit velocities depend on the gas density model and the error contours are large due to degeneracies, but our detection level is stable, and our data allow us to exclude null velocity at 5.1σ for subcluster B, and 3.4σ for subcluster C. Subcluster B is likely to be falling onto C (the most massive of the subclusters, e.g., Limousin et al. 2016). We note that the projected merger axis (northwest/southeast) coincides with the filament observed on the southeast of subcluster C (e.g., Medezinski et al. 2013). If the filament axis has a non-zero line-of-sight component, it could contribute to an explanation for the high velocity of B, via its gravitational interaction with the subcluster. While the main merger scenario supported by our data is relatively clear, the detail of each subcluster remains difficult to interpret because of complex structure within the subclusters themselves.

The root mean square of cluster peculiar velocities is expected to be of the order of 250 km s^-1 (e.g., Hernández-Monteagudo & Sunyaev 2010). The velocities we measure in MACS J0717.5+3745 deviate significantly from the bulk velocity flow. However, the free fall of massive subclusters onto one another from large distances is expected to lead to relative velocities of several thousands of km s^-1 (Sarazin 2002). The velocities measured in MACS J0717.5+3745 are high, but in line with expectations. Our results suggest that other well known mergers could be ideal targets for kinematic studies using kSZ mapping.

7. Summary and conclusions

The merging cluster, MACS J0717.5+3745, was observed using the NIKA camera at the IRAM 30 m telescope. The diffuse SZ signal is detected in the two NIKA bands at 150 and 260 GHz. We report the detection of radio and submillimeter point sources that contaminates our data. Using radio data from the literature and submillimeter data from SPIRE, in addition to the NIKA constraints themselves, we model the contaminant sources and extrapolate them to the NIKA bands to clean our data. One of the submillimeter sources detected with NIKA coincides with a high redshift lensed galaxy detected by HST. If confirmed, this would illustrate the strength of NIKA SZ observations in detecting high redshift sources lensed by clusters of galaxies by combining NIKA and Herschel data (see also Adam et al. 2015, 2016). In addition to point source removal, we take advantage of the high resolution of the NIKA observations to mask possible residuals.

To effectively use our data, we account for the different sources of noise including uncorrelated instrumental noise, atmospheric and electronic correlated noise residuals and a significant contribution from the CIB. While the non-astrophysical noise reduces by increasing the integration time, this is not the case for the CIB. Since our kSZ mapping is limited by statistics, itself dominated by the noise in our 260 GHz band, the CIB is already non negligible for the data we present (13.1 h on target). This work thus shows that the CIB is likely to be a limiting factor for future high angular resolution kSZ mapping.

By combining the two NIKA bands, under the assumption that the SZ signal is the sum of a tSZ and a kSZ contribution, we extract the first map of the kSZ effect toward a cluster of galaxies. We account for relativistic corrections using an X-ray temperature map derived from XMM-Newton data. The kSZ map is dominated by a dipolar structure with peak signal in the direction of regions B (negative) and C (positive). This is in agreement with the merger scenario in which subcluster B is a compact core falling onto the main cluster, C. However, our data indicate that subcluster B and C are themselves likely to be made of subcomponents. Our results complement that of Sayers et al. (2013) at different scales. Apart from the relativistic corrections, which assume a constant temperature along the line-of-sight, the kSZ map we present in this paper is model independent.

We use the NIKA surface brightness maps in addition to X-ray temperature estimates of each subcluster to constrain the gas density and velocity distribution in MACS J0717.5+3745. The gas density model assumes that each subcluster is described by a β-model with constant line-of-sight velocity and temperature. We perform a MCMC fit and provide constraints on the gas velocity toward each subcluster. While subcluster D is consistent with null line-of-sight velocity, subcluster B and C are moving away from us and toward us in the CMB reference frame, respectively. The relative significance of our measurement is stable and model independent, but the absolute velocity depends on our assumptions in the structure of the gas of MACS J0717.5+3745. Our constraints are limited by degeneracies between the optical depth and velocity of the subclusters, and by the simplicity of our model, which is a strong idealization for a system known to be an extremely complex cluster of galaxies.

The results presented in this paper allow us to study the details of the gas motion in the cluster and complement other observations of MACS J0717.5+3745. This shows the feasibility to study the distribution of the gas velocity within galaxy clusters, which is essential to our understanding of the dynamics of mergers at the substructure level.

Acknowledgments

We are thankful to the anonymous referee for useful comments that helped improve the quality of the paper. We are grateful to Jack Sayers for helping us obtain the Herschel SPIRE catalog we use in this paper. We would like to thank the IRAM staff for their support during the campaigns. The NIKA dilution cryostat has been designed and built at the Institut Néel. In particular, we acknowledge the crucial contribution of the Cryogenics Group, and in particular Gregory Garde, Henri Rodenas, Jean Paul Leggeri, Philippe Camus. This work has been partially funded by the Foundation Nanoscience Grenoble, the LabEx FOCUS ANR-11-LABX-0013 and the ANR under the contracts “MKIDS”, “NIKA” and ANR-15-CE31-0017. This work has benefited from the support of the European Research Council Advanced Grant ORISTARS under the European Union’s Seventh Framework Programme (Grant Agreement No. 291294). We acknowledge fundings from the ENIGMASS French LabEx (B.C. and F.R.), the CNES post-doctoral fellowship program (R.A.), the CNES doctoral fellowship program (A.R.) and the FOCUS French LabEx doctoral fellowship program (A.R.). E.P. acknowledges the support of the French Agence Nationale de la Recherche under grant ANR-11-BS56-015. This research has received funding from the European Research Council under the European Unions Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement No. 340519.

References

All Tables

All Figures

Fig. 1 Top: NIKA maps of MACS J0717.5+3745 at 150 GHz (left) and 260 GHz (right). The gray contours show the significance in units of standard deviation. They are multiples of 2σ at 150 GHz and 1σ at 260 GHz, starting at ±2σ. Both maps have been smoothed to have the same effective resolution of 22 arcsec FWHM, as represented by the white circle on the bottom left corner of the maps. The regions defined in Table 2 are shown as gray dot-dashed circles. Bottom: same as top-panel figures, in the case of half difference maps computed from two equivalent sub-samples. The steps between contours is 1σ for both frequencies.

In the text

Fig. 2 Point source model at 150 GHz (left) and 260 GHz (middle), and point source mask constructed from the 150 and 260 GHz point source models (right). The color scale of the models is the same as in Fig. 1. We also provide the labels of each source on the figure, in dashed circles for radio sources and dotted circles for submillimeter sources. The point sources subtracted maps can be seen in Fig. 7.

In the text

Fig. 3 XMM-Newton and Chandra observations of MACS J0717.5+3745. Top left panel: exposure- and vignetting-corrected Chandra photon count image in the [0.5−2.5] keV band. The image has a resolution of 3.1′′. Top right image: wavelet de-noised map of the Chandra image shown in the top left panel. Iso-contours are shown with black solid lines to enhance substructures. The four regions in Table 2 are highlighted with black circles. Bottom left image: XMM-Newton wavelet de-noised temperature map, obtained using a threshold cut of 1σ. Contours from the wavelet map image are overlaid using white contours. Bottom right image: same as bottom left panel except that we show the temperature map computed using the Chandra dataset. All the maps shown are centered on the X-ray peak.

In the text

Fig. 4 Constraints on the kSZ and the tSZ amplitudes toward the regions associated with the four subclusters (from A to D). The black data points provide the average brightness measured by NIKA in the 22 arcsec radius regions with center coordinates given in Table 2 and shown in Fig. 1. The shaded areas give the 68% confidence level constraints on the kSZ (red), the tSZ (blue) and the sum of the two (gray). Possible point source residuals were masked when computing the average surface brightness.

In the text

Fig. 5 Map of the thermal SZ effect, ytSZ (left), and the kinetic SZ effect, ykSZ (right). Gray contours are multiples of 1σ, starting at ±2σ. The map effective resolution, 22 arcsec, is shown as a white circle on the bottom left corner. Subcluster regions are represented in gray. The positions of the identified point sources (see Sect. 3), from which mis-subtraction residuals remain a potential local contaminant to the signal, are represented in light gray.

In the text

Fig. 6 Comparison between the kSZ signal and other probes of the cluster components, shown on saturated linear scales. The astrometries are the same as those in Fig. 5, and the green contours reproduce the kSZ signal-to-noise ratio, given as a solid line if positive (+2 and +3σ) and dashed line if negative (−2, −3, −4 and −5σ). Our reference regions are provided as gray circles. The positions of B1 and B2 are also shown as gray diamonds on the left image. Left: Chandra X-ray photon count image, smoothed to 7 arcsec FWHM. Middle: HST optical combined image from the CLASH data (Postman et al. 2012). Right: lensing mass (κ) from Zitrin et al. (2011).

In the text

	Fig. 7 Comparison between the data and the best-fit model. From left to right, the figure shows the input data subtracted from point sources, the best-fit model F1, the residual F1, the best-fit model F2, and the residual F2. The top row provides the 260 GHz data and the bottom row the 150 GHz data. Contours are multiples of 2σ at 150 GHz and 1σ at 260 GHz, starting at ±2σ.
In the text

Fig. 8 Constraints of the model on the ICM distribution in the plane line-of-sight velocity – central optical depth, ${\mathit{v}}_{\mathit{z}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}$ – ${\mathit{\tau }}_{\mathrm{0}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}$, for subclusters B–D. The purple, dark blue and light blue are constraints contours at 68, 95 and 99% confidence limit in the case of our baseline fit (F1, without X-ray imaging prior on the gas density). For subcluster B, the red contours are similar to the blue ones in the case of the fit F2 (with the X-ray imaging prior). The marginalized probability density distribution are also given by the histograms in the case of F1. The dashed line give the zero location of the velocity axis. As discussed in the text, given the limitation of the models to describe the cluster, we stress that the error contours do not necessarily reflect the true underlying uncertainties.

In the text

Fig. 9 Top: map of the best-fit optical depth model. Bottom: projected gas line-of-sight velocity map. The resolution of the map is 22 arcsec FWHM as provided by the white circle in the bottom left corner. The left panels correspond to fit F1 and the right panels to fit F2. Regions where the rms of the noise is larger than 3 times the minimum value have been masked for display purpose. Contours show the detection level of the velocity (thus proportional to that of the kSZ map of Fig. 5) and are spaced by 1σ, starting at ±2σ. We stress that they do not reflect the true uncertainty on the velocity because they do not account for the statistical and systematic effect of the optical depth models. The field shown in here is smaller than the other ones displayed in this paper.

In the text