The NuSTAR View of Perseus: The Intracluster Medium and a Peculiar Hard Excess

As the brightest galaxy cluster in the X-ray sky, Perseus is an excellent target for studying the intracluster medium (ICM), but until recently, its active galactic nucleus (AGN) made studies of the diffuse emission near its center nearly impossible to accomplish with NuSTAR, due to the extended wings of NuSTAR's point-spread function. The development of a new open-source software package—nucrossarf—now allows the contribution from point and diffuse sources to be modeled so that scattered light from the AGN can be accounted for. Using this technique, we present an analysis of diffuse hard X-ray (3–25 keV) emission from the ICM using three archival NuSTAR observations of the Perseus cluster. We find a ∼10% excess of emission beyond 20 keV not describable by purely thermal models. By performing similar analyses of AGNs in archival observations, we have characterized the systematic uncertainty of the modeled AGN contribution to be 3.4%. However, in order to explain the excess, the total scattered AGN emission would have to be 39% stronger than we have measured. We test physical explanations for the excess, such as diffuse inverse-Compton emission potentially originating from the radio mini-halo, but we determine that none of the models are compelling. An upper limit on the inverse-Compton flux (≤1.5 × 10−11 erg s−1 cm−2) and a corresponding lower limit on the global magnetic field strength (≥0.35 μG) are derived. We discuss the potential origin and implications of the excess and present our characterization of the nucrossarf systematic uncertainty, which should be useful for future work.


INTRODUCTION
Galaxy clusters are the largest virialized structures in the universe, and their galaxies are immersed in a hot, ionized plasma-the intracluster medium (ICM)-that emits at X-ray energies via thermal bremsstrahlung.The brightest galaxy cluster in the X-ray sky is Perseus (Figure 1) at z = 0.018 (Struble & Rood 1999), and this proximity has allowed the investigation of the ICM and its structure in great depth (e.g., Sanders & Fabian 2007;Fabian et al. 2006Fabian et al. , 2011Fabian et al. , 2015)).One of the most important findings of the previous studies is the cluster's prominent cool core and bubbles blown by past Active Galactic Nuclei (AGN) activity.Cool cores are typically formed in dynamically relaxed clusters (i.e., clusters that have not recently undergone a major disruption to their ICM, such as a merger event) as the dense, cen-tral plasma cools radiatively more quickly than the rest of the gas in the halo.
Radio observations reveal diffuse and faint synchrotron emission associated with the cool core-region of some relaxed galaxy clusters (e.g., Gendron-Marsolais et al. 2017;Timmerman et al. 2021) known as minihalos.Mini-halos are treated separately from the largescale, merger induced extended radio sources-radio halos and radio relics-that are not directly associated with cluster galaxies (e.g., Gitti et al. 2002).Unlike radio halos and relics, these mini-halos are confined to cluster cores, often bounded by cold fronts, in clusters without other obvious merger signatures.
It is thought that the relativistic electrons forming radio halos are accelerated by the turbulence that develops in the ICM during major merger events (e.g., Tribble 1993;Brunetti et al. 2001).Similarly, the particles in radio relics (Ballarati et al. 1981) are thought to get ac-celerated by shock waves (e.g., Enßlin & Brüggen 2002).Unlike radio halos, these relics are elongated, irregularly shaped structures that appear in the periphery of the ICM.
Barring re-heating processes, radio relics and giant radio halos should not be observable after ∼ 10 8 yr (Gitti et al. 2002).However, some relaxed clusters with cool cores host a mini-halo that extends up to ∼500 kpc beyond the radio galaxy at the cluster's center.These objects are distinct from giant radio halos; they only appear in relaxed clusters with cool cores.The first of these objects was discovered in the Perseus cluster (Burns et al. 1992), and Perseus has since been a prime specimen for studying the physics of mini-halos.
The same pool of relativistic electrons responsible for synchrotron emission in radio structures is expected to uplift cosmic microwave background (CMB) photons to X-ray energies via inverse Compton (IC) scattering.The strength of IC emission depends on the number density of the relativistic electron populations.In the case of radio mini-halos, the synchrotron emission is notoriously faint, which suggests that the electron density is low, especially given the expectation that the magnetic field is strongest in cluster centers.
The synchrotron flux from a population of relativistic electrons is related to the magnetic field strength of the environment, but attempts to estimate the magnetic field using only synchrotron emission rely on the assumptions of equipartition.However, if a detectable IC signal is present, one can derive the average magnetic field strength from the ratio of the IC and synchrotron luminosity without evoking these assumptions: where u B is the magnetic field energy density and u phot is the energy density of the CMB, which is the seed photon field for IC scattering.Unlike equipartition or Faraday rotation methods [see Govoni & Feretti (2004) for an overview of these topics], using IC to measure magnetic field strength is unambiguous and does not rely on assumptions (as in equipartition) or an incomplete knowledge of gas structure (as in Faraday rotation).However, previous searches for IC-which have largely focused on merging and disturbed clusters with bright radio structures-have been inconclusive, so only lower limits to magnetic field strength have been found so far (e.g., Nevalainen et al. 2004;Wik et al. 2009Wik et al. , 2014;;Rojas Bolivar et al. 2021, 2023;Mirakhor et al. 2022;Tümer et al. 2023).
Detecting IC is tricky because the thermal ICM mostly dominates the X-ray spectra until ≳15-20 keV.
NuSTAR is the only existing imaging X-ray telescope capable of detecting a signal at these hard energies; the sensitive region for diffuse sources is ∼3-30 keV.Unfortunately, the instrumental background at ≳20 keV has rendered many past IC searches inconclusive.
Since mini-halos are much fainter than large-scale halos, the relativistic electron density is expected to be much lower than that of large-scale halos.Therefore, mini halos are not of great interest for IC searches in ICM.However, the Perseus cluster 's mini-halo (12.64 Jy at 230-470 MHz;Gendron-Marsolais et al. 2017) is brighter by ∼ 3 − 4 orders of magnitude when compared to average mini-halos (Giacintucci et al. 2014).Thus, Perseus is the best candidate for searching for IC emission from a mini-halo.
This study reports on NuSTAR observations of the core of the Perseus Galaxy Cluster, which is at redshift z = 0.018 (Struble & Rood 1999) and has an X-ray luminosity of 2.19×10 45 ergs s −1 (2-10 keV, Mushotzky et al. 1978), making it the brightest X-ray cluster in the sky.The virial mass and radius derived from X-ray measurements are M ⊙ = 6.65 × 10 14 M ⊙ and r 200 = 1.79 Mpc, respectively (Simionescu et al. 2011).Perseus is a cool core cluster with a central temperature of ∼3 keV that rises to ∼5.5 keV at 120 kpc (Schmidt et al. 2002) with a central X-ray emitting AGN and a radio mini-halo.
Throughout this paper, we assume a ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 , Ω M = 0.27, and Ω Λ = 0.73.Under these assumptions, the angular scale at the redshift of the Perseus Cluster is ∼0.37 kpc/ ′′ .All uncertainties are quoted at the 3σ confidence levels unless otherwise stated.
In this study, we report on our search for an IC signal in the ICM of the Perseus Galaxy Cluster.In Section 2, we describe observations and data reduction processes.In Section 3, we discuss the spectral analysis, our treatment of systematic uncertainties, and the characteristics of the hard excess.In Sections 4 and 5, we discuss our results and summarize our findings.Lastly, derivation of the systematic uncertainty associated with scattered emission from the AGN is described in an Appendix.

OBSERVATIONS AND DATA REDUCTION
In this work, we use data from three archival NuSTAR observations of the Perseus cluster.Two of these observations (ObsIDs 90202046002 and 90202046004) are centered on 3C 84, the radio counterpart to the bright central galaxy of the cluster (NGC 1275).The third observation (ObsID 60061361002) is offset to the north by 2. ′ 101.Together, these three observations yield 70.4 ks of data covering the central r ≲ 6 ′ of the cluster.They  are summarized in Table 1 and their combined data are shown in Figure 1.These data were processed using NuSTAR's standard pipelines (HEASoft v6.30.1 and NUSTARDAS 1.9.7).From archival-provided cleaned event files, we manually screened data for periods of excess counts caused by high solar activity or high particle background, as can be caused by passages through the South Atlantic Anomaly.Light curves in 1.6-20 keV (solar) and 50-160 keV (particle) bands with 100 s bins were constructed and periods with excess count rates were manually identified and removed from the Good Time Interval (GTI) files.This process was repeated within the solar and particle bands for light-curves within the central r ≲ 6 ′ of the cluster.NuSTAR's stable background allows this simple procedure to quickly identify periods of high rates, which appear as clear outliers relative to typical variations during an observation.This method typically removes less exposure time than conservative filtering with the SAAMODE=STRICT and TENTACLE=yes keywords set, which removes GTIs based on NuSTAR's geographic location and only catches periods of high particle background.Additionally, this manual process improves the quality of GTI filtering; high background periods occasionally fall outside of the defined geographic regions excluded by the conservative filtering method, and in the opposite case, the background may be at nominal levels even while the telescope passes through those geographic regions.The new GTIs were used to reprocess the data and create the cleaned event files with nupipeline, which have been used to create all further data products.The resulting exposure times are shown in Table 1.
Using nuproducts, we extracted higher level products from the cleaned event files.Vignetted exposure maps were generated at 10 keV using nuexpomap in order to produce exposure-corrected 4-25 keV band images (e.g., Fig. 1).
To minimize contamination from the central AGN and isolate the central cool core from the hotter gas at larger radii, higher level data products-including Auxillary Response files (ARFs), spectra, and light curves-were extracted from a circular region (r < 1 ′ ) and a concentric annulus spanning the cluster gas (1 ′ < r < 7 ′ ; see Figure 1).The spectra from both focal plane modules (FPMA and FPMB) of all three observations were jointly analyzed for this study.

Background treatment
NuSTAR's background spectrum can be separated into four distinct components that are each constrained within a reconstructed uncertainty.These components are described in detail by Wik et al. (2014).Here, we briefly describe the components and list their 1σ uncertainties.
The internal background (σ sys = 3%) originates from two sources: fluorescence (22-32 keV) and activation (all energies, most prominent above 30 keV) lines within the telescope, and a flat component that represents instrumental noise and unresolved lines.This component dominates at energies E > 20 keV.
The aperture stray light (σ sys = 8%) component is a byproduct of NuSTAR's open mast that separates the optics from the focal plane.Light from the cosmic  The normalization of this model should be within ∼10% of the expected CXB level.Outside this range, the model is unphysical and is being confused with other background, foreground, or source components.When the background for observations 60061361002 and 90202046004 were being fit, the free aperture CXB normalization dropped to ∼60% of the expected flux, so we instead fixed it to the nominal value.Because the cluster emission is brighter than the aperture CXB flux, this choice has little impact on ICM measurements.
The focused cosmic X-ray background (σ sys ∼ 40%) is the cumulative light of unresolved X-ray sources-primarily AGN-focused within the FOV.This component is relatively uniform across the FOV, although the number of sources in typical regions drives the large systematic uncertainty.The variance for this component is derived empirically from the Extended Chandra Deep Field South (ECDFS), and then adjusted conservatively for the smaller solid angle and higher flux limit sampled by galaxy cluster observations.Lastly, scattered light from the Sun (and potentially from other sources) forms the so-called solar (σ sys = 10%) background component.Wik et al. (2014) obtains the uncertainty by fitting the background of the ECDFS survey, and they find that this component is accurately reconstructed with a standard deviation of 10% after statistical fluctuations are accounted for.This same exercise is used to find the uncertainty of the instrumental background discussed above.
In order to fit these background components to our data, spectra were extracted from the regions shown in Figure 2.Even though we used exclusion regions to block most of the cluster gas, the extracted spectra are still contaminated by significant thermal emission from the ICM.In order to account for this, we include an apec component (Smith et al. 2001) in fits to the spectra.As a demonstration, the fitted background for observation 90202046004 is shown in Figure 3.
At energies ≲ 15 keV, the cluster emission far dominates the background across the entire field of view, so the background model is not well constrained and residuals are left behind.These residuals are likely due to imperfect modeling of background and cluster lines, which may suffer from small gain shifts (see Section 3.2).These residuals are not particularly concerning: Because the thermal emission far dominates in this energy band, a slightly imperfect background at ≲ 15 keV will have negligible effects on model fitting.
Near ∼20 keV, the background becomes the dominant source of emission and it can be constrained more accurately with the background model.Just before this background-dominated regime-near ∼15 keV-there is a slight excess that cannot be reasonably fit by any components in the current background model.This excess will be explored in greater detail in Section 3.
Figure 1 shows the background-subtracted, exposurecorrected image combining all three observations from both focal plane modules.The AGN-which acts as a point source at the center of the cluster-has been subtracted using PSF models from nucrossarf, which will be discussed in more detail in Section 3. Aside from the AGN, no point sources were apparent in soft (3-20 keV) or hard-band (20-50 keV) images of the cluster.There was also no evidence of stray light leaking into the FOV from bright sources 1-3 • off-axis, and the NuSTAR stray-light evaluation tool reported no issues.

Spectral analysis and results
nucrossarf1 is a software suite that generates auxiliary response files (ARFs) to model cross-talk between different sources of emission in all regions of a NuSTAR observation.This accounts for NuSTAR's large PSF wings, which cause contamination between sources (see Appendix A for more details).Tümer et al. (2023) explains this set of IDL routines in detail, and we followed their procedures to extract spectra of the AGN emission in both an inner circular region and a global annulus encompassing the cluster gas (see Figure 1).The AGN is modeled as a point source, while diffuse emission from these two regions (a circular region with r < 1 ′ and a concentric annulus with 1 ′ < r < 7 ′ ) is modeled as extended sources following the spatial distribution of emission from a NuSTAR image.The data were grouped using grppha to a minimum of three counts per bin.Using PyXspec v2.1.1,we fit the spectra from 4-25 keV using Cash statistics (C-stat; Cash 1979).We did not include the 3-4 keV band-to avoid noise from the solar background-or data above 25 keV where the source falls below the systematic limit of the modeled background.An example fit is shown in Figure 4, where the data are additionally binned for display purposes only (S/N = 10, N bins ≤ 50).The contribution from the AGN is modeled with a simple power law spectrum, and the thermal emission from both extraction regions is described by single-temperature apec models.Some residuals appear near the iron line complex at 6 keV, and above ∼ 15 keV (yellow box), a subtle excess is apparent.As discussed in Section 3.4, we find that this excess is statistically significant.

Spectral modeling
This analysis follows a similar process to past searches for IC (e.g., Wik et al. 2014Wik et al. , 2009;;Rojas Bolivar et al. 2021, 2023;Tümer et al. 2023).To search for nonthermal contribution to the spectra, we first characterize the dominant thermal emission by fitting the data with XSPEC's apec model, which models singletemperature (1T) emission.Then, to asses the possible contribution from a non-thermal component describing IC emission, we introduce a powerlaw component (1T+IC).To ensure the 1T+IC model isn't just representing multi-temperature emission, we also consider a two-temperature model (2T), which can adequately describe multi-temperature ICM emission and is well modeled by two apec components.We adopt these three models (1T, 1T+IC, 2T) for the following analysis.
All of these models are convolved with by Galactic absorption (wabs), a gain response (gain), and normalization constants (constant).The Galactic absorption was modeled by a wabs2 component that was fixed to n H = 1.4 × 10 21 cm −2 , matching Rani et al. (2018) (calculated using the nH FTOOL from HEASoft).Given that we only include events with energies E > 4 keV, the exact amount of the absorption and how it is modeled is negligible and could have been completely ignored.The gain response accounts for the slight offset between actual photon energy and channel energy, which is necessary to achieve a consistent redshift with past measurements; the redshift is left fixed to this value (see Section 1).We fixed the gain slope to 1.0, while the offset was left as a free parameter as in Rojas Bolivar et al. (2023).The normalization constant allows the model to accommodate for any flux or instrumental differences between the three observations and instruments (FPMA and FPMB).Unless explicitly stated otherwise, Galactic absorption, gain response, and normalization constants were included for the entirety of our analysis.

Modeling the AGN
In order to model the AGN contribution to the extracted spectra, we add a point source model at the AGN's location in nucrossarf.Assuming that the AGN emission is point-like, nucrossarf generates ARFs that model the expected counts from the AGN in our two regions, accounting for the energydependent scatter of the PSF.The ARF for the smaller central region is essentially identical to what would be produced by nuproducts, and the "cross-ARF" for the larger annulus accounts for the emission scattered into that region by the PSF wings.These ARFs are used to convolve a powerlaw spectrum to model the AGN counts in each region.Given perfectly-realized response functions and a non-variable AGN flux over time, the power law index and normalization would then be applied across both telescopes and the three observations; to allow for time variability and imperfect energyindependent calibration, we allow the cross-calibration constant to float between telescopes and observations (but not between the regions of a single dataset).Linking the photon index across datasets, we find a bestfit value of Γ = 1.82 +0.13 −0.01 , which is consistent with Rani et al. (2018), who finds Γ = 1.85 +0.12 −0.14 (observation 90202046002) and Γ = 1.90 +0.13  −0.11 (observation 90202046004).However, the C-statistic space was not entirely well-behaved, especially on the lower end of the parameter range.In an effort to minimize degeneracy with the thermal cluster component, we tried changing the region size from r < 1 ′ to r < 30 ′′ , but this did not improve the behavior, demonstrating the difficulty of simultaneously fitting the AGN and cluster emission without other constraints or assumptions.We only aim to accurately model the AGN emission so that it does not bias ICM measurements, so in light of this issue, we fixed the powerlaw index to Γ = 1.82 for the remainder of the analysis.This value is consistent with Rani et al. (2018), whose focus was on the AGN, thus supporting our choice.Within the outer region, we find that scattered light from the AGN constitutes ∼ 26% of the 20-25 keV emission and ∼ 4% of the 4-25 keV emission in our ICM spectra.

Modeling the ICM
Our goal is to describe the nature of the emission in the outer (1 ′ < r < 7 ′ ) region that dominates the data and is less contaminated by the central AGN.For the central region, we model the gas with an independent single temperature thermal component and, when applicable, non-thermal components with independent normalizations but with indices tied to that in the outer region.This procedure allows for extensive flexibility in modeling the inner region spectrum.While it reduces the constraining power of the inner region models, this procedure also removes as much as possible constraints on model components that describe the outer region spectrum.This conservative approach allows the data to constrain our models as realistically and accurately as possible.In the following, we focus on the model parameters found for the outer region, treating the inner region values more like nuisance parameters.

Single-Temperature Model
Up to ∼ 15 keV, the global spectrum of Perseus's ICM is adequately described by a single temperature (1T) model, with C-stat/Degrees of Freedom (hereafter C-stat/ν) = 6537.35/5983.However, this model assumes that the ICM has a uniform temperature, which is only a first approach to the system.Following Potter et al. (2023), we obtain NuSTAR and Chandra temperature maps (Figure 5) by fitting single-apec models to coarse spectra that are created by extracting background-subtracted, exposure-corrected counts from three narrow energy bands (3−5, 6−10, and 10−20 keV) in regions (r ≥ 5 pixels) spaced every fifth pixel across the FOV.The center of each region is assigned its corresponding best-fit temperature, and the resulting grid of temperatures is interpolated to produce an image.The maps in Figure 5 hint to a considerable amount of temperature variation, and a 1T model, only able to estimate an average temperature in the region, may be too simplistic.
Figure 4 shows the 1T fit to the data, and residuals appear around the iron line.We have seen such features in other NuSTAR observations of relaxed clusters; it may be that line abundances at this location in the spectrum are not well captured by our assumptions-single temperature, solar abundance for all elements-or it may indicate some very small calibration issue with the response around the Fe complex.The residuals cannot be attributed to problems with redshift or abundance, and given how weak and localized they are in the spectrum, we do not expect them to affect our broadband analysis.
Above 15 keV, a slight (∼ 10%) excess becomes apparent (see Figure 4 and 6).We explore this hard excess in more depth by testing several different models.

Two-Temperature Model
As shown in Figure 5, the Perseus Cluster has a multitemperature structure that should not be fully consistent with a 1T model.As explained in Rojas Bolivar et al. ( 2021), the thermal continua described by APEC models have a lack of strong features at high energies that ultimately allows multi-temperature structures to reduce to a two-temperature (2T) scenario.
While a 2T model describes the data with a lower fit statistic than the 1T model (C-stat/ν = 6482.52/5981),the best fit parameters for the second APEC component are unconstrained, and its temperature reaches the default upper limit of 64 keV.At these high temperatures, the turn-over of the bremsstrahlung is not visible within the analyzed energy band, so it resembles a powerlaw describing a non-thermal emission mechanism.Additionally, 64 keV is not a physically realistic temperature for the ICM, especially in a relaxed, cool-core cluster like Perseus.Hence, the hard excess does not describe a physical thermal component.
When a physically-motivated upper limit of kT ≤ 9 is applied to match the highest temperature in the global region (seen in the NuSTAR temperature map in Figure 5), the C-stat/ν becomes 6497.29/5981.

Non-Thermal Model
In order to assess the contribution of possible IC emission, we included a powerlaw component to account for possible non-thermal processes (1T + IC).The powerlaw (C-stat/ν = 6475.78/5981)was strongly preferred over a 1T (C-stat/ν=6537.35/5983)or 2T (C-stat/ν = 6482.52/5981)model.However, its best-fit results are inconsistent with the expectations from a diffuse IC component.When left free to vary, the best-fit photon index (Γ) dropped to 0.254.This indicates that the excess emission -if described by a powerlaw -is extremely flat.Gendron-Marsolais et al. (2021) found that the radio index of the Perseus Cluster to be α = 1.2, which gives an expected IC photon index Γ IC = α + 1 = 2.2.When the photon index is fixed to 2.2 (1T + IC), the resulting statistics (C-stat/ν = 6498.18/5981) is improved with respect to a 1T model, but it is worse than either a 2T model or a 1T + IC model with a free Γ parameter.
When Γ is fixed to the expected value of 2.2, we obtain an upper limit on the 4 − 25 keV IC flux of 1.5 × 10 −11 erg s −1 cm −2 .As discussed in section 4.5, this constrains the magnetic field to be B ≥ 0.35 µG.

Testing other possibilities
None of the standard models discussed thus far (1T, 2T, and 1T +IC) can both adequately describe the data and be physically motivated.To further investigate the origin of the hard excess, we tested three additional models.1T * Comptonization: Some plasmas may be able to Comptonize their own thermal bremsstrahlung emission, which could produce a hard excess.This phenomenon is modeled by convolving an apec with XSPEC's simpl model (Steiner et al. 2009).While we don't expect to see this phenomenon to occur in the diffuse ICM, a Comptonization model can mimic a non-Maxwellian scenario.
In models of collisionally ionized plasma, the base simplifying assumption is that collisions result in a Maxwellian distribution of electrons.However, simulations by Akahori & Yoshikawa (2010) show that cluster mergers can disturb the thermal equilibrium of the plasma.Relaxed clusters such as Perseus are not known to exhibit non-Maxwellian behavior, but we tested whether a similar model would be able to explain our excess.
We applied this Comptonization model to both the inner and outer regions of diffuse emission.We fixed the UpScOnly parameter to 0, which specifies that photons are only up-scattered.This choice reflects the fact that the relativistic electrons responsible for the Comptonization in this scenario are at much higher energies than the thermal photons they are scattering, and so we approximate that no photons are down-scattered.We linked the other parameters (the fraction of scattered photons and Γ, which is the photon index of the scattered emission) across both regions.With C-stat/ν = 6477.73/5982,this self-Comptonization model slightly improved the statistics with respect to the 2T model, but it was disfavored over 1T+IC with a free Γ.

Non-Maxwellian:
Because the Comptonization model-which mimics a non-Maxwellian scenario-fit the data fairly well, we implemented a more physically accurate non-Maxwellian model in search of a more realistic description of the data.PyKappa3 , which can be imported into PyXSPEC as a custom model, implements the findings of Hahn & Savin (2015) to model a non-Maxwellian equivalent of the APEC.
Fitting PyKappa to data is computationally expensive, so we fit a single apec to the central thermal region and a kappa model to the global annulus excluding the inner circle (Figure 1) in order to speed up fitting.This model slightly improved the statistics with respect to the 1T model but was otherwise disfavored compared to other models, with a C-stat/ν = 6536.28/5983.The kappa parameter-which describes the degree of non-Maxwellianness (see Cui et al. 2019)-had a value of 23.5 +4.7 −5.1 , which implies only a slight deviation from Maxwellianness.1T + Gaussian: We find that the hard excess was best described by a wide Gaussian (gauss) added to the 1T model.Applying the Gaussian to the model describing the AGN emission (powerlaw + gauss; Cstat/ν = 6527.02/5980)only slightly improves the overall fit compared to a 1T model.However, adding the Gaussian component to the model describing the ICM (apec + gauss, C-stat/ν = 6471.20/5980)made a significant improvement, and it was statistically favored over all of our other models.The Gaussian was initialized at 20 keV with a width of 1 keV to roughly match the observed properties of the hard excess.After fitting, the centroid became 17.7 keV, and the width became 3.6 keV.

Applying Crossarf and Background Systematics
As described in section 3.1, we use nucrossarf to extract spectra and generate cross-ARFs for three sources: a point-source at the center of the cluster (representing the AGN), an extended source in a small (r = 1') circular region, and a second extended source in a concentric annulus (r in = 1' and r out = 7', see Figure 1).
In the Appendix, we find that nucrossarf is able to model PSF leakage from point sources with a 3.4% systematic uncertainty within similarly sized regions and for a similar bandpass.This uncertainty governs how well we can model the AGN emission being scattered into the outer region, where it could be confused for diffuse non-thermal emission.Because the modeled AGN light scattered into the outer region contributes 25% of the total 20-25 keV ICM emission, it would need to be 39% stronger than has been predicted by nucrossarf in order to explain the 10% excess, which far exceeds the 3.4% uncertainty.However, to thoroughly account for this uncertainty in our models, we followed a similar process that Wik et al. (2014) used to characterize the systematic error of the background components; we applied this process to the AGN cross-ARF as well as to the background spectra.
Following this method, we generated 1000 realizations of the AGN cross-ARF.We varied the normalization of the cross-ARFs with a scatter of 3.4%.Similarly, we also simultaneously generated 1000 realizations of the background model, assuming a Gaussian distribution (with σ equal to the systematic error) about the nominal value for each component.For each iteration, we fit four of our models [1T, 2T, 1T+IC (free Γ), and 1T+IC (fixed Γ)] with the modified background and AGN cross-ARF, finding new values of C-stat for each model, that we compare relative to each other for each iteration.The results are shown in Figure 7.

Spatial distribution of the hard excess
In order to characterize the spatial distribution of the hard excess, we performed the following analysis for each NuSTAR pixel within a box with the range 50.107 • < α < 49.796 • , 41.395 • < δ < 41.628 • , which contains the entire circular global region shown in Figure 1.
First, we centered a circular region (r = 10 pixels) around the pixel.Using the soft-band (3-10 keV) image (background and AGN subtracted), we calculated the total number of soft counts within the region.We repeated this process for the hard-band (15-25 keV) image.Meanwhile, we placed the circle region in the NuSTAR temperature map (Figure 5) and took the mean pixel value as the average temperature within that region.
For each pixel, we then used the NuSTAR response files (for both FPMA and FPMB ) and the nearest ARFs to that pixel for each observation (found in the nucrossarf libraries) to simulate an APEC model with kT set to the the average temperature found within the r = 10 pixel region.
From the APEC model, we extracted the flux count ratio between the soft (3-10 keV) and hard (15-25 keV, where the excess is strongest) energy bands.If all the hard-band counts from the ICM are due to thermal emission, this modeled soft-to-hard ratio will be equal to the actual soft-to-hard count ratio of the data.So, by dividing the soft-band counts in the data by the modeled soft-to-hard ratio, we predicted the thermal counts in the hard (15-25 keV) band.We then subtracted this number from the actual hard-band data to obtain a map of the hard excess.These steps are summarized by the images in Figure 8.By performing the above process for each pixel, we obtained an image of the hard excess and calculated its signal-to-noise ratio using the following equation: Using this procedure, we also mapped the spatial distribution of the hard excess in each individual observation.In Figure 9, we show the discrepancy between these three observations, which have roll angles of 89.7 • (90202046002 and 90202046004) and 323.5 • (60061361002).The spatial distribution of the excess  varies considerably amongst the pointings, and correcting for the roll angle does not eliminate the discrepancy significantly.
In the light of this variation, we then fit the spectra from each individual observation to a 1T + Gauss model in order to check for the presence and significance of the excess on an observation-by-observation basis.As shown in Figure 10, the excess is only significant at the 3σ level in observation 90202046004, and is not detected at even the 90% level in observation 60061361002.

Spectral modeling
As described in section 3.4, we fit several models to our data in order to unveil the physical nature of the ICM in the Perseus cluster.The statistics and best-fit parameters of these models are summarized in Table 2.
Fitting the ICM spectra to a single-temperature APEC (1T) yields a global temperature of 5.45 +0.04 −0.03 keV.This is consistent with the temperature found using Chandra (∼5.5 keV at 120 kpc.See Figure 5 and Sanders & Fabian 2007;Schmidt et al. 2002).This 1T model describes the data well, but when additional components are added to the spectra, their best-fit parameters reveal a hard excess at ∼ 15−25 keV.
When a second APEC component is added to describe the multi-temperature nature of the ICM, the best fit second temperature reaches nonphysically high values and hits the default ceiling for the parameter (64 keV).At these high temperatures, the turn-over for bremsstrahlung emission occurs outside of the analyzed energy band, causing the second APEC component to mimic a powerlaw.We conclude that the hard excess can not be described by thermal emission.
When we included a powerlaw component to account for possible non-thermal contribution to the ICM spectra (1T + IC), the best-fit model was not consistent with IC.From the radio synchrotron index found by Gendron-Marsolais et al. ( 2021), the IC photon index should be Γ = 2.2.But the best-fit powerlaw is extremely flat (Γ = 0.254), and fixing the powerlaw index to the expected value (2.2) was statistically disfavored over a powerlaw with a free Γ.This implies that the hard excess is not IC emission, which is consistent with expectations.The Perseus Cluster has a radio minihalo (discovered by Burns et al. 1992), implying that the relativistic electron population is diffuse in radio mini-halos.The radio synchrotron flux of these objects is notoriously faint, so the IC signal from these relativistic electrons should also be relatively weak.With a synchrotron flux density of 12.64 Jy (230-470 MHz, see Gendron-Marsolais et al. 2017), the Perseus Cluster's mini-halo is brighter than most by 3-4 orders of magnitude.However, if the global magnetic field is ∼ 50µG, as estimated by Taylor et al. (2006), the IC flux would be 3 × 10 −16 erg s −1 cm −2 , which is five orders of magnitude fainter than the best-fit IC model (where Γ is fixed).So, although the hard excess is not IC, that finding is consistent with our expectations.
We fit the hard excess to three additional models that could describe a hard excess: a PyKappa model, which is the Non-Maxwellian equivalent of the APEC, a Comptonization model (expressed by a SIMPL convolution), and a Gaussian hump.Of these three, only the Gaussian was favored over the 1T, 2T, and 1T+IC models.
Originally, the Gaussian model was motivated by the Compton Hump that occurs in obscured AGN spectra, but Rani et al. (2018) find no evidence of reflection components (including the Compton Hump) in the AGN of the Perseus Cluster, so we did not expect to find a significant improvement.As expected, applying the Gaussian to the model that describes the AGN is only slightly favored over the 1T model presented in Table 2, where the AGN is described only by a powerlaw.However, adding the Gaussian to the model describing the ICM made a significant improvement.Although this component does not describe a known physical phenomenon, it provides useful information about the nature of the excess.Unlike a 1T+IC or 2T model, the Gaussian centered at hard energies has negligible contribution to soft-band emission.Similarly, the best-fit flat powerlaw and 64 keV 2T models are trying to fit the excess while minimizing their effect on the soft-band spectrum.This is demonstrated in Figure 6, which compares the two models that best describe the excess (1T + Gaussian Where the lower bound dropped to zero, the normalization is represented as an upper limit.and 1T + IC with a free Γ) with the three main physically motivated models (1T, 2T, and 1T + IC with Γ fixed).The components that attempt to describe the excess converge near ∼ 17 keV, demonstrating that all four of these models are fitting to the same feature in the hard end of the spectrum.These results suggest that the source of the hard excess only makes meaningful contributions to the emission around 15-20 keV, while the lower-energy spectrum is well described by pure thermal emission.
The wide Gaussian that best fits our data resembles an unaccounted-for effect in the background.For example, a more recent analysis of archival data has revealed emission lines in the instrumental background that are not included in the current background model (Roach et al. 2023;Rossland 2022).Two of these emission lines are described by Lorentzians with µ = 13.07,11.40 keV and σ = 20.00,19.37 keV.Neither of these lines align perfectly with the best-fit Gaussian for the hard excess-their locations are softer and their σ values are wider than our Gaussian component-but fitting the excess to the harder of these two Lorentitzans (C-stat/d.o.f = 6494.68/5982) is statistically favored over the 1T and 1T + IC (fixed Γ) models.It is disfavored over the 2T, 1T + IC (free Γ), and 1T + Gauss models, so we cannot conclude that this new background feature explains the excess.
In addition to the Perseus cluster, we find evidence of features similar to the hard excess in other relaxed clusters, although the significance of the excess varies.Investigating these additional clusters is beyond the scope of this paper, but for future work, investigations of the hard excess would benefit from future X-ray telescopes with broader energy capabilities and lower background, such as the proposed High Energy X-ray Probe (HEX-P Madsen et al. 2019); reliable data above 25 keV would allow for better characterization of the hard excess and test if it is truly represented by a Gaussian.Additionally, the low background of HEX-P would better constrain the nature of the excess and whether it is an effect of the sky, an intrinsic property of the ICM, or an instrumental effect of NuSTAR.

Implications of cross-ARF and background systematics
In Figure 7, we show the results of applying the systematic errors in the nucrossarf and the NuSTAR background models to our physically motivated models of the ICM (1T, 2T, and 1T + IC with both fixed and freed Γ).In this analysis, we generated 1000 iterations of the background models and AGN cross-ARFs, varying their contribution to the spectra within their systematic errors.The right-most plot in the figure shows the difference in the C statistic (∆ C-stat) between a 2T and a given model (for ∆ C > 0, the 2T model is favored).This plot demonstrates that the 2T model is consistently favored over a 1T or physically realistic IC model (Γ fixed at 2.2), while it did worse than a powerlaw with a free Γ.Additionally, the second temperature in the 2T model hit the 64 keV ceiling in almost every case.
As shown in the middle panel of Figure 7, ∆ C-stat does not correlate with the strength of the AGN's contribution (within the 3.4% systematic error).However, there is an obvious correlation with the aperture stray light component of the background (S d ).Increasing the value of the S d component caused our models to converge to a similar C-stat, while decreasing the contribution of S d widened the disparity between models and increased the significance of the hard excess.At the same time, increasing S d beyond its nominal value worsened the fit for all models, which is strong evidence that increasing the S d does not describe the actual background.
We conclude that the hard excess cannot be explained by systematic uncertainties in either the background or AGN contribution, and its characteristics remain consistent regardless of changes to the background or AGN contribution.

Variation of the hard excess
Figure 8 shows the AGN and background-subtracted image of excess photons in the 15 − 25 keV energy band; these photons are not well described by thermal emission, and we take them to be responsible for the hard excess.There are regions of excess photons with signalto-noise > 3, which suggests that the hard excess may be localized.However, the spatial distribution of the excess photons differ across observations, as shown in Figure 9.These variations occur on kpc scales; if the excess were an attribute of the ICM, it should remain static over the time these observations were taken.If the excess is due to an unaccounted-for detector effect, it may be correlated with the location of the detector chips.However, the roll angle correction shown in the bottom-right plot of Figure 9 does not significantly improve the discrepancy between observations.
In Figure 10, we show that the strength of the excess (when represented by a Gaussian) is inconsistent between observations at the 90% level.It is undetected at the 90% level in observation 60061361002, which has the lowest exposure time of the three observations, and is only detected at the 3σ level in observation 90202046004, which has the highest exposure time.To investigate the chance that the hard excess is present-simply not detectable-in observation 60061361002, we present the signal-to-noise ratio of the excess (SNR excess = N excess / √ N total , where N excess is the modeled number of photons in the Gaussian model component and N total is the total number of events present, both in the 15-25 keV band) of our observations in Table 3.While SNR appears to increase with exposure time, the difference in SNR excess between observations 90202046002 and 90202046004 far exceeds what we would expect given the 0.3 ks difference in exposure times.While this does not completely rule out the possibility that the excess is a physical property of the ICM, it is more in line with expectations for a transient background feature or other instrumental effect.

Thermal distribution
Figure 5 shows the NuSTAR and Chandra temperature map of the ICM of the Perseus Cluster.The comparison shows that the temperature structure is consis-tent between both telescopes, with a cool core and hot (∼ 9 keV) gas outside the core.
The main discrepancy between the maps is that NuS-TAR measures higher temperatures across the cluster.It is a long-standing issue that different instrumentsincluding NuSTAR and Chandra-measure different temperatures when observing the same ICM.The calibration differences between these telescopes are not well understood, though work is being done to quantify them (Potter et al. 2023).In the simplest form, the discrepancy that we see between the measurements in Figure 5 is consistent with our expectations; NuSTAR's harder bandpass biases it towards measuring higher temperatures when a multi-temperature gas is present in a spectrum being fit to a single apec model.
Both NuSTAR and Chandra clearly observe the Perseus Cluster's cool core, and the cold front-the sharp discontinuity between the cool core and the hot outer ICM (first discovered by Fabian et al. 2011)-is apparent in both temperature maps.

Magnetic field calculation
In Section 3.4.3,we obtained an upper limit to the IC flux of F IC(4−25 keV) ≤ 1.5 × 10 −11 erg s −1 cm −2 .From this, we can constrain the magnetic field using the following equation as derived by Govoni & Feretti (2004): where the energy range of F IC(E1−E2) is 4-25 keV, α is the spectral index of the radio synchrotron emission, F sync(νr) is the synchrotron flux at frequency ν r , and h(α) is tabulated in Table 2 of Govoni & Feretti (2004).Gendron-Marsolais et al. (2021) report that the synchrotron flux of the radio mini-halo is 12.64 Jy (230-470 MHz; for the following calculation, we use the midpoint frequency of ν r = 350 MHz) and find the sychrotron index to be α = 1.2, which gives h(α) = −1.37 × 10 −14 .Folding these values into Equation 3, we constrain the magnetic field to be ≥ 0.35 µG.Although this measurement is consistent with previous estimates of the magnetic field in the Perseus cluster, it is an order of magnitude lower than the lower limit found by Aleksić et al. (2012) and two orders of magnitude lower than the the estimate by Taylor et al. (2006).Our lower limit is consistent both with these estimates and with the expectation that the magnetic field in the cool cores of clusters should be stronger than the globally averaged value of B; NuSTAR measurements of other clusters yield similar lower limits (Rojas Bolivar et al. 2023), so a detection of IC emission in the core of Perseus at this level would be unexpected.

CONCLUSIONS
• We studied the thermal structure of the Perseus cluster using three NuSTAR observations, and the results are consistent (within known systematic effects) with Chandra temperature measurements, as shown in Figure 5.
• We find evidence of a hard excess in two of the three observations that cannot be well modeled by thermal emission or explained by a known systematic error in the background or AGN cross-ARF modeling.The excess is best described by a wide Gaussian at 17 keV, and it is not consistent with a realistic inverse Compton model.
• The strength and spatial distribution of the hard excess is inconsistent across observations, which possibly points to a non-astrophysical origin.
• By applying a physically motivated IC model to the excess, we constrain the cluster's global magnetic field to be ≳ 0.35 µG, which is 1-2 orders of magnitude lower than (and consistent with) previous measurements made by Aleksić et al. (2012) and Taylor et al. (2006).
• We report that nucrossarf models scattered light from point sources with an accuracy of 3.4%, which we report as its systematic uncertainty (see Appendix A for details).a concentric annulus from 30 ′′ to 3 ′ .After obtaining ARF and RMF files for each region, we loaded the data from the inner region into XSpec and fit the AGN with all model parameters freed.
In most cases, we fit the AGN to a powerlaw model.An exception was the blazar Mrk 421.For observations of this well-studied object, we instead fit the spectra to a log-parabola model with a pivot energy fixed at 5 keV, as described by Baloković et al. (2016).
Once the inner region's spectrum was fit, we loaded the spectrum and response files for the outer region.Ideally, the extrapolated model in the outer annulus should match the data without re-fitting.Discrepancies would stem from inaccurate cross-ARFs, which would be due to systematic uncertainties in the PSF modeling employed by nucrossarf.
We calculated the systematic uncertainty using the equation where E N is the normalized error and N model and N data are the modeled and actual number of photons in the outer annulus, respectively.If all sources of error are accounted for, the distribution of E N for many observations should be a Gaussian with µ = 0 and σ = 1.We considered two sources of error in our observations: statistical error (σ stat = √ N data ) and systematic error (σ sys ).We solved this equation on our data by finding a value of σ sys that produced a distribution of E N with µ = 0 and σ = 1 (Figure 12).To do this, we used the number of backgroundsubtracted photons from the outer annulus (N data ) and the number of photons predicted by the model (N model ) and made an initial estimate for a value of σ sys to calculate E N for each source.Once the distribution of E N was obtained, we fit it to a Gaussian model.When the standard deviation of this distribution was greater/smaller than 1.0, we increased/decreased the value of σ sys .Following this treatment, we used the new σ sys to repeat the process until the standard deviation of E N was within 1.000 ± 0.005.Once E N = 1.000 ± 0.005, we extracted the corresponding value of σ sys .Here, we report the systematic uncertainty as a relative standard error σ/N data × 100.However, depending on the methods used to fit our data, we obtained two different values of the systematic error.In method one, we froze all model parameters after fitting to the inner region, and we obtained σ sys based on the pure extrapolation into the outer region.This yielded a systematic uncertainty of σ sys /N × 100 = 6.1%.
In method two, the normalization parameter was freed (and tied between the inner and outer regions) to vary while other parameters were fixed, and re-fit the model to our outer and inner regions simultaneously, yielding a systematic uncertainty of σ sys /N × 100 = 3.4%.This method more accurately reflects the way that nucrossarf results are modeled, so we chose to use 3.4% as the systematic error for the this analysis.This reflects the accuracy of nucrossarf to consistently model a point-like source between a typical circular region and a large outer annulus that captures emission from the wings of the PSF.
The E N distribution we obtained with σ sys /N × 100 = 3.4% is shown in Figure 12.The distribution does not have significant skewdness (0.21) or kurtosis (-0.34).
The exact level of uncertainty should be assessed on a case by case basis, as it may depend on the size and shape of spectral extraction regions and the nature of the source emission (point-like versus diffuse).Additionally, future work is required to characterize the uncertainties in different energy bands.

Figure 1 .
Figure 1.The combined background-subtracted, exposurecorrected, and AGN-subtracted image of Perseus in the 4-25 keV band, smoothed with a Gaussian (σ = 4 pixels).The black circles are the spectral extraction regions.

Figure 2 .
Figure 2. The combined FPMA and FPMB images from each of the three observations of the Perseus cluster.Shown in white are the extraction regions used to fit the background across the FOV with labels denoting the corresponding detector chip, and the black hashed circle is the exclusion region.

Figure 3 .
Figure3.The background data and best-fit model for observation 90202046004.The spectra and models corresponding with FPMA detectors 1,2, and 3 (see Figure2) are black, red, and green, respectively, and the three detectors on FPMB are blue, cyan, and magenta, respectively.

Figure 4 .
Figure 4. Single temperature fits to spectra from the two regions shown (inset), including an power law AGN model at the center of the inner region.Top: The data (crosses, with a different color for each individual spectrum) plotted with the different components of the model (solid lines, color-matched to the corresponding data) with components labeled.The higher-flux data belongs to the global outer annulus, and the lower flux data belong to the inner circular region (Fig. 1).Bottom: The ratio of the data compared to the total model.The yellow box indicates the appearance of the hard excess.

Figure 5 .
Figure 5.Comparison between the temperature map produced using NuSTAR data (left) and Chandra data (right, from Sanders & Fabian 2007).NuSTAR temperature contours are overlaid in both plots for easy comparison.

Figure 6 .
Figure 6.Top panel: Comparison of the components of four different models that could describe the hard excess in the Perseus cluster (solid lines are the total model, and the dashed lines of the same color are the individual components of the model): 1T, 2T, 1T + IC (free photon index), 1T + IC (fixed photon index), and a Gaussian.The total model is plotted by solid lines while components are dashed.The dotted black line gives the AGN's contribution to the model.Bottom five panels: the ratio of the models to the data, demonstrating the goodness of fit across the energy range.The C-stats and best-fit parameters are shown in Table 2 photons keV −1 cm −2 s −1 at 1 keV)

Figure 7 .
Figure 7.The difference in the C-statistic (∆ C-stat) between the 2T model and 3 other models: a 1T (cyan; triangles), an IC model with a fixed photon index (pink; circles), and an IC model with a free powerlaw index (purple; squares) for 1000 realizations of the background model and scattered AGN contribution.The 2T model is favored when ∆ C-stat > 0. Left: The histograms of ∆ C-stat, where the dashed line marks ∆ C-stat = 0. Middle: ∆ C-stat as a function of the AGN scattering factor, which shifts the contribution of the AGN relative to the nominal value (the AGN contribution is strengthened for a factor > 1).Right: ∆ C-stat as a function of the S d (aperture stray light) factor, which shifts the strength of the S d background component relative to its nominal value (S d is strengthened for a factor > 1).See section 3.5 for discussion.

Figure 8 .
Figure 8. Left: The AGN and background-subtracted counts from 15-25 keV for the three summed observations.Middle-left: The modeled thermal emission in the hard (15-25 keV) band based on the procedure outlined in Section 3.6.Middle-right:The thermal-subtracted hard-band counts.Right: The SNR between the data, representing the local significance of the hard excess.All images are smoothed by the binning process described in section 3.6, and the radius of the smoothing kernel is r = 10 pixels.

Figure 9 .
Figure9.Top: The excess in each observation (from left to right: 60061361002, 90202046004, and 90202046002) in units of counts, where the rounded extrusions along the edge of the observations are artifacts of the smoothing function.The locations of the four detector chips are divided by the black cross-hatch and labeled with white numbers.The color of the circle surrounding the AGN gives the summed excess in counts/pixel within a sector centered on that angular bin.If the excesses originated from the cluster emission, it should be correlated between observations.Bottom: The excess as a function of angle for all three observations, either with respect to the sky (left) or the detector plane (right).The excess with respect to angle on the sky shows clear discrepancies between observations, as it also does with respect to roll angle.

Figure 10 .
Figure10.The normalization of the Gaussian component, which is defined as the total photons cm −2 s −1 contained within the line.This component represents the hard excess in the spectra, with 90% (purple) and 3σ (black) error bars.Where the lower bound dropped to zero, the normalization is represented as an upper limit.

Figure 11 .
Figure 11.The observations used to estimate the accuracy of cross-ARF modeling.The observation IDs are listed along the top of the background-subtracted images; see table 4 for more information.

Figure 12 .
Figure12.The normalized error EN (see Eq. A1) for 25 observations in both FPMA and FPMB (50 total) of point sources in the NuSTAR archive.The solid black line shows an illustrative probability distribution function of a Gaussian with µ = 0 (dashed black line) and σ = 1, which is the distribution of EN when all sources of error are accounted for.The pink histogram is the distribution of EN when there is no systematic error accounted for.The purple histogram includes a 3.4% systematic error, which produces the desired distribution of EN.The data appear to be negatively skewed, which would indicate that nucrossarf tends to over-model point source emission, but the skewdness (0.21) is not significant.

Table 1 .
Observations of the Perseus Cluster

Table 3 .
Signal to Noise Ratio of Excess

Table 4 .
Observations Log For Crossarf Calculation Errors