<figure><img src="image_1.png" /></figure>

# Spectral Evolution of the Ultraluminous X-Ray Sources M82 X-1 and X-2

Murray Brightman¹, Dominic J. Walton² ID, Yanjun Xu¹ ID, Hannah P. Earnshaw¹ ID, Fiona A. Harrison¹, Daniel Stern³ ID, and Didier Barret⁴,⁵ ID

¹Cahill Center for Astrophysics, California Institute of Technology, 1216 East California Boulevard, Pasadena, CA 91125, USA

²Institute of Astronomy, Madingley Road, Cambridge CB3 OHA, UK

³Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA

⁴Université de Toulouse, UPS-OMP, IRAP, Toulouse, France

⁵CNRS, IRAP, 9 Av. colonel Roche, BP 44346, F-31028 Toulouse cedex 4, France

Received 2019 September 18; revised 2019 December 5; accepted 2019 December 15; published 2020 January 27

## Abstract

M82 hosts two well-known ultraluminous X-ray sources (ULXs). X-1, an intermediate-mass black hole (IMBH) candidate, and X-2, an ultraluminous X-ray pulsar. Here, we present a broadband X-ray spectral analysis of both sources based on 10 simultaneous observations made with *Chandra* and *NuSTAR*. *Chandra* provides the high spatial resolution to resolve the crowded field in the 0.5–8 keV band, and *NuSTAR* provides the sensitive hard X-ray spectral data, extending the bandpass of our study above 10 keV. These observations, taken in 2015–2016, cover a period of flaring from X-1, allowing us to study the spectral evolution of this source with luminosity. During four of these observations, X-2 was found to be at a low flux level, allowing an unambiguous view of the emission from X-1. We find that the broadband X-ray emission from X-1 is consistent with that seen in other ULXs observed in detail with *NuSTAR*, with a spectrum that includes a broadened disk-like component and a high-energy tail. We find that the luminosity of the disk scales with inner disk temperature as $L \propto T^{-3/2}$, contrary to expectations of a standard accretion disk and previous results. These findings rule out a thermal state for sub-Eddington accretion, and therefore do not support M82 X-1 as an IMBH candidate. We also find evidence that the neutral column density of the material in the line of sight increases with $L_X$, perhaps due to an increased mass outflow with accretion rate. For X-2, we do not find any significant spectral evolution, but we find the spectral parameters of the phase-averaged broadband emission are consistent with the pulsed emission at the highest X-ray luminosities.

*Unified Astronomy Thesaurus concepts: Astrophysical black holes* (98); *Intermediate-mass black holes* (816); *Stellar mass black holes* (1611); *X-ray point sources* (1270); *X-ray sources* (1822); *Pulsars* (1306); *Neutron stars* (1108)

## 1. Introduction

Ultraluminous X-ray sources (ULXs), first discovered by the *Einstein X-ray observatory* (Giacconi et al. 1979; Fabbiano 1989), are observed as bright sources of X-rays that appear to exceed the Eddington limit of the typical $10 M_{\odot}$ black holes found in our own Galaxy. ULXs are not coincidental with the nuclei of their host galaxies, so they are not be powered by the supermassive black holes that power active galactic nuclei (AGNs), although ~25% of candidate ULXs are likely to be background AGNs (e.g., Walton et al. 2011; Sutton et al. 2015; Earnshaw et al. 2019). At first, ULXs were promising intermediate-mass black hole (IMBH) candidates, as larger black holes can radiate at higher luminosities due to the Eddington limit scaling with mass (e.g., Colbert & Mushotzky 1999; Miller et al. 2003). But as more observations were made, the data did not appear consistent with this scenario, and instead, lower-mass black holes accreting at super-Eddington luminosities were favored for most sources (e.g., Mizuno et al. 1999), with a few IMBH candidates still remaining (e.g., ESO 243-49 HLX-1, Farrell et al. 2009).

M82 (the “cigar galaxy,” or NGC 3034; Watson et al. 1984) hosts two well-known ULXs that, as they are separated by only 5" on the sky, were only first resolved by *Chandra* (Matsumoto et al. 2001). The history of studies of M82 X-1 (CXOU J095550.2+694047), the brightest ULX in M82, embodies the above narrative. It has long been considered one of the best IMBH candidates because of its high luminosity, which can

reach ~$10^{41}$ erg s$^{-1}$ (e.g., Ptak & Griffiths 1999; Rephaeli & Gruber 2002; Kaaret et al. 2006); detection of low-frequency quasi-periodic oscillations (QPOs) in the power spectrum (54 mHz, Strohmayer & Mushotzky 2003; Dewangan et al. 2006; Mucciarelli et al. 2006), indicative of a compact, unbeamed source; as well as twin-peaked QPOs at 3.3 and 5.1 Hz, which lead to a mass estimate of 400 $M_{\odot}$ using scaling laws between the QPOs frequencies and mass used for stellar-mass black holes (Pasham et al. 2014). Additionally, Feng & Kaaret (2010) (hereafter F10) observed the source with *XMM-Newton* and *Chandra* over the course of a flaring episode, and fitted the spectra with the standard thin accretion disk model. They found that the luminosity of the disk, $L$, scaled with inner temperature as $L \propto T^4$, which is expected from a thin accretion disk with a constant inner radius. From this, they inferred a black hole mass in the range 300–810 $M_{\odot}$, assuming that the black hole is rapidly spinning to avoid extreme violations of the Eddington limit, therefore adding support to the IMBH scenario.

The luminosity argument no longer stands, however, as the ULX NGC 5907 ULX1, which also reaches ~$10^{41}$ erg s$^{-1}$ (Sutton et al. 2013; Fürst et al. 2017), and was once considered an IMBH candidate, was discovered to be powered by a neutron star with only 1–2 $M_{\odot}$ (Israel et al. 2017) from the detection of coherent pulsations. The mass measurement from the twin QPOs is ambiguous too, as it is not known where these originate, and it is not clear if the scaling relationship used
extends to the IMBH range. In Brightman et al. (2016b, hereafter B16), we presented a combined spectral analysis of X-1, also during a flaring episode, using simultaneous observations with *Chandra*, *NuSTAR* and *Swift*. With broader-band data than Feng & Kaaret (2010), we found that the temperature profile as a function of disk radius ($T(r) \propto r^{-p}$) is significantly flatter than expected for a standard thin accretion disk as implied by F10, and instead characteristic of a slim disk that is expected at high accretion rates. As only one observation was analyzed, the $L \propto T^4$ relationship could not be tested. Nevertheless, the mass estimates inferred from the inner disk radius for this model imply a stellar-remnant black hole ($M_{BH} = 26_{-6}^{+9} M_{\odot}$), when assuming zero spin, or an IMBH ($M_{BH} = 125_{-30}^{+45} M_{\odot}$) when assuming maximal spin.

M82 X-2 (CXOM82 J095551.1+694045) is typically the second brightest X-ray source in the galaxy, and was the first ultraluminous X-ray pulsar (ULXP) discovered (Bachetti et al. 2014) using NuSTAR (Harrison et al. 2013). With only a 5" angular separation from X-1, studying the spectral properties of this source in detail has been limited to *Chandra* observations. We conducted a study of the spectral and temporal properties of M82 X-2 in Brightman et al. (2016a), finding that the source's luminosity varies over two orders of magnitude over the range $10^{38}-10^{40}$ erg s$^{-1}$. Using timing analyses, we were able to isolate the pulsed emission from this source with NuSTAR, finding that it was well described by a cutoff power law. In a follow-up study, we found evidence that the variations in luminosity are modulated on a ~60 days period, and that as the orbit of the neutron star and its companion is known to be 2.5 days, the ~60 days period must be super-orbital in origin (Brightman et al. 2019).

In this paper, we present analysis of a systematic monitoring campaign on M82 by *Chandra* that took place in 2016. The primary goals of this campaign were to perform a temporal analysis of X-1 and X-2 to search for orbital and super-orbital modulations; to perform spectroscopic studies of the ULXs; and to study the nature of the other binary systems in M82. We report the temporal analysis in Brightman et al. (2019), and here we focus on the second objective: to present the most comprehensive X-ray spectral analysis of the two ULXs, M82 X-1 and X-2, to date. We combine simultaneous observations with *Chandra* to spatially resolve the two sources from each other and the other sources in M82, and simultaneous observations with NuSTAR to extend the spectral coverage up to 79 keV. We also present results from *Swift* monitoring observations of M82 that have been ongoing since 2012. Throughout this paper, we assume a distance to M82 of 3.3 Mpc (Foley et al. 2014), which is inferred from the light curve of SN2014J.

## 2. Observations and Data Reduction

All *Chandra* and NuSTAR observations studied here were taken simultaneously, or quasi-simultaneously (within 24 hr of each other). Table 1 provides a description of the observational data. The following sections describe the individual observations and data reduction.

### 2.1. Chandra

As the angular separation of X-1 and X-2 is only 5", only *Chandra* (Weisskopf 1999) can spatially resolve the emission from these two sources. The majority of the *Chandra* data

<table><caption>**Table 1**
Observational Data</caption><thead><tr><th>Observatory</th><th>Start Date</th><th>Start Time (UT)</th><th>ObsID</th><th>Exposure (s)</th></tr></thead><tbody><tr><td>NuSTAR</td><td>2015 Jan 15</td><td>21:41:07</td><td>50002019002</td><td>31249</td></tr><tr><td>Chandra</td><td>2015 Jan 16</td><td>13:40:00</td><td>17578</td><td>10070</td></tr><tr><td>NuSTAR</td><td>2015 Jun 20</td><td>14:21:07</td><td>90101005002</td><td>37409</td></tr><tr><td>Chandra</td><td>2015 Jun 21</td><td>02:45:09</td><td>17678</td><td>10100</td></tr><tr><td>Chandra</td><td>2016 Jan 26</td><td>19:44:49</td><td>18062</td><td>25100</td></tr><tr><td>NuSTAR</td><td>2016 Jan 26</td><td>20:06:08</td><td>80202020002</td><td>36136</td></tr><tr><td>NuSTAR</td><td>2016 Feb 23</td><td>17:01:08</td><td>80202020004</td><td>31670</td></tr><tr><td>Chandra</td><td>2016 Feb 24</td><td>00:37:06</td><td>18063</td><td>25100</td></tr><tr><td>NuSTAR</td><td>2016 Apr 5</td><td>09:41:08</td><td>80202020006</td><td>30505</td></tr><tr><td>Chandra</td><td>2016 Apr 5</td><td>16:04:41</td><td>18064</td><td>25090</td></tr><tr><td>NuSTAR</td><td>2016 Apr 24</td><td>19:16:08</td><td>80202020008</td><td>40357</td></tr><tr><td>Chandra</td><td>2016 Apr 24</td><td>20:02:13</td><td>18068</td><td>25100</td></tr><tr><td>NuSTAR</td><td>2016 Jun 3</td><td>20:41:08</td><td>30202022002</td><td>39022</td></tr><tr><td>Chandra</td><td>2016 Jun 3</td><td>22:10:58</td><td>18069</td><td>25100</td></tr><tr><td>NuSTAR</td><td>2016 Jul 1</td><td>16:26:08</td><td>30202022004</td><td>47043</td></tr><tr><td>Chandra</td><td>2016 Jul 1</td><td>23:18:08</td><td>18067</td><td>26100</td></tr><tr><td>Chandra</td><td>2016 Jul 29</td><td>07:50:01</td><td>18065</td><td>25100</td></tr><tr><td>NuSTAR</td><td>2016 Jul 29</td><td>23:06:08</td><td>30202022008</td><td>42846</td></tr><tr><td>NuSTAR</td><td>2016 Oct 7</td><td>20:21:08</td><td>90202038002</td><td>45476</td></tr><tr><td>Chandra</td><td>2016 Oct 8</td><td>00:39:30</td><td>18070</td><td>23200</td></tr></tbody></table>

analyzed here were taken during 2016 (Cycle 17) as part of a large program aimed at systematic monitoring of binaries in M82. The program consisted of 12 individual observations taken at ~monthly intervals, but only eight having simultaneous NuSTAR observations were used in this work. We additionally use two observations taken in 2015 that also have a simultaneous NuSTAR observation. Full details are listed in Table 1. All 2016 observations were taken with ACIS-I at the optical axis using a 1/8 subarray of pixels on chip I1 or I3, depending on the roll angle. The ULXs at the center of M82 were placed 3.5 off-axis to smear out the point-spread function (PSF) enough to reduce the effects of pileup, but not so much as to cause significant blending of the PSFs. The subarray of pixels was used to decrease the readout time of the detector to 0.4 s, further reducing the effects of pileup.

We proceeded to extract the *Chandra* spectra in the same way as described in B16, using the CIAO (v4.7, CALDB v4.6.5) tool SPECEXTRACT. For the point sources, spectra were extracted from elliptical regions drawn by eye to encompass the shape of the off-axis *Chandra* PSF. For X-1, we used an ellipse with a semimajor axis of 2"-3" and a semiminor axis of 1"-2". For X-2, the major and minor axes were 2" and 1", respectively. A small rectangular region close by was used for background subtraction. Figure 1 shows examples of the extraction regions used.

We used the CIAO tool PILEUP_MAP to give an indication of the level of pileup in each observation. The output, which is in counts per frame, ranges from 0.06 to 0.35 at maximum, which occurs at the position of X-1. These numbers correspond to pileup fractions of <5% for <0.1 counts per frame to >10% for >0.2 counts per frame. For three observations, obsIDs 17678, 18065, and 18067, X-1 has a pileup fraction of >10%. We will discuss the potential effects of this on our results later in the paper.

For X-ray emission from M82 that does not come from X-1 or X-2, but that contributes to the emission seen by NuSTAR, such as from the fainter point sources and the diffuse emission, we extract spectra from a 49" radius circular region centered on
<figure><img src="image_2.png" /><figcaption>**Figure 1.** *Chandra* image of M82 from obsID 17678 showing examples of the various extraction regions used, including the background for *Chandra* analysis shown as a small rectangular region. Small ellipses are used for X-1 and X-2, whereas a large circle with a radius of 49" is used to extract events from the rest of the galaxy, and the *NuSTAR* and *Swift/XRT* events. The brightest point sources within this region are labeled, however X-1 and X-2 dominate. North is up, east is left, as indicated by the arrows in the upper right corner, which are 10" long.</figcaption></figure>

X-1, but with the X-1 and X-2 regions as described above masked out. A larger background region external to the galaxy was extracted to assess the background coming from the cosmic X-ray and particle backgrounds.

## 2.2. NuSTAR

The raw *NuSTAR* data were reduced using the NUSTARDAS software package version 1.4.1 and CALDB version 20150316. The events were cleaned and filtered using the *nupipeline* script with standard parameters. The *nuproducts* task was used to generate the spectra and the corresponding response files. Spectra were extracted from a circular aperture of radius 49" centered on the peak of the emission. The background spectra were extracted from a circular region encompassing the same detector chip as the source, with a radius of 118", excluding the source extraction region and avoiding the wings of the PSF as much as possible. Data from both focal plane modules (FPMA and FPMB) were used for simultaneous fitting, without coadding.

## 2.3. Swift

*Swift* conducted monitoring of M82 with a typical cadence of a few days between 2012 and 2018. As this monitoring ran contemporaneously with our *Chandra* and *NuSTAR* observations, the well-sampled light curve provides context for our study. A total of 113 observations, mostly consisting of obsIDs 00032503099–154 and 00092202001–051, have been made of the galaxy over the 2015–2016 period which we use here to calculate a long-term light curve.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** 0.5–8 keV light curve of the total observed X-ray flux from M82 using *Swift/XRT*. The highly variable X-ray emission is caused by X-1. The times of the simultaneous *Chandra* and *NuSTAR* observations we use here to study the spectral evolution of X-1 are shown with dashed lines.</figcaption></figure>

We calculate the fluxes via spectral fitting. We use the HEASOFT (v 6.16) tool XSELECT to filter events from a 49" radius region centered on the ULXs and to extract the spectrum. This extraction region encloses all sources of X-ray emission in the galaxy. Background events were extracted from a nearby circular region of the same size. We group the spectra with a minimum of one count per bin using the HEASOFT tool GRPPHA. We conduct spectral fitting in the range 0.2–10 keV. We fit the spectra with a simple power law subjected to absorption intrinsic to M82 at $z = 0.00067$ (zwabs*powerlaw* in XSPEC) with the Cash statistic (Cash 1979), which uses a Poisson likelihood function and is hence most suitable for low numbers of counts per bin. From this model we calculate the observed flux in the 0.5–8 keV range, equivalent to the *Chandra* band. The light curve is presented in Figure 2.

The *Swift* light curve shows that our observations took place during a period of flaring activity from M82, which we found to be due to increased activity from X-1 (B16). The first *Chandra/NuSTAR* observations took place before the increase in activity, and the remaining observations tracked the activity over the next two years.

## 3. Spectral Analysis

All *Chandra* and *NuSTAR* spectra were grouped with a minimum of 20 counts per bin using the HEASOFT tool *grppha*. Spectral fitting was carried out using XSPEC v12.8.2 (Arnaud 1996) and the $\chi^2$ statistic was used for spectral fitting to background subtracted spectra. All uncertainties quoted are 90%. We present the *Chandra* spectra of X-1, X-2, and the additional X-ray emission from M82, and the *NuSTAR* spectra of the entire galaxy for each of the 10 observational epochs in Figure 3.

While the presence of X-2 in a bright state introduced some ambiguity to the results on the spectral analysis of X-1 in B16, in four observations here, *Chandra* obsIDs 18062, 18065, 18069, and 18070, X-2 was observed in a low state, allowing an unambiguous view of X-1 for the first time using *NuSTAR*. For these four observations, we neglect the emission from X-2 in our spectral fits, whereas for the rest of the observations we included it.

We proceed to conduct a spectral analysis for all 10 observations, where we fitted the spectra in the same way as described in B16. A cross-calibration constant was applied to
<figure><img src="image_4.png" /><figcaption>**Figure 3.** *Chandra* spectra of X-1 (red), X-2 (green), and the diffuse emission (purple), and the *NuSTAR* FPMA spectra of all sources (blue) for each of the 10 observational epochs listed in Table 1. The spectra have been heavily rebinned for clarity.</figcaption></figure>

each spectral data set to allow for absolute differences in normalization, and allowed to vary by ±10% (Madsen et al. 2015). For the diffuse emission from M82, we use a combination of three absorbed *zwabs*+*apec* models with the temperatures and abundances fixed to the values found in B16. We allow the normalizations, both with respect to B16 and to each other, to vary here to account for small differences in the detector responses. For the spectrum of X-1, we use the *zwabs*+*diskpbb* model, where *zwabs* is a redshifted neutral absorption component and *diskpbb* is a model representing emission from a multicolour accretion disk, with a variable radial temperature profile. This model combination was found to best represent the emission from X-1 with regards to other disk models in B16. Additionally, as X-1 is a bright source, and despite measures taken to reduce pileup, the source still suffers from pileup. We account for this in spectral fitting using the *pileup* convolution model, with frame time set to 0.4 s (Davis 2001).

In the top panel of Figure 4, we show the residuals to this set of models which reveal strong residuals above 10 keV, especially in obsID 18062, which is reminiscent of the hard tail seen in the *NuSTAR* spectra of several other ULXs such as Holmberg II X-1 (e.g., Walton et al. 2015), including those already identified as neutron star accretors such as NGC 7793

P13 (Walton et al. 2018a) and NGC 5907 ULX (Fürst et al. 2017; Walton et al. 2018b). This was also seen in our analysis in B16, but as X-2 was bright during that observation, there was ambiguity regarding its origin. Here it is clear that it originates from X-1.

We test two models to fit this component, *simpl*, which describes the power-law emission from the Comptonization of a disk spectrum (Steiner et al. 2009), and a *cutoffpl* model used to model the pulsed emission from ULXPs (Brightman et al. 2016a; Walton et al. 2018a, 2018b). In Walton et al. (2018b), where no pulsations from a source had been detected, as is the case for M82 X-1, Γ was fixed at 0.5 for the *cutoffpl* model which is the average for the pulsed emission from ULXPs. Therefore we do the same when using the *cutoffpl* model and fix Γ = 0.5. The energy of the cutoff was allowed to vary. We find that the *simpl* model provides a better description of the data, with Δχ² = 70–200 for the addition of two free parameters from *simpl* and Δχ² = 3–90 for the addition of two free parameters from *cutoffpl* as shown in the middle and bottom panels of Figure 4. We do not find a significant improvement for the *cutoffpl* model if we allow Γ to vary. Therefore we use the *simpl* model to fit the high-energy emission from X-1 for the full data set.

We proceed to fit the *Chandra* and *NuSTAR* spectra of all 10 observational epochs with the *pileup*+*zwabs*+*diskpbb*+*simpl* model combination to describe the emission from X-1. For the six observations where X-2 is bright, we model the emission from this source with an absorbed cutoff power-law model, *zwabs*+*cutoffpl*, which was used in B16 to model the pulsed emission. Figure 5 presents the unfolded spectra with the different model components shown. Table 2 presents the best-fit spectral parameters for X-1, and Table 3 presents the best-fit spectral parameters for X-2.

To calculate the intrinsic luminosity of each source, we use the *cflux* model component in *xspec* placed *after* the absorption component with the normalization of the main model fixed. We calculate the intrinsic flux over the range 0.5–30 keV, and present these values with their uncertainties in Tables 2 and 3. We calculate the intrinsic luminosities over these ranges assuming a distance to M82 of 3.3 Mpc (Foley et al. 2014).

## 4. Spectral Evolution of M82 X-1

Our main goal is to explore the spectral evolution of M82 X-1. For each pair of spectral parameters, we compute Spearman's rank correlation, using the IDL tool `r_corre-late.pro`, to assess the presence of any correlation and its significance. This assesses how well the relationship between two variables can be described using a monotonic function. The results from this test, being the rank correlation coefficient, *ρ*, and the two-sided significance of its deviation from zero, *p*, are presented in Table 4.

The significance is a value in the interval [0.0, 1.0] and a small value indicates a significant correlation. The commonly used threshold to judge that a correlation is significant is *p* < 0.05. We find that this criterion is met for the following pairs of parameters: $L_X$ and $N_{\text{H}}$; $L_X$ and $T_{\text{in}}$; $L_X$ and $\log_{10}$ norm; $T_{\text{in}}$ and $\log_{10}$ norm; $L_X$ and $f_{\text{scatt}}$; $T_{\text{in}}$ and $f_{\text{scatt}}$; and Γ and $f_{\text{scatt}}$.

However, for many parameters, especially those of the *simpl* model, the uncertainties are large, which the correlation test does not account for. First, to account for this, we conduct 1000 Monte Carlo simulations where we randomly draw a
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Spectral residuals for a fit to the 18062 data set on M82 X-1 where X-2 is off and a clear view of the source is seen. A prominent hard excess is seen with NuSTAR when fitted with a pileup\*zwabs\*diskpbb model (top), which we account for with a simpl model (middle) and cutoffpl model (bottom).</figcaption></figure>

value from the 90% confidence interval for each pair of parameters. We calculate from how many of the 1000 simulations we find a correlation with a *p*-value less than 0.05. For all correlations involving parameters of the simpl model, we find that in less than 25% of the simulations, a *p*-value less than 0.05 is recovered. Therefore, we do not have confidence in these correlations being real, due to the large uncertainties in the parameters.

Furthermore, for the correlations where the uncertainties in the parameters are smaller, some of these correlations can be expected due to degeneracies between parameters. To determine if these correlations are driven by degeneracies, we explore the two-dimensional $\chi^2$ space around the best fit using the XSPEC command *steppar*, and overplot the $3\sigma$ contours on Figure 6. We do this only for one observation, *Chandra* obsID 18069, as it is computationally expensive, but this should be adequate to reveal any spectral degeneracies. This data set is a typical observation where X-2 is at low fluxes, pileup for X-1 is low, and the measured parameters are in the middle of the distributions.

The contours show that there is slight degeneracy between $L_X$ and $N_H$, but that does not appear to be strong enough to induce the correlation. There is no apparent degeneracy between $L_X$ and $\log_{10}$ norm is expected, and a clear degeneracy between $T_{in}$ and $\log_{10}$ norm is seen.

As the background is higher in the NuSTAR data than in the *Chandra* data, and the signal to noise lower, we investigated whether binning the NuSTAR spectra with more counts would

<figure><img src="image_6.png" /><figcaption>**Figure 5.** *Chandra* spectra of X-1 (red), X-2 (green) and the diffuse emission (purple), and the NuSTAR FPMA spectra of all sources (blue) unfolded through the spectral responses with the assumed spectral models.</figcaption></figure>

affect our results. We grouped the NuSTAR spectra with a minimum of 60 counts per bin rather than 20, and refit the spectra. We looked at the temperature of the diskpbb component, which is most sensitive to the NuSTAR data and the one of the key parameters involved in our results. We found the average difference in temperature between the stronger binning and the weaker one to be -0.02 keV, which is much smaller than the typical uncertainty on the parameter. We therefore conclude that the spectral binning does not affect our results.

Taking into account the uncertainties in parameters and degeneracies between them, we can only say with confidence that there is a correlation and therefore physical link between the neutral column density and the X-ray luminosity, and the inner disk temperature and X-ray luminosity.

We test the apparent dependence of $N_H$ on $L_X$ further by fixing $N_H$ at the approximate mean value of $1.3 \times 10^{22} \text{ cm}^{-2}$ in all fits and note the difference in $\chi^2$. In all cases, $\chi^2$ is worse or the same as when $N_H$ is a free parameter. For *Chandra* obsID 18062, $L_X$ is high and the $N_H$ measured is particularly high at $1.8 \times 10^{22} \text{ cm}^{-2}$. When $N_H$ is fixed at the mean value, the difference in $\chi^2$ is 80. For obsID 18068, $L_X$ is low and $N_H = 9 \times 10^{21} \text{ cm}^{-2}$. When fixed at $1.3 \times 10^{22} \text{ cm}^{-2}$, the $\chi^2$ increases by 150, supporting the finding that $N_H$ does indeed depend on $L_X$.
<table><caption>**Table 2**
Spectral Fitting Results for X-1</caption><thead><tr><th>Chandra<br>ObsID</th><th>N<sub>H</sub><br>(10<sup>22</sup> cm<sup>-2</sup>)</th><th>T<sub>in</sub><br>(keV)</th><th>p</th><th>log<sub>10</sub> norm</th><th>Γ</th><th>f<sub>scatt</sub></th><th>χ<sup>2</sup>/DoF</th><th>F<sub>X</sub><br>(10<sup>-11</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>L<sub>X</sub><br>(10<sup>40</sup> erg s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>17578</td><td>0.98<sup>+0.15</sup><sub>-0.33</sub></td><td>2.8<sup>+2.7</sup><sub>-1.9</sub></td><td>0.59<sup>+0.13</sup><sub>-0.03</sub></td><td>-2.6<sup>+2.2</sup><sub>-1.2</sub></td><td>3.0<sup>+2.0</sup><sub>-2.0</sub></td><td>0.98<sup>+0.02</sup><sub>-0.93</sub></td><td>941.5/840</td><td>0.97<sup>+0.04</sup><sub>-0.04</sub></td><td>1.3<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>17678</td><td>1.35<sup>+0.17</sup><sub>-0.08</sub></td><td>1.7<sup>+0.4</sup><sub>-0.4</sub></td><td>0.64<sup>+0.08</sup><sub>-0.03</sub></td><td>-0.8<sup>+0.7</sup><sub>-0.4</sub></td><td>3.7<sup>+0.7</sup><sub>-0.5</sub></td><td>0.73<sup>+0.27</sup><sub>-0.40</sub></td><td>1984.5/1259</td><td>4.26<sup>+0.19</sup><sub>-0.13</sub></td><td>5.6<sup>+0.3</sup><sub>-0.2</sub></td></tr><tr><td>18062</td><td>1.81<sup>+0.13</sup><sub>-0.13</sub></td><td>2.3<sup>+0.8</sup><sub>-0.7</sub></td><td>0.52<sup>+0.03</sup><sub>-0.02</sub></td><td>-1.8<sup>+0.6</sup><sub>-0.5</sub></td><td>3.3<sup>+0.4</sup><sub>-0.8</sub></td><td>0.60<sup>+0.40</sup><sub>-0.47</sub></td><td>2070.3/1246</td><td>3.73<sup>+0.08</sup><sub>-0.15</sub></td><td>4.9<sup>+0.1</sup><sub>-0.2</sub></td></tr><tr><td>18063</td><td>1.25<sup>+0.07</sup><sub>-0.11</sub></td><td>6.7<sup>+0.6</sup><sub>-3.1</sub></td><td>0.57<sup>+0.01</sup><sub>-0.01</sub></td><td>-3.7<sup>+0.7</sup><sub>-0.2</sub></td><td>4.5<sup>+0.5</sup><sub>-4.5</sub></td><td>0.00<sup>+0.26</sup><sub>-0.00</sub></td><td>1445.8/1282</td><td>1.82<sup>+0.05</sup><sub>-0.05</sub></td><td>2.4<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>18064</td><td>1.09<sup>+0.16</sup><sub>-0.15</sub></td><td>4.5<sup>+1.7</sup><sub>-2.3</sub></td><td>0.58<sup>+0.04</sup><sub>-0.03</sub></td><td>-3.3<sup>+0.9</sup><sub>-0.6</sub></td><td>2.7<sup>+2.4</sup><sub>-2.7</sub></td><td>0.23<sup>+0.77</sup><sub>-0.23</sub></td><td>1486.0/1066</td><td>1.04<sup>+0.03</sup><sub>-0.03</sub></td><td>1.4<sup>+0.0</sup><sub>-0.0</sub></td></tr><tr><td>18065</td><td>1.25<sup>+0.08</sup><sub>-0.08</sub></td><td>2.3<sup>+0.7</sup><sub>-0.3</sub></td><td>0.57<sup>+0.02</sup><sub>-0.02</sub></td><td>-1.6<sup>+0.2</sup><sub>-0.4</sub></td><td>3.7<sup>+1.3</sup><sub>-0.6</sub></td><td>0.92<sup>+0.08</sup><sub>-0.60</sub></td><td>2069.2/1404</td><td>4.21<sup>+0.16</sup><sub>-0.10</sub></td><td>5.5<sup>+0.2</sup><sub>-0.1</sub></td></tr><tr><td>18067</td><td>1.38<sup>+0.11</sup><sub>-0.10</sub></td><td>2.0<sup>+0.5</sup><sub>-0.6</sub></td><td>0.56<sup>+0.03</sup><sub>-0.02</sub></td><td>-1.4<sup>+0.5</sup><sub>-0.2</sub></td><td>3.5<sup>+1.1</sup><sub>-0.7</sub></td><td>0.58<sup>+0.42</sup><sub>-0.36</sub></td><td>2410.8/1567</td><td>3.10<sup>+0.12</sup><sub>-0.11</sub></td><td>4.0<sup>+0.2</sup><sub>-0.1</sub></td></tr><tr><td>18068</td><td>0.94<sup>+0.12</sup><sub>-0.14</sub></td><td>3.8<sup>+0.9</sup><sub>-1.7</sub></td><td>0.64<sup>+0.06</sup><sub>-0.02</sub></td><td>-2.7<sup>+1.0</sup><sub>-0.4</sub></td><td>2.4<sup>+1.1</sup><sub>-1.4</sub></td><td>0.31<sup>+0.58</sup><sub>-0.24</sub></td><td>1502.7/1366</td><td>0.94<sup>+0.07</sup><sub>-0.05</sub></td><td>1.2<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>18069</td><td>1.25<sup>+0.10</sup><sub>-0.09</sub></td><td>3.6<sup>+1.3</sup><sub>-1.0</sub></td><td>0.56<sup>+0.02</sup><sub>-0.02</sub></td><td>-2.7<sup>+0.5</sup><sub>-0.3</sub></td><td>2.7<sup>+1.1</sup><sub>-1.4</sub></td><td>0.28<sup>+0.33</sup><sub>-0.26</sub></td><td>1529.0/1275</td><td>2.31<sup>+0.09</sup><sub>-0.06</sub></td><td>3.0<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>18070</td><td>1.30<sup>+0.10</sup><sub>-0.11</sub></td><td>3.1<sup>+0.6</sup><sub>-0.5</sub></td><td>0.56<sup>+0.02</sup><sub>-0.02</sub></td><td>-2.5<sup>+0.3</sup><sub>-0.1</sub></td><td>4.0<sup>+1.0</sup><sub>-1.4</sub></td><td>0.91<sup>+0.09</sup><sub>-0.86</sub></td><td>1598.2/1252</td><td>2.02<sup>+0.05</sup><sub>-0.05</sub></td><td>2.6<sup>+0.1</sup><sub>-0.1</sub></td></tr></tbody></table>

**Note.** Best-fit parameters for the *diskpbb* model fitted to X-1. Fluxes and luminosities are given in the 0.5–30 keV range, and are corrected for absorption and pileup. Luminosities are calculated assuming a distance of 3.3 Mpc to M82.

<table><caption>**Table 3**
Spectral Fitting Results for X-2</caption><thead><tr><th>Chandra<br>ObsID</th><th>N<sub>H</sub><br>(10<sup>22</sup> cm<sup>-2</sup>)</th><th>Γ</th><th>E<sub>C</sub><br>(keV)</th><th>log<sub>10</sub> norm</th><th>χ<sup>2</sup>/DoF</th><th>F<sub>X</sub><br>(10<sup>-11</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>L<sub>X</sub><br>(10<sup>40</sup> erg s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>17578</td><td>1.7<sup>+2.4</sup><sub>-1.4</sub></td><td>-1.3<sup>+3.4</sup><sub>-3.6</sub></td><td>1.5<sup>+22</sup><sub>-0.4</sub></td><td>-3.6<sup>+0.9</sup><sub>-0.5</sub></td><td>899.9/840</td><td>0.4<sup>+0.3</sup><sub>-0.1</sub></td><td>0.5<sup>+0.4</sup><sub>-0.1</sub></td></tr><tr><td>17678</td><td>3.8<sup>+1.2</sup><sub>-1.1</sub></td><td>1.7<sup>+0.7</sup><sub>-0.7</sub></td><td>15<sup>+77</sup><sub>-16</sub></td><td>-2.7<sup>+0.6</sup><sub>-0.4</sub></td><td>1378.4/1259</td><td>1.4<sup>+0.5</sup><sub>-0.3</sub></td><td>1.8<sup>+0.6</sup><sub>-0.4</sub></td></tr><tr><td>18063</td><td>1.9<sup>+0.4</sup><sub>-0.4</sub></td><td>-1.2<sup>+0.4</sup><sub>-0.4</sub></td><td>2.0<sup>+0.4</sup><sub>-0.3</sub></td><td>-3.6<sup>+0.2</sup><sub>-0.1</sub></td><td>1370.8/1282</td><td>1.0<sup>+0.1</sup><sub>-0.1</sub></td><td>1.2<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>18064</td><td>2.3<sup>+1.3</sup><sub>-1.0</sub></td><td>-1.3<sup>+1.3</sup><sub>-1.6</sub></td><td>1.4<sup>+0.9</sup><sub>-0.4</sub></td><td>-3.4<sup>+0.4</sup><sub>-0.4</sub></td><td>1233.9/1066</td><td>0.5<sup>+0.2</sup><sub>-0.1</sub></td><td>0.6<sup>+0.3</sup><sub>-0.1</sub></td></tr><tr><td>18067</td><td>3.0<sup>+0.7</sup><sub>-0.6</sub></td><td>1.3<sup>+0.5</sup><sub>-0.6</sub></td><td>16<sup>+25</sup><sub>-10</sub></td><td>-3.0<sup>+0.3</sup><sub>-0.2</sub></td><td>1725.4/1567</td><td>1.0<sup>+0.1</sup><sub>-0.2</sub></td><td>1.3<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>18068</td><td>0.8<sup>+0.2</sup><sub>-0.8</sub></td><td>-2.9<sup>+6.9</sup><sub>-2.0</sub></td><td>1.2<sup>+0.0</sup><sub>-0.1</sub></td><td>-4.0<sup>+0.0</sup><sub>-0.0</sub></td><td>1478.3/1366</td><td>0.7<sup>+0.0</sup><sub>-0.0</sub></td><td>1.0<sup>+0.1</sup><sub>-0.0</sub></td></tr></tbody></table>

**Note.** Best-fit parameters for the *cutoffp1* model fitted to X-2. Fluxes and luminosities are given in the 0.5–30 keV range, and are corrected for absorption. Luminosities are calculated assuming a distance of 3.3 Mpc to M82.

<table><caption>**Table 4**
Spearman's Rank Correlation Results for X-1</caption><thead><tr><th></th><th>L<sub>X</sub></th><th>N<sub>H</sub></th><th>T<sub>in</sub></th><th>p</th><th>log<sub>10</sub>norm</th><th>Γ</th></tr></thead><tbody><tr><td>N<sub>H</sub></td><td>0.76</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>0.0111</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>T<sub>in</sub></td><td>-0.76</td><td>-0.56</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>0.011</td><td>0.090</td><td></td><td></td><td></td><td></td></tr><tr><td>p</td><td>-0.33</td><td>-0.62</td><td>0.03</td><td></td><td></td><td></td></tr><tr><td></td><td>0.347</td><td>0.054</td><td>0.934</td><td></td><td></td><td></td></tr><tr><td>log<sub>10</sub>norm</td><td>0.79</td><td>0.62</td><td>-0.99</td><td>-0.10</td><td></td><td></td></tr><tr><td></td><td>0.006</td><td>0.054</td><td>0.000</td><td>0.777</td><td></td><td></td></tr><tr><td>Γ</td><td>0.48</td><td>0.58</td><td>-0.22</td><td>-0.36</td><td>0.28</td><td></td></tr><tr><td></td><td>0.162</td><td>0.082</td><td>0.533</td><td>0.310</td><td>0.425</td><td></td></tr><tr><td>f<sub>scatt</sub></td><td>0.28</td><td>0.09</td><td>-0.67</td><td>0.05</td><td>0.65</td><td>0.21</td></tr><tr><td></td><td>0.425</td><td>0.803</td><td>0.033</td><td>0.881</td><td>0.043</td><td>0.043</td></tr></tbody></table>

**Note.** For each pair of parameters we list the rank correlation coefficient (top) and the two-sided significance of its deviation from zero (bottom).

The anticorrelation between $L_X$ and $T_{in}$ that we find here is in contrast with the correlation found between these two parameters for M82 X-1 by Feng & Kaaret (2010), where the apparent $L_X \propto T^4$ relationship lead the authors to conclude that the source was observed in the thermal dominant state. To determine the exponent on the relationship, we take the logarithm of both $L_X$ and $T_{in}$, such that $\log_{10}L_X = \alpha\log_{10}T_{in}+\beta$ and conduct a linear regression analysis. We use the IDL tool *linmix_err.pro*, which takes into account uncertainties in both parameters. We find that $\log_{10}L_X$

= (-1.47 ± 1.00)log<sub>10</sub>T<sub>in</sub> + (41.16 ± 0.47), where the uncertainties are 1σ. Therefore, the exponent is -1.47 ± 1.00, which is a >5σ difference from the results from F10, where the exponent was implied to be 4, although the uncertainty in this parameter was not listed. If we run the same linear regression analysis on the results presented in F10, we find that the relationship is essentially unconstrained.

The key differences between our analysis and that of F10 are that they did not have accompanying *NuSTAR* data, so their analysis was restricted to the narrow bandpass of 0.7–7 keV. The spectral models used are also different, whereby F10 used the *diskpbb* model, which is related to *diskpbb* that we use when *p* is fixed at 0.75. However, as we showed in B16, and confirm here, the data are inconsistent with the *p* = 0.75 that describes a geometrically thin accretion disk. Furthermore, without *NuSTAR* data, the high-energy excess that we detect and fit with *simpl* was not observable by F10. Finally, the uncertainties related to degeneracies in the *pileup* model are reduced when data from an instrument without pileup such as *NuSTAR* is used.

We investigate what effect pileup has on our results by exploring the dependence of $L_X$, $N_H$ and $T_{in}$ on the *pileup* model parameter, $\alpha$. We find that none of these parameters show any dependence on $\alpha$, leading us to conclude that this model component does not drive our results. Furthermore, as mentioned earlier, for three observations, odsID 17678, 18065, and 18067, the pileup fraction for X-1 is greater than 10%, and therefore the *pileup* model used in the spectral fitting may not be able to reliably account for this effect. We test the
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Relationship between the spectral parameters of the *zwabs*\*simpl\*diskpbb model for X-1 where a significant correlation was found. Red data points indicate observations where X-2 is at low fluxes and thus the view of the emission from X-1 is unambiguous. We also overplot the $3\sigma$ $\chi^2$ contour in blue, shifted for clarity, from *Chandra* obsID 18069 to demonstrate if the correlation is driven by degeneracy between the parameters, which appears to be the case for $T_{in}$ and $\log_{10}\text{norm}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Relationship between the spectral parameters of the *zwabs*\*cutoffpl model for X-2 where evidence for a correlation has been found. Red data points indicate the best-fit parameters of the pulsed emission. We also overplot the $3\sigma$ $\chi^2$ contour in blue, shifted for clarity, from *Chandra* obsID 18063 to demonstrate if the correlation is driven by degeneracy between the parameters, which appears to be the case for all parameter pairs.</figcaption></figure>

dependence of our results on these observations by removing them from our analysis. In doing so, we still find a significant correlation between $L_X$ and $N_H$, and indeed the significance increases, but for $L_X$ and $T_{in}$ the correlation is no longer significant. A fit with a linear relationship reveals $\log_{10}L_X = (-1.23 \pm 2.29)\log_{10}T_{in}+(41.03 \pm 1.28)$ and therefore we can only rule out a $L_X \propto T^4$ dependency at $\sim 2\sigma$.

Finally, we note that some of our spectral fits have large $\chi^2$ relative to the number of degrees of freedom (DoF), indicating an unacceptable fit. This is likely due to the fits being quite complex, with many DoFs and data sets being fitted simultaneously. However, we note the presence of a potential excess of counts in the *NuSTAR* data at 3–4 keV that has been found in bright X-ray binaries and is a known calibration issue. We find that if we ignore data from *NuSTAR* below 5 keV, the fits improve and the result of this is to systematically reduce the temperature of the *diskpbb* component. This reduction is consistent with the uncertainties on this component, however, and does not alter our result that there is an anticorrelation between $L_X$ and $T_{in}$ because the effect is systematic across all observations.

These findings therefore rule out a thermal state for sub-Eddington accretion and therefore do not support M82 X-1 as an IMBH candidate.

## 5. New Insights into M82 X-2

We conduct the same analysis for X-1 on X-2, where the spectral parameters are plotted against each other in Figure 7.

<table><caption>**Table 5**
Spearman's Rank Correlation Results for X-2</caption><thead><tr><th></th><th>L<sub>X</sub></th><th>N<sub>H</sub></th><th>E<sub>C</sub></th><th>Γ</th></tr></thead><tbody><tr><td>N<sub>H</sub></td><td>0.71</td><td></td><td></td><td></td></tr><tr><td></td><td>0.111</td><td></td><td></td><td></td></tr><tr><td>E<sub>C</sub></td><td>0.71</td><td>0.77</td><td></td><td></td></tr><tr><td></td><td>0.111</td><td>0.072</td><td></td><td></td></tr><tr><td>Γ</td><td>0.77</td><td>0.83</td><td>0.94</td><td></td></tr><tr><td></td><td>0.072</td><td>0.042</td><td>0.005</td><td></td></tr><tr><td>log<sub>10</sub>norm</td><td>0.71</td><td>1.00</td><td>0.77</td><td>0.83</td></tr><tr><td></td><td>0.111</td><td>0.000</td><td>0.072</td><td>0.042</td></tr></tbody></table>

**Note.** For each pair of parameters, we list the rank correlation coefficient (top) and the two-sided significance of its deviation from zero (bottom).

and the correlation analysis is shown in Table 5. Here our analysis suggests correlations between $N_H$ and $\Gamma$, $E_C$ and $\Gamma$, $N_H$ and $\log_{10}\text{norm}$, and $\Gamma$ and $\log_{10}\text{norm}$. However, the $\chi^2$ contours show that it is likely that strong degeneracies between these parameters drive these apparent correlations (Figure 7). We therefore do not find evidence for any significant spectral evolution in X-2.

While the emission from X-2 is difficult to disentangle from the other sources of emission in M82, we were able to isolate the pulsed emission in the *NuSTAR* band from this source in Brightman et al. (2016a). We found that the pulsed emission is best fit by a power law with a high-energy cutoff, where $\Gamma = 0.6 \pm 0.3$ and $E_C = 14_{-3}^{+5}$ keV. In Figure 7, we show the
parameters of the pulsed emission as a separate data point for comparison. We see that the values for $\Gamma$ and $E_C$ from our broadband fits are consistent with the pulsed emission when X-2 is at its highest luminosities, $L_X > 10^{40}$ erg s$^{-1}$, indicating that at these times the pulsations are most likely to be detected. Bachetti et al. (2019) have recently detected pulsations again from a NuSTAR observation taken on 2016 September 10, which unfortunately did not have any simultaneous Chandra observations, so we did not include it in our analysis here.

# 6. Discussion and Implications

## 6.1. M82 X-1 as an Intermediate-mass Black Hole Candidate

M82 X-1 has been claimed to be an intermediate-mass black hole candidate based on its high X-ray luminosity, twin QPOs, and $L_X \propto T^4$ scaling, all of which put the mass of the black hole at $\sim 10^2 M_\odot$. Here, we find that we can rule out a $L_X \propto T^4$ scaling. This scaling relationship was expected from a standard accretion disk which exists at moderate accretion rates, and allows an estimate of the black hole mass from measurements of the inner edge of the accretion disk. Without these pieces of evidence, the status of M82 X-1 as an IMBH accretor is less certain.

## 6.2. M82 X-1 as a Super-Eddington Accreting Stellar Remnant

The X-ray properties of M82 X-1 may be explained by it harboring a super-Eddington accreting stellar-remnant black hole or neutron star. The spectral shape of M82 X-1, consisting of a broadened disk with a high-energy tail is very similar to all other ULXs with high-quality broadband spectral data from NuSTAR (e.g., Walton et al. 2018b). This sample includes the known neutron star accretors NGC 5907 ULX1 and NGC 7793 P13, and at first glance their spectral shapes are not dissimilar from the rest of the sample. This implies that these ULXs, including M82 X-1, are also super-Eddington accretors, although it is still not known whether they are powered by neutron stars or black holes.

One popular model to explain the spectral evolution of ULXs is that they are stellar-remnant black holes accreting at super-Eddington rates. In this model, a powerful wind is radiatively driven from the accretion disk (e.g., Poutanen et al. 2007), as recently revealed through the detection of highly ionized material in the high-resolution X-ray spectra of NGC 1313 X-1 (Pinto et al. 2016) among others (Pinto et al. 2017; Kosec et al. 2018) and the detection of blueshifted iron-K absorption (Walton et al. 2016). Regarding the link between $N_H$ and $L_X$ as seen in M82 X-1, if this source is a stellar-mass black hole accreting at super-Eddington rates, as the mass accretion rate increases an increase in the X-ray luminosity follows, which drives further outflow of material from the system. Depending on the line of sight, this can cause an increase in the line of sight absorption. However, this material is expected to be highly ionized at small radii. Middleton et al. (2015) also found for NGC 1313 X-1 that the neutral column density anticorrelates with spectral hardness, suggesting that at large radii, there is a cool, neutral component of the outflow, where the column density is linked to mass loss via increase mass accretion rate, as predicted by Poutanen et al. (2007).

A negative relationship between $L_X$ and $T_{in}$ has been observed in other ULXs. For example, Luangtip et al. (2016) studied the spectral evolution of Holmberg IX X-1, finding that the peak of the spectrum decreases with luminosity, suggesting

an anticorrelation between $L_X$ and $T$ (see also Walton et al. 2017). Kajava & Poutanen (2009) also explored the relationship between $L_X$ and $T$ in a sample of ULXs with XMM-Newton and Chandra observations. While they found that for several sources fitted with a multicolor disk model follow the $L_X \propto T^4$ scaling, sources with a disk plus power-law spectral shape show a negative $L_X \propto T^{-3.5}$ scaling, similar to what we find here and predicted by Poutanen et al. (2007).

# 7. Summary and Conclusions

We conducted a comprehensive investigation into the spectral evolution of the ultraluminous X-ray sources M82 X-1 and X-2 using 10 simultaneous Chandra and NuSTAR observations. The Chandra data allowed us to spatially resolve the sources, separated by only 5", while the NuSTAR data allowed a broadband X-ray spectral analysis. We found that for X-1, the luminosity of the disk scales with the inner temperature as $L \propto T^{-3/2}$, which is contrary to previous findings of a $L \propto T^4$ scaling that supported a standard accretion disk powering the system. We furthermore find evidence that the neutral column density of the material in the line of sight increases with $L_X$, perhaps due to an increased mass outflow with accretion rate. For X-2, we do not find any significant spectral evolution, but we can constrain the spectral parameters, and find the broadband emission to be consistent with the pulsed emission at the highest X-ray luminosities.

This research has made use of data obtained with NuSTAR, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by NASA. We thank the NuSTAR Operations, Software and Calibration teams for support with the execution and analysis of these observations. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA). Support for this work was provided by the National Aeronautics and Space Administration through Chandra Award Number GO6-17080X issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. This research has made use of software provided by the Chandra X-ray Center (CXC) in the application package CIAO. We also acknowledge the use of public data from the Swift data archive. D.J.W. acknowledges support from an STFC Ernest Rutherford Fellowship.

Facilities: Chandra (ACIS), NuSTAR, Swift (XRT).

# ORCID iDs

Dominic J. Walton https://orcid.org/0000-0001-5819-3552

Yanjun Xu https://orcid.org/0000-0003-2443-3698

Hannah P. Earnshaw https://orcid.org/0000-0001-5857-5622

Daniel Stern https://orcid.org/0000-0003-2686-9241

Didier Barret https://orcid.org/0000-0002-0393-9190

# References

Arnaud, K. A. 1996, in ASP Conf. Ser. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 17

Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, Natur, 514, 202

Bachetti, M., Maccarone, T. J., Brightman, M., et al. 2019, arXiv:1905.06423
Brightman, M., Harrison, F., Walton, D. J., et al. 2016a, ApJ, 816, 60

Brightman, M., Harrison, F. A., Bachetti, M., et al. 2019, ApJ, 873, 115

Brightman, M., Harrison, F. A., Barret, D., et al. 2016b, ApJ, 829, 28

Cash, W. 1979, ApJ, 228, 939

Colbert, E. J. M., & Mushotzky, R. F. 1999, ApJ, 519, 89

Davis, J. E. 2001, ApJ, 562, 575

Dewangan, G. C., Titarchuk, L., & Griffiths, R. E. 2006, ApJL, 637, L21

Earnshaw, H. P., Roberts, T. P., Middleton, M. J., Walton, D. J., & Mateos, S. 2019, MNRAS, 483, 5554

Fabbiano, G. 1989, ARA&A, 27, 87

Farrell, S. A., Webb, N. A., Barret, D., Godet, O., & Rodrigues, J. M. 2009, Natur, 460, 73

Feng, H., & Kaaret, P. 2010, ApJL, 712, L169

Foley, R. J., Fox, O. D., McCully, C., et al. 2014, MNRAS, 443, 2887

Fürst, F., Walton, D. J., Stern, D., et al. 2017, ApJ, 834, 77

Giacconi, R., Branduardi, G., Briel, U., et al. 1979, ApJ, 230, 540

Harrison, F. A., Craig, W. W., Christensen, F. E., Hailey, C. J., & Zhang, W. W. 2013, ApJ, 770, 103

Israel, G. L., Belfiore, A., Stella, L., et al. 2017, Sci, 355, 817

Kaaret, P., Simet, M. G., & Lang, C. C. 2006, ApJ, 646, 174

Kajava, J. J. E., & Poutanen, J. 2009, MNRAS, 398, 1450

Kosec, P., Pinto, C., Walton, D. J., et al. 2018, MNRAS, 479, 3978

Luangtip, W., Roberts, T. P., & Done, C. 2016, MNRAS, 460, 4417

Madsen, K. K., Harrison, F. A., Markwardt, C. B., et al. 2015, ApJS, 220, 8

Matsumoto, H., Tsuru, T. G., Koyama, K., et al. 2001, ApJL, 547, L25

Middleton, M. J., Walton, D. J., Fabian, A., et al. 2015, MNRAS, 454, 3134

Miller, J. M., Fabbiano, G., Miller, M. C., & Fabian, A. C. 2003, ApJL, 585, L37

Mizuno, T., Ohnishi, T., Kubota, A., Makishima, K., & Tashiro, M. 1999, PASJ, 51, 663

Mucciarelli, P., Casella, P., Belloni, T., Zampieri, L., & Ranalli, P. 2006, MNRAS, 365, 1123

Pasham, D. R., Strohmayer, T. E., & Mushotzky, R. F. 2014, Natur, 513, 74

Pinto, C., Alston, W., Soria, R., et al. 2017, MNRAS, 468, 2865

Pinto, C., Middleton, M. J., & Fabian, A. C. 2016, Natur, 533, 64

Poutanen, J., Lipunova, G., Fabrika, S., Butkevich, A. G., & Abolmasov, P. 2007, MNRAS, 377, 1187

Ptak, A., & Griffiths, R. 1999, ApJL, 517, L85

Rephaeli, Y., & Gruber, D. 2002, A&A, 389, 752

Steiner, J. F., Narayan, R., McClintock, J. E., & Ebisawa, K. 2009, PASP, 121, 1279

Strohmayer, T. E., & Mushotzky, R. F. 2003, ApJL, 586, L61

Sutton, A. D., Roberts, T. P., Gladstone, J. C., & Walton, D. J. 2015, MNRAS, 450, 787

Sutton, A. D., Roberts, T. P., Gladstone, J. C., et al. 2013, MNRAS, 434, 1702

Walton, D. J., Fürst, F., Harrison, F. A., et al. 2017, ApJ, 839, 105

Walton, D. J., Fürst, F., Harrison, F. A., et al. 2018a, MNRAS, 473, 4360

Walton, D. J., Fürst, F., Heida, M., et al. 2018b, ApJ, 856, 128

Walton, D. J., Middleton, M. J., Pinto, C., et al. 2016, ApJL, 826, L26

Walton, D. J., Middleton, M. J., Rana, V., et al. 2015, ApJ, 806, 65

Walton, D. J., Roberts, T. P., Mateos, S., & Heard, V. 2011, MNRAS, 416, 1844

Watson, M. G., Stanger, V., & Griffiths, R. E. 1984, ApJ, 286, 144

Weisskopf, M. C. 1999, arXiv:astro-ph/9912097