ZWCL 1856.8: A Rare Double Radio Relic System Captured within NuSTAR and Chandra Field of View

Observations of galaxy cluster mergers provide insights on the particle acceleration and heating mechanisms taking place within the intracluster medium. Mergers form shocks that propagate through the plasma, which result in shock/cold fronts in the X-ray, and radio halos and/or relics in the radio regime. The connection between these tracers and the mechanisms driving non-thermal processes, such as inverse Compton, are not well understood. ZWCL 1856.8 is one of the few known double radio relic systems that originate from nearly head-on collisions observed close to the plane of the sky. For the first time, we study NuSTAR and Chandra observations of this system, which contains both relics within their field of view. The spectro-imaging analyses results of the system suggest weak shock fronts with M numbers within 2σ of the radio-derived values, and provide evidence of inverse Compton emission at both relic sites. Our findings have great uncertainties due to the shallow exposure times available. Deeper NuSTAR and Chandra data are crucial for studying the connection of the radio and X-ray emission features and for constraining the thermal versus non-thermal emission contributions in this system. We also present methods and approaches on how to investigate X-ray properties of double relic systems by taking full advantage of the complementary properties of NuSTAR and Chandra missions.


INTRODUCTION
According to large scale structure formation scenarios, clusters of galaxies are hierarchically formed by the merger of smaller scale structures.Galaxy cluster mergers are the most energetic (10 64−65 ergs) processes in the universe, which drive shocks and turbulence into the intracluster medium (ICM).As dark matter halos and galaxy velocity distributions revirialize relatively slowly, the kinetic energy of the merger is transformed more quickly within the collisional hot gas, driving shock and cold fronts that heat and mix the thermal gas while also (re)accelerating electrons and cosmic rays through Fermi-like acceleration processes (see, e.g.Brunetti & Jones 2014).
Shock fronts in the X-ray regime are evidenced by sharp surface brightness (SB) discontinuities accompa-nied by a temperature drop from the upstream to the downstream of the front.Merger shocks are characterized by a Mach number of < 3 (see, e.g.Gabici & Blasi 2003;Ryu et al. 2003).Cold fronts, on the other hand, are due to the motion of cooler subcluster in the ambient gas with higher entropy with respect to the subcluster.Cold fronts exhibit similar SB discontinuity behaviour like the shock fronts, but the temperature drop direction is reversed; the temperature ahead of the cold front is higher than the temperature behind (Markevitch & Vikhlinin 2007).
In addition to shock and cold fronts, mergers of galaxy clusters form extended radio structures in the form of radio halos and relics.Galaxy clusters hosting these structures provide information on the astrophysical particle acceleration mechanisms as well as non-thermal processes that take place within the ICM (Finner et al. 2021).Radio halos are thought to be formed by the turbulence introduced by merger shocks, and are found in the central regions of these systems.Radio relics (radio shocks), on the other hand, are anticipated to trace large-scale shock waves propagating through the ICM (de Gasperin et al. 2014).Relics usually reside in galaxy cluster periphery, since half of the virial radius from the cluster center is where the kinetic energy dissipated in merger shocks peaks (Vazza et al. 2012).Further details on the particle acceleration mechanisms forming extended radio structures and their efficiencies are provided by Bonafede et al. (2012); van Weeren et al. (2019).Double radio relic formation, which are nearly equidistant from the cluster center, are quite rare.These structures are thought to have formed as a result of a single merger event and merely 12 double relics are identified in literature (van Weeren et al. 2019).This scarcity is due to the merger geometry; i.e., the merger has to take place near the plane of the sky for both relics to be visible to the observer, in addition to the requirement of nearly head-on collisions of similar mass systems.Since each of the double relic pair trace the same single merger event, they provide the advantage of putting tighter constraints on the merger scenarios (Finner et al. 2021) and particle acceleration mechanisms.
The coincidence of X-ray detected shock fronts and radio detected relics in merger systems is evidenced by multiple studies, but require further confirmations (van Weeren et al. 2019).To confirm shock fronts at relic locations, where the S/N is low due to the large distances from the X-ray peaks, deep observations with high spatial resolution are required to locate the surface brightness discontinuities, as well as wide X-ray bandwidth to constrain the temperature drops across these edges.These confirmations are extremely important to shed light on the particle acceleration and ICM heating mechanisms.
Wide X-ray band also provide key insights into the conditions and mechanisms that are operating in the ICM with the capability of disentangle the thermal and non-thermal emission contributions.The same relativistic electrons producing the synchrotron-powered radio halos also upscatter cosmic microwave background (CMB) photons to X-ray energies via inverse Compton (IC) process.The ratio of fluxes between the synchrotron and this IC emission is simply the ratio of the energy densities of their respective radiation fields, i.e.; the magnetic field strength B and the CMB, respectively.Since the latter is well known, an IC detection or upper limit leads to an estimate or lower limit on the volume-averaged value of B, a quantity that is poorly constrained in galaxy clusters in general, which has significant implications on the assumptions used for galaxy cluster mass estimates (see, e.g.Biffi et al. 2016;Sarazin et al. 2016).IC emission in galaxy clusters requires better constraints to confirm that the dynamical role of magnetic fields in clusters are negligible, as is currently assumed in simulations and mass scaling relations used in cosmology (e.g., (Vikhlinin et al. 2009)).
To constrain the IC emission within the ICM, hard Xray data are required that enable disentangling the thermal bremsstrahlung and relatively weak non-thermal emission features in the ICM.Narrow X-ray band is incapable of disentangling the tail of IC emission that extends beyond the thermal bremsstrahlung turn over at higher energies (Wik et al. 2014).
Chandra and XMM-Newton data show a double peaked cluster structure driven by the merging subclusters (Finner et al. 2021;Jones et al. 2021).Using Monte Carlo Merging Analysis Code (MCMAC; (Dawson 2013)) Finner et al. (2021) suggest that ZWCL 1856.8 is at an early stage merger; returning from the first apocenter.A weak-lensing study by Finner et al. (2021) report M 200 =2.4×10 14 M ⊙ for the system, and defines the system as a major merger with a near 1:1 mass ratio.They assign masses of M 200 =1.2×10 14 M ⊙ and M 200 =1.0×10 14 M ⊙ to the north and south subclusters, respectively.While they report a global temperature of kT =3.7 keV and kT =4.3 keV for the northern subcluster using XMM-Newton data, they cannot constrain the temperature for the southern subcluster.
Full-polarization WSRT radio observation study at 1.4 GHz by de Gasperin et al. (2014) report the linear sizes of the relics to be ∼0.9 and ∼1.4 Mpc for the northern and the southern relics, respectively in addition to the hints of a central radio halo.They cannot report radio spectral indices due to the insufficient data.Jones et al. (2021) report the sizes as ∼0.9 and ∼1.5 Mpc using LOFAR data at 140 MHz.They report integrated (injection) radio spectral indices of α int(inj) =-0.95 (-0.87) and α int(inj) =-1.17 (-0.97) for the northern and southern relics respectively.Assuming the shock fronts lie at the outer edge of the relics, Jones et al. (2021) find spectral steepening towards the cluster center due to the downstream synchrotron and IC losses.They report that the Mach number derived from the integrated spectral indices, M int , is unconstrained for the northern relic and is M int = 3.6 ± 0.7 for the southern relic.Whereas Mach numbers derived from the injection spectral indices, M int , are M inj = 2.5 ± 0.2 and M inj = 2.3 ± 0.2 for the northern and southern relic, respectively.Guided by the low-resolution LOFAR image, they also add that the extended emission lying to the east of the cluster in the region between the relics, may be a candidate radio halo.In addition, strong polarization in the northern relic and weak polarization in the southern relic points to a slightly tilted view from the plane of the sky (de Gasperin et al. 2014).
At the cluster redshift, NuSTAR and Chandra were able to observe both of the radio relics within one field of view (FOV).This was the main motivation behind requesting NuSTAR follow-up observations for obtaining hard X-ray information on the system.To elaborate; when both the regions of interest and the bright galaxy cluster centers do not fall within the same NuS-TAR FOV, scattered light becomes difficult to tackle with, since bright sources near but out of the NuS-TAR FOV causes energy and position dependant emission contamination within the FOV (Tümer et al. 2022).Hence, this NuSTAR observation is unique in a way that provides hard X-ray double relic data without scattered light contamination.
In this paper, we investigate the X-ray emission features of this merger via spectro-imaging analyses of NuSTAR and Chandra data.The paper is organized as follows: observations, data reduction processes and the background assessment of the NuSTAR and Chandra data are presented in Section 2. In Section 3, the spectro-imaging analyses are described and their results are presented.We discuss our findings in Section 4, and conclude our work in Section 5 including the plan for future work.
Throughout this paper, we assume the ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.27, Ω Λ = 0.73.According to these assumptions, at the cluster redshift, a projected intracluster distance of 100 kpc corresponds to an angular separation of ∼22 ′′ .All uncertainties are quoted at the 68% confidence levels unless otherwise stated.Each color represents a region selected for the background fit.For plotting purposes, adjacent bins are grouped until they have a significant detection at least as large as 8σ, with maximum 12 bins.
In order to filter the data, standard pipeline processing using HEASoft (v.6.28) and NuSTARDAS (v.2.0.0.) tools were used.To clean the event files, the stage 1 and 2 of the NuSTARDAS pipeline processing script nupipeline were utilized.Regarding the cleaning of the event files for the passages through the South Atlantic Anomaly (SAA) and a "tentacle"-like region of higher activity near part of the SAA, instead of using SAAMODE=STRICT and TENTACLE=yes calls, we produced light curves and manually removed increased count rate intervals using 100 s time bins to create good time intervals (GTIs) without fully discarding the passage intervals.

Tümer et al.
This new set of GTIs were then reprocessed with nupipeline stages 1 and 2, and cleaned images were generated at different energy bands with XSELECT.We used nuexpomap to create exposure maps in the 3.0 -8.0 keV and 8.0 -15.0 keV bands accounting for vignetting.To produce the corresponding spectra for the regions of interest as well as the corresponding Response Matrix Files (RMFs) and Ancillary Response Files (ARFs), stage 3 of nuproducts pipeline were used.
We modeled the NuSTAR background using a set of IDL routines called nuskybgd that has become a standard for background assessment of extended sources in the last decade.The treatment of the background and its components are explained in detail in Wik et al. (2014).We applied nuskybgd to the spectra extracted from cluster emission free regions at each detector but still accounting for a residual ICM emission with a single temperature model.The joint fit of the background models as well as the contribution from the residual ICM emission for each FPM are shown in Figure 1.This global background model was used to create background images, which were then subtracted from the images and were corrected by the corresponding exposure maps.Background subtracted, exposure corrected images in the 3.0 -8.0 keV and 8.0 -15.0 keV bands are presented in the left and middle panels of Figure 2.

Chandra
Chandra observed ZWCL 1856.8 in 2018 with Chandra Advanced CCD Imaging Spectrometer Imaging (ACIS-I) array with Observation ID : 19752 for 42.4 ks of exposure time.For data reduction and spectral extraction we utilized Chandra Interactive Analysis of Observations (CIAO) version 4.14 package with CALDB 4.9.8.In order to create level-2 files, we used the standard CIAO script chandra repro.Then we cleaned the data from solar flares using lc clean.
wavdetect algorithm was used to detect point sources followed by a visual inspection.We used acis bkgrnd lookup script to obtain the blank-sky background file from the CALDB that matched the observation epoch.These blank-sky background event files were used to obtain the flux images.Background subtracted, exposure corrected images in various energy bands are presented in Figure 3.

Imaging Analysis
In order to search for substructures via surface brightness variations in the cluster, we applied Gaussian Gradient Magnitude (GGM) filter (Sanders et al. 2016) to the background-subtracted, exposure corrected Chan-dra photon image in the broad (0.5 -7.0 keV) band.The GGM filter calculates the gradient of an image assuming Gaussian derivatives with a width of σ (pixels).This method captures gradients at different scales depending on the selection of the width, i.e., small σ values are used at the central regions where there are many counts, and large σ values are better at capturing gradients at cluster outskirts that have lower photon counts.To create the GGM filtered image, the image itself is convolved with the gradient of a 1D Gaussian function for two axes, then these two resulting images are combined for the 2D gradient image.
In Figure 4 we present GGM filtered images with σ = 2, 4, 8, and 16 pixels, where we overlaid LOFAR contours of the radio relics.Due to the lack of deeper data, GGM σ = 1 image is mainly dominated by the noise, therefore we excluded it from this figure.This also can be said for σ = 2, however, it is somewhat important to show where gradients begin to appear.
A strong, asymmetrical gradient connecting the two subclusters appears with σ = 4.In the same image, we see that northern subcluster is more pronounced than the southern subcluster.GGM filter with σ = 8 hints at a surface brightness edge within the southern relic polygon region, where a similar edge seems to reach the inner part of northern relic.At σ = 16, GGM shows a relatively symmetrical gradient distribution, yet sharper edges are present in the northern regions.
Since the main motivations of this work are to search for any shock features coinciding with the relics, as well as any radio relic -IC connection, we defined the regions of interest as explained in the following section.

Temperature distribution within the ICM
In line with the motivation of this work, we selected seven regions of interest."Center" region, defined as an ellipse, is selected to assess the subcluster features of the system enclosing the extended emission pronounced in photon images (Fig 2 and Fig 3).We then defined two polygons tracing the northern (NR) and southern (SR) relics from the LOFAR data, to study the thermal and non-thermal emission contributions.Four sectors are selected to measure the temperature of the inner and outer regions of the northern and southern relics, i.e.NR in , NR out , SR in and SR out .These regions are shown in Figure 5.
NuSTAR spectra were extracted using nuproducts, followed by background extraction using nuskybgd that creates spectra from the model defined across the FOV as defined in Section 2. For analyzing the Chandra data, we extracted the source spectra using specextract task,     which also creates the weighted RMFs and ARFs.A local background extraction approach from the same observation was chosen given the redshift of the cluster.With this specextract task, we simultaneously extracted a background spectrum from a region with r = 3 ′ in the field of view where no cluster emission is expected.For spectral fitting, we use XSPEC1 .We apply maximum likelihood-based statistic (hereafter, C-stat) that is appropriate for Poisson data as proposed by Cash (1979).Photon counts used in spectral analysis for both NuSTAR and Chandra spectra are grouped to have at least 3 counts in each bin.We could only use 3.0 -10.0 keV band of NuSTAR's bandpass due to low S/N at high energies where background is dominant, and used 0.5 -7.0 keV band for Chandra.While the hard energy information beyond 10.0 keV is extremely important to differentiate the thermal vs. non-thermal emission components, slightly lower global temperature of this cluster with respect to a typical violent merger encourages the spectral analysis with the data at hand even in the 3.0 -10.0 keV band.

476.64/489
To begin with, we fitted the spectra of NuSTAR and Chandra individually.The initial assessment of the temperatures were achieved using a single temperature apec model whose redshift and metal abundance parameters were fixed at z = 0.304, and Z 1 = 0.3 Z ⊙ , respectively, due to low photon statistics at hand.This model was then convolved with the foreground hydrogen column density modeled by phabs, where the N H parameter was fixed at 4.45×10 20 cm −2 (based on Leiden/Argentine/Bonn Galactic HI survey (Kalberla et al. 2005)).We used the abundance table of Wilms (Wilms et al. 2000).This single temperature model assumes an isothermal ICM.
We then fit the NuSTAR and Chandra spectra jointly.During this procedure, the apec temperature, abundance and redshift parameters were linked within the instruments.We accounted for the cross-calibration differences of detectors with the constant model, where the normalization of the models of different instruments was tied.The joint-fit model hence becomes; constant × phabs × apec.For the Center region, where we have the highest S/N, we could allow the cross calibration constant parameter be free.Namely, one instrument constant is fixed to 1, whereas the other instrument constant parameters are left free to vary.The cross-calibration constants of Chandra (0.94) and NuSTAR FPMs are within 6%, and FPMA constant (fixed to 1) and FPMB (0.92) constant are within 8% of each other.The temperature resulting from this fit is 6.30 +0.37  −0.35 keV with C-stat/d.o.f.(C/ν) 434.56/514.The NuSTAR and Chandra joint spectral fit for the Center region using constant × phabs × apec is shown in Figure 6 as an example of our fitting procedure.
As we move away from the central region to the cluster outskirts, where the S/N decreases and the background begins to dominate, the cross-correlation constants become incongruent, pointing to discrepancies in NuSTAR and Chandra fluxes and/or systematics in the background subtraction, resulting in unconstrained and unrealistic temperatures.Hence, we removed the crosscorrelation constant parameter and let normalizations of all instruments be free and kept temperature, abundance and redshift parameters still tied amongst the instruments.Only then the joint-fit temperature values reflected what was obtained from the individual spectral fits.
We note that although this approach gives new information on the temperature parameter, it limits any interpretations of the densities that are obtained from the normalization values of these joint fits since all normalization are free and the which instrument represent the correct flux is not obvious.We present the fit parameter values of individual fits as well as the joint fits with free instrument normalizations in Table 1.The temperature results are plotted in Figure 7 with regions ordered from North to South to provide a visualization of a cluster radial profile.
We also fitted the Center region with two other models; a two temperature thermal ICM (2T; constant × phabs × [apec + apec]) and a thermal plus non-thermal model (1T+IC; constant × phabs × [apec + powerlaw]).We let the cross correlation constants free for these procedures.
The resulting temperature values from 2T model temperatures are 6.75 keV and 1.45 keV, both unconstrained.The statistics yield C/ν = 434.01/510,and constant parameters are 1.00 (fixed), 0.91 and 0.94 for NuSTAR FPMA, NuSTAR FPMB, and Chandra, respectively.Only when these constant parameters are fixed to their best-fit values, we were able to find a lower limit to the high temperature component of 6.21 keV.
The 1T+IC model results in a apec temperature value of 6.45 +0.81  −1.57keV and a powerlaw photon index (Γ) value of 2.41 (unconstrained) with C/ν = 434.43/510and same constant values found for 2T model.When we fixed the constant parameter, Γ was still unconstrained and the temperature and the corresponding errors become 6.45 +0.67  −1.47 keV.Finally, we removed the constant component from the 2T and 1T+IC models, leaving all instrument normalizations free and we present the results in Table 2.
To quantify the systematics of the background characterization, we applied a generous ±15% change to the total background of the central region spectrum of Chan-dra and NuSTAR spectral fits, and reran the individual fits.The temperature value changes were within 1.1 σ for Chandra and 0.9 σ of the statistical errors on the temperature of that region.Given that the background level is unlikely to be off by this amount, we conclude that the measurements are not particularly sensitive to it.

Relics
For the relic regions, we also investigated the possible non-thermal emission components.This was achieved in the similar manner that was applied to the Center region, leaving the cross normalizations of instruments untied and no additional constant parameter.We applied a two temperature thermal ICM (2T), and a thermal plus non-thermal model (1T+IC) to the joint NuS-TAR and Chandra spectra.2T model does not only test the possibility of multi-temperature structure, but also is a good judgement of the validity of the 1T+IC model.Mainly, if the statistics favour 1T+IC and 2T models equally well, then the conclusion of a dominant IC emission within the ICM requires further investigations.When the low and high temperature components of the 2T model have close values, this is a good confirmation that the ICM can be described well by a single temperature model.On the other hand, 2T model also emphasizes idiosyncrasies of the data (low S/N, background assessment issues) when it results in unconstrained or unphysical temperature values.Therefore, it is a good practice to apply 2T model to the data prior to making conclusions about the contribution of non-thermal emission.
Although thermal emission should be present at the relic sites we also fit each relic region spectra jointly with NuSTAR and Chandra spectra using a single IC model.This approach provides a higher limit to the detection significance with the assumption that ICM emission is entirely dominated by the non-thermal IC emission.
Northern Relic: 2T model fit to the joint NuS-TAR and Chandra NR region spectra results in a lower temperature component of kT = 0.40 +0.29  −0.27 keV and a higher temperature component of 3.25 keV, with a lower limit of 1.55 keV.The resulting statistics is C/ν = 86.72/98.
When we reintroduced the constant parameter with values fixed to Center region values (hence tied the normalizations), one of the temperatures hit the 64 keV hard limit, and the other temperature became kT = 1.11 +0.54  −0.31 keV with C/ν = 88.28/100.For comparison, 1T model fit to the joint NuSTAR and Chandra spectra of the same region is C/ν = 89.54/102(constant parameter.)1T+IC implementation of this region results in unconstrained temperature and photon index values.When we repeat this fit by fixing the photon index to Γ = 1.87 obtained from α inj reported by Jones et al. (2021) (Γ = α+1), we get kT = 1.02 +0.39  −0.61 keV with C/ν = 87.27/99.This is based on the assumption that the synchrotron emitting electron population is also responsible for the IC emission.
Single powerlaw fit to the NR results in a Γ = 4.31 +0.91  −0.79 with C/ν = 86.14/102.This fit implies a significance of σ = 2.1.We also repeated this fit by fixing the Γ to 1.87 implied by the radio spectral index found by Jones et al. (2021).With this approach, the detection significance becomes; σ = 1.9 with fit statistics C/ν = 96.21/103 Southern Relic: When applied these same procedure steps to the southern relic.2T model results in a higher and lower temperature components of kT = 9.84 keV(unconstrained) and kT =3.40 +7.23  −2.45 keV, respectively, both unconstrained.The resulting statistics for this fit is C/ν = 245.79/261,where 1T model fit to the joint NuSTAR and Chandra spectra of the same region is C/ν = 245.78/263.Again, we introduced the constant parameter values fixed to the Center values, and the temperatures for both models become kT = 8 keV (unconstrained) with C/ν = 250.46/263.
1T+IC modeling of this region results in an extremely volatile photon and temperature values.We repeated the same procedure applied to the northern relic 1T+IC fit by fixing the Γ to 1.97 using α inj value Jones et al. (2021), we find kT = 2.30 +1.05 −1.50 keV with C/ν = 245.15/260.
When we applied the single powerlaw model, we find a photon index value of Γ = 2.49 +0.51  −0.44 with C/ν = 247.06/263.With the aforementioned assumption of the origin of the emission and model use, we put of detection limit of at least 5σ.When we use Γ = 1.97, the significance becomes 6σ.
In Table 2, we provide a comparison table of the findings from 1T, 2T, and 1T+IC model implementations on the Center, NR and SR region spectra and additional IC model results to the NR and SR regions.

Cross-talk Analysis
Even single bright sources in the NuSTAR FOV cause a scatter of photons due to the large (∼1 ′ Half Power Diameter, ∼18 ′′ Full Width at Half Maximum), slightly energy-dependent point spread function (PSF).This results in a cross-contamination, namely crosstalk, of multiple emission features in the regions of in-regions Center, NR and SR, with the additional IC model results of NR and SR spectral fits.Abundances are fixed to Z1 = 0.3 Z⊙ and instrument normalizations are left free to vary.terest, although the emission may not purely originate from the selected region.Standard NuSTAR pipeline, nuproducts, produces ARFs for point or diffuse sources inside the user defined extraction regions.However, it does not account for the ARFs for other sources whose emission originates outside these extraction regions that contaminate the spectra of those regions, which will be referred to as cross-ARFs.
A set of IDL routines created to account for this contamination that divorces the contamination and source emission is nucrossarf 2 .This method is explained, applied and tested in detail by Tümer et al. (2023).
Given the lack of good photon statistics, we were forced to simply our regions of interest for the cross-ARF analysis to increase the photon statistic.We kept the Center, NR, and SR regions to assess the central temperature more accurately, as well as to constrain the relic emission better to study the implications of a possible IC emission.We also selected a circular region encircling but excluding the aforementioned regions, which we will refer to as; "Circle".This way, the contribution from both the dominant Center emission and the cooler gas beyond the central region is separated.
While running the nucrossarf an order is assigned to each region.In our treatment, the order from 1-4 is designated to Center, NR, SR, and the Circle regions, respectively.To assess the cross-talk, a model is assigned to each source distribution.Whereas we expect a dominant thermal ICM emission (apec) from the Center and Circle regions, we applied a powerlaw model to represent a non-thermal (IC) emission coming from the relics to assess the statistical significance assuming the relic emission is purely from the IC.Although, nucrossarf allows assigning multiple models to each region emission for a more accurate representation of complex structures, the lack of good photon statistics of the current data dictates a simplistic approach.Namely; we expect a considerable amount of photons to leak inside the relic regions from the more dominant thermal ICM of the Center region, hence this will decrease the number of photons originating from the relic regions.Since the initial fits of the relic regions already resulted in large uncertainties, especially in the SR region (Table 1), this will worsen with fewer photons once the cross-talk is applied.
apec model redshift and metal abundance values were fixed to z = 0.304, and Z 1 = 0.3 Z ⊙ , respectively as in the original fits without the cross-talk analysis.
The visual 4 × 4 matrix given in Figure 8 presents the mapping of the contribution from the local photons within a region and the cross-talk, where the selected regions are overlaid.In this M × N matrix, M represents the local source PSF and N represents the region number, i.e., the box at 2 × 4 shows the PSF leakage from the source situated at region 2 into Region 4. In other words, i × i boxes represent the PSF distribution confined with the source region, whereas i × j shows the extend of that PSF with other regions of interest.
The spectral distribution of the source emission, the individual source ARFs and the cross-ARFs are shown in in Figure 9.This plot is the joint fit result of all 4 source distributions with their respective models, as well as the model contributions from other regions.
The resulting fit gives ICM temperature of kT = 6.27 keV with a lower limit of kT = 5.98 keV for the Center region, and kT = 1.99 +0.73  −0.28 for the Circle region.
With the data at hand, we find that they are consistent with IC-like emission at the 2.3σ and 1.1σ significance levels for the northern and southern relics, respectively, where Γ is fixed at 2 as suggested by the radio spectral index found by Jones et al. (2021).When we free the Γ parameter, the fit reaches to the lower limit of Γ=1.More data are needed to conclusively determine the presence of IC emission in these regions.

DISCUSSION
Full band Chandra photon image and surface brightness contours clearly show X-ray peaks located near the BCGs of the subclusters and bright emission in a region connecting these halos.Chandra images do not indicate the inverse-S shape seen by XMM-Newton images, which is interpreted as the merger having a non-zero impact parameter by Finner et al. (2021).However, GGM filter of the Chandra image (Figure 4, upper right panel) clearly shows an strong arc-like gradient to the east, supporting this claim.Moreover, another strong gradient to the north (Figure 4, upper right and lower left panel) may indicate that the northern halo has stripped material from the southern halo during the first core passage.NuSTAR image in the soft band (Figure 2, left panel) do not reveal any important information, yet the lack of strong emission in the hard band (Figure 2, middle panel) points to the lack of hot ICM or strong nonthermal emission, yet deeper data is required to make a definite claim.
The spectra extracted from the central region, Center, is best described by a single temperature model indicating a lack of a cool core remnant as well as a robust non-thermal emission.
In addition, to compare our results with the literature (Finner et al. 2021), we also used a central r = 3 ′ region for spectral fitting with Chandra data to compare XMM-Newton and Chandra results.Our temperature results of Chandra fit agree with XMM-Newton results within 1σ for this same region selection.While Chandra gives a temperature of kT = 4.00 +0.43  −0.36 keV using N H = 9.90×10 20 cm −2 , the column density value used by Finner et al. (2021), XMM-Newton temperature is kT = 3.7 +0.6 −0.5 keV.We applied four different models to the NR and SR region spectra to investigate single (1T) and two temperature (2T), thermal plus non-thermal (1T+IC) and pure non-thermal emission by jointly fitting the NuSTAR and Chandra data.Application of 2T model was aimed at testing the statistical significance of a thermal plus nonthermal model 1T+IC) over a two temperature model rather than searching for a physical multi-temperature plasma structure.
1T+IC modeling of both relic regions results in unconstrained temperature and photon index values as well, which can be attributed to the lack of S/N that prohibits the fit to find the best-fit values since 2T model implementation is also unstable and problematic.Assuming normalizations hence higher detection significance of any model that would be applied.
It is also important to note that NuSTAR and Chandra temperature results are in a very good agreement (Figure 7) where we have moderate S/N; namely NR, NR in , Center, SR in , and SR regions for the ∼2 -7 keV range.

CONCLUSION AND FUTURE WORK
In this work, we analyzed 30 ks raw exposure of NuS-TAR and 43 ks raw exposure of Chandra data of a rare double radio relic system ZWCL 1856.8.We found an indication of a shock front coinciding with the northern radio shock (relic) with both NuSTAR and Chandra spectra, and evidence for weak IC emission at the relic sites using NuSTAR data.
Combining the distinct superiorities of these two instruments provides the assessment of a more complete description of the X-ray emission through wider combined bandpass, better capability to constrain multiple temperatures and non-thermal contribution to the signal.
We have established a ground work for the NuS-TAR follow up observation of this system that is scheduled for 240 ks during Cycle-9.In addition, Chandra will also observe the system with a deep exposure (∼200 ks) for Cycle-24.These two sets of data together will provide strong claims on the existence and locations of the shock/cold fronts as well as a robust detection or refutation of IC emission at the relic sites.This will be possible not only by increasing the photon statistics of both satellite data on the cluster to constrain the model parameters better, but with higher photon statistics, we will be able to use NuSTAR's bandpass well beyond 10.0 keV (to at least 15.0 keV) that will enhance the confidence of the possible detection limits of the non-thermal components in the ICM.Furthermore, by obtaining better constraints on the parameters of the 1T+IC models, reliable upper limits on the IC, hence lower limits on the magnetic field strength at relic sites will be driven.

Figure 1 .
Figure1.Joint-fit of background and cluster emission of NuSTAR FPMA (upper panel ) and FPMB (lower panel ).Each color represents a region selected for the background fit.For plotting purposes, adjacent bins are grouped until they have a significant detection at least as large as 8σ, with maximum 12 bins. e+00

Figure 3 .
Figure 3. Background-subtracted, exposure-corrected, Chandra flux images of ZWCL 1856.8 in broad, soft, medium and hard bands smoothed with a Gaussian of r=5 pixels (1 pixel subtending 0.492 ′′ ).Zoomed in images focusing on the BCGs are presented on the lower panels with Chandra X-ray contours.

Figure 4 .Figure 5 .
Figure 4. Chandra GGM filtered background subtracted, exposure corrected images in the 0.5 -7.0 keV band with σ = 2, 4, 8, and 16 pixels from; upper left, upper right, lower left and lower right panels, respectively.White polygons roughly trace the radio relics seen by LOFAR.The scales emphasize the regions with largest SB gradients and does not accurately reflect the magnitude of the jumps (Sanders et al. 2016).Large gradients such as edges, are indicated with lighter colors and flatter SB areas are darker colors.

Figure 6 .
Figure 6.Central region spectra showing NuSTAR (black and red) and Chandra (blue) joint fits using a single temperature model.The corresponding blue, black, and red circles below the source spectra indicate the background level for each instrument.

Table 1 .
Parameter values obtained from individual and joint fits to the NuSTAR and Chandra spectra using single temperature model from the regions show in Figure 5 and plotted in Figure 7. NH is kept fixed at 4.45 × 10 20

Figure 7 .
Figure7.Best-fit temperatures in each region.Chandra temperature is unconstrained at NRout hence is not shown.NuSTAR temperature upper limit of ∼42 keV (Table 1) in addition to NuSTAR and Chandra joint fit temperature upper limit of ∼58 keV are cropped from the top to enable better visualization of other error bars.

2Figure 8 .
Figure 8. nucrossarf analysis images showing the relative emission scattered by the PSF from one region into the others.Row numbers on the left indicate the region from which source emission originates, and the column numbers on the top indicate the region number into which the source photons leak.The diagonal images show the actual source image for each of the four regions.The scale of the logarithmic color bar has arbitrary units.

Figure 9 .
Figure 9. Joint cross-talk spectral fit from all regions shown in Fig 8, comprising of four source spectra, four source models and the 3×4 cross-talk models.

Table 2 .
This table represents the results of 1T, 2T, and 1T+IC model fits of the joint NuSTAR and Chandra