Eclipsing ULXs in M 51

Two eclipsing ultraluminous X-ray sources in M 51

R. Urquhart11affiliationmark: and R. Soria11affiliationmark: 22affiliationmark: International Centre for Radio Astronomy Research, Curtin University, GPO Box U1987, Perth, WA 6845, Australia
Sydney Institute for Astronomy, School of Physics A28, The University of Sydney, Sydney, NSW 2006, Australia
Abstract

We present the discovery, from archival Chandra and XMM-Newton data, of X-ray eclipses in two ultraluminous X-ray sources (ULXs), located in the same region of the galaxy M 51: CXOM51 J132940.0471237 (ULX-1, for simplicity) and CXOM51 J132939.5471244 (ULX-2). Three eclipses were detected for ULX-1, two for ULX-2. The presence of eclipses puts strong constraints on the viewing angle, suggesting that both ULXs are seen almost edge-on and are certainly not beamed towards us. Despite the similar viewing angles and luminosities ( erg s in the 0.3–8 keV band for both sources), their X-ray properties are different. ULX-1 has a soft spectrum, well fitted by Comptonization emission from a medium with electron temperature keV. ULX-2 is harder, well fitted by a slim disk with –1.8 keV and normalization consistent with a 10  black hole. ULX-1 has a significant contribution from multi-temperature thermal plasma emission ( erg s); about 10% of this emission remains visible during the eclipses, proving that the emitting gas comes from a region slightly more extended than the size of the donor star. From the sequence and duration of the Chandra observations in and out of eclipse, we constrain the binary period of ULX-1 to be either 6.3 days, or 12.5–13 days. If the donor star fills its Roche lobe (a plausible assumption for ULXs), both cases require an evolved donor; most likely a blue supergiant, given the young age of the stellar population in that galactic environment.

accretion, accretion disks – stars: black holes – X-rays: binaries

1. Introduction

Ultraluminous X-ray sources (ULXs) are the high-luminosity end of the X-ray binary population, with X-ray luminosities erg s, which is the approximate peak luminosity of Galactic stellar-mass black holes (BHs). The most likely explanation for the vast majority of ULXs is that they are stellar-mass BHs (or neutron stars; Bachetti et al., 2014) accreting well above the critical accretion rate and their luminosity is a few times the classical Eddington limit of spherical accretion. Another possibility is that ULXs are powered by accreting BHs up to (Belczynski et al., 2010), several times more massive than typical Galactic stellar-mass BHs (: Kreidberg et al. 2012). In addition, ULXs may appear more luminous because their X-ray emission is partly collimated along our line of sight. This may happen at super-critical accretion rates, because of the predicted formation of a dense radiatively-driven disk outflow and a lower-density polar funnel, along which more photons can escape (Ohsuga & Mineshige, 2011; Jiang et al., 2014; Sa̧dowski & Narayan, 2015). Finally, some of the brightest ULXs may contain a population of intermediate-mass BHs (; Farrell et al., 2009; Zolotukhin et al., 2016). It is difficult to determine the relative contribution of those three factors (mass, accretion rate and viewing angle), and therefore also determine the true isotropic luminosity and accretion rate of ULXs, without at least a direct constraint on their viewing angle.

There is already indirect evidence that ULXs are not strongly beamed. Modelling of the optical light curve from the irradiated donor star in NGC 7793-P13 showed (Motch et al., 2014) that the source is viewed at an angle 20 and more likely much higher; thus, in that case, super-Eddington accretion is the reason for the high luminosity, not a heavier BH or a down-the-funnel view. Studies of large (100 pc) photo-ionized and/or shock-ionized plasma bubbles around ULXs (Pakull & Mirioni, 2002; Pakull & Grisé, 2008) provide other clues about viewing angles: if fast-accreting BHs appeared as ULXs only for a narrow range of face-on inclinations, we would see many more of those large ionized bubbles without a ULX inside; moreover, the true X-ray luminosity of a beamed source (much lower than the apparent luminosity) would not be enough to explain the strong He II emission observed from some of the photo-ionized ULX bubbles (Pakull et al., 2006). Both the fact that most ULX bubbles do contain a bright, central X-ray source, and the fact that (in photo-ionized bubbles) the apparent X-ray photon flux from the central source is consistent with the He II photon flux from the bubble, suggest that, statistically, ULXs are seen over a broad range of viewing angles. X-ray spectroscopic studies can also be used to qualitatively constrain ULX viewing angles: it was suggested (Sutton et al., 2013) that ULXs seen at lower inclination (down the polar funnel) have harder X-ray spectra while those seen at higher inclination (through the disk wind) have softer spectra with a lower-energy downturn, due to a higher degree of Compton scattering in the wind. This interpretation is consistent with the presence of absorption and emission features (interpreted as signatures of the outflow) in the X-ray spectra of ULXs with softer spectra (Middleton et al., 2014, 2015b). It is also in agreement with a higher degree of short-term variability (interpreted as the imprint of a clumpy wind) in sources with softer spectra (Middleton et al., 2015a).

Apart from those indirect or statistical arguments, until recently there was no bright extragalactic stellar-mass BH for which the viewing angle could be directly pinned down. We have now discovered two such sources, both located in the same spiral arm of the spiral galaxy M 51; in fact, surprisingly, they appear projected in the sky within only of each other (see Figure 1). Both sources have X-ray luminosities 10 erg s, and crucially, they both show sharp X-ray drops and rebrightenings, which we interpret as eclipses by their donor stars, occulting the inner region of the disk. The presence of eclipses places a lower limit on the inclination angle () as we must be viewing the X-ray sources near edge-on. In this paper, we present the eclipse discovery and the main X-ray timing and spectral properties of the two sources. We will also briefly discuss more general implications and opportunities provided by the detection of eclipses, for our modelling of these systems. In a companion paper (Soria et al., in prep.) we will present a study of the optical counterparts and other interesting, newly-discovered properties of those same two ULXs, which show optical and radio evidence of jets and outflows.

2. Targets of our study

M 51, also known as the Whirlpool Galaxy, is an interacting face-on spiral at a distance of Mpc (Bose & Kumar, 2014). The two eclipsing sources discussed in this paper are those catalogued as CXOM51 J132940.0471237 (henceforth, ULX-1) and CXOM51 J132939.5471244 (henceforth, ULX-2) in Terashima & Wilson (2004). We re-estimated their positions using all the Chandra data available to-date, and obtained RA (J2000) , Dec. (J2000) for ULX-1, and RA (J2000) , Dec. (J2000) for ULX-2. Both positions are subject to the standard uncertainty in the absolute astrometry of Chandra pointings, 0”.6 at the 90% confidence level. A more precise determination of their positions is left to a follow-up study (Soria et al., in prep.) of their optical and radio counterparts.

ULX-1 and ULX-2 were first discovered as a single unresolved source by the Einstein Observatory (Palumbo et al., 1985). This was followed up with observations with ROSAT (source C in Ehle et al. 1995; source R7 in Marston et al. 1995). The higher spatial resolution of Chandra’s Advanced CCD Imaging Spectrometer (ACIS) finally led to the two sources being resolved (source 6 and source 5 in Terashima & Wilson, 2004). ULX-1 was found to be a relatively soft source, with very few counts above . ULX-2 was found to be variable, decreasing in luminosity by a factor of between observations (Terashima & Wilson, 2004). Further spectral studies of the two sources, based on a 2003 XMM-Newton observation, were carried out by Dewangan et al. (2005). With a much larger database of Chandra and XMM-Newton observation available since then, we have now studied the two sources in more detail, and found more intriguing properties.

3. Data Analysis

M 51 was observed by Chandra/ACIS-S fourteen times between 2000 and 2012: two of those observations were too short (2 ks) to be useful, and another two did not include our sources in the field of view; the other 10 observations are listed in Table 1. (See Kuntz et al. 2016 for a full catalog and discussion of all the Chandra sources in M 51.) We downloaded the Chandra data from the public archives and re-processed them using standard tasks within the Chandra Interactive Analysis of Observations (CIAO) Version 4.7 software package (Fruscione et al., 2006). Any intervals with high particle backgrounds were filtered out. We extracted spectra and light-curves for ULX-1 and ULX-2 using circular regions of 4″ radii and local background regions three times as large as the source regions. For each observation, background-subtracted light curves were created with the CIAO task dmextract. Spectra were extracted with specextract, and were then grouped to a minimum of 15 counts per bin, for fitting.

M 51 was also observed by XMM-Newton nine times between 2003 and 2011, although no data were recorded on three occasions due to strong background flaring (Table 1). We downloaded the XMM-Newton data from NASA’s High Energy Astrophysics Science Archive Research Centre (HEASARC) archive. We used the European Photon Imaging Camera (EPIC) observations and re-processed them using standard tasks in the Science Analysis System (SAS) version 14.0.0 software package; we filtered out high particle background exposure intervals. Due to the lower spatial resolution of XMM-Newton/EPIC, the ULXs cannot be entirely visually resolved, although the elongated appearance of the EPIC source is consistent with the two separate Chandra sources (as discussed in Section 4.1). For each observation we extracted a single background-subtracted light-curve and spectrum for both sources combined, using a circular extraction region of 20″ radius, and a local background region that is at least three times larger, does not fall onto any chip gap and is of similar distance to the readout nodes as the source region. Standard flagging criteria #XMMEA_EP and #XMMEA_EM were used for pn and MOS respectively, along with FLAG=0. We also selected patterns 0–4 for pn and 0–12 for MOS. For our timing study, we extracted light-curves with the SAS tasks evselect and epiclccorr. For our spectral study, we extracted individual pn, MOS1 and MOS2 spectra with standard xmmselect tasks; whenever possible, we combined the pn, MOS1 and MOS2 spectra of each observation with epicspeccombine, to create a weighted-average EPIC spectrum. In some observations, the pn data were not usable because the source falls onto a chip gap; in those cases, we used only the MOS1 and MOS2 data in epicspeccombine. Finally, we grouped the spectra to a minimum of 20 counts per bin so that we could use Gaussian statistics.

For both Chandra and XMM-Newton data, spectral fitting was performed with XSPEC version 12.8.2 (Arnaud, 1996). Timing analysis was conducted with standard FTOOLS tasks (Blackburn, 1995), such as lcurve, efsearch and statistics. Imaging analysis was done with HEASARC’s DS9 visualization package, and adaptive image smoothing with CIAO’s csmooth routine.

4. Results

4.1. Eclipses

4.1.1 ULX-1 eclipses and dips in the Chandra data

From our inspection of the Chandra light-curves, we have discovered 3 epochs (ObsIDs 1622, 13813 and 13814) in which the flux of ULX-1 is strongly reduced for at least part of the observation (Table 2 and Figures 2, 3, 4). The transition between the long-term-average flux level and the lower level occurs too quickly ( s) to be explained by a state transition in the inflow, or a change in the mass accretion rate. Our identification of the low state in ObsID 1622 as a true stellar eclipse rather than a dip may be debatable, given that the flux drop happens right at the start of the observation; however, the presence of eclipses is very clear in ObsIDs 13813 (2012 September 9) and 13814 (2012 September 20), which show a low-to-high and a high-to-low transition, respectively. We also checked that ULX-1 is not at the edge of the chip, there are no instrumental glitches, and no other source in the field has a count-rate step change at the same time. We conclude that the simplest and most logical explanation is an eclipse of the X-ray emitting region by the donor star. The flux during the eclipse is not exactly zero: by stacking the time intervals during eclipses, we can find a faint but statistically significant residual emission, softer than the emission outside eclipses. We will discuss the spectrum of the residual emission in Section 4.4.

The way ULX-1 enters the eclipse in ObsID 13814 (Figure 4) is also interesting. The transition to eclipse in the soft band (0.3–1.2 keV) appears less sharp than the transition in the hard band (1.2–7.0 keV): the soft-band count rate drops to effectively zero in 4 ks, while the same transition happens in 1 ks for the hard band. This can be explained if the softer X-ray photons come from a more extended region that takes longer to be completely occulted than the effectively point-like central region responsible for the harder X-ray photons (Bonnet-Bidaud et al., 2001); for example, the softer emission may have contributions from the outer, cooler parts of an outflow. However, we cannot rule out that the discrepancy is simply due to small-number statistics.

Finally, we find a deep dip in the Chandra light curve of ULX-1 during ObsID 13812 (Figure 5). The count rate drops to zero and then recovers to the pre-dip level, just like during an eclipse. However, the short duration (20 ks) and double-dipping substructure of this phase suggest that this occultation is not due to the companion star; we suggest that it is more likely the result of lumps or other inhomogeneities in the thick outer rim of the disk, or is caused by the accretion stream overshooting the point of impact in the outer disk and covering our view of the inner regions (Frank et al., 1987; Armitage & Livio, 1996). Analogous X-ray dips are seen in several Galactic X-ray binaries (e.g., White & Swank, 1982; Bonnet-Bidaud et al., 2001; Homan et al., 2003; Díaz Trigo et al., 2006) and are interpreted as evidence of a high viewing angle. Assuming that the occultation is produced by a geometrically thick structure in Keplerian rotation, we can estimate the angular extent of this feature by scaling the duration of the dipping phase to the binary period of ULX-1. If the period is 6 d (see Section 4.2), the occulting structure spans ; for a 13 d period, .

4.1.2 ULX-2 eclipse in the Chandra data

In the same set of Chandra observations, we also discovered one eclipse in ULX-2, in observation 13813 (Figure 3, bottom panel). The abrupt nature of the transition from low to high count rates once again suggests that we are looking at an occultation by the companion star. Remarkably, the egress from the eclipse of ULX-2 happens only 8ks later than the egress from the ULX-1 eclipse, at MJD 56180.30 and 56180.20, respectively (cf. bottom and top panels of Figure 3). The small but significant time difference guarantees that the two count-rate jumps seen in the two ULXs are not instrumental anomalies but real physical events. Moreover, we did extensive checks on other bright X-ray sources in the same ACIS-S3 chip, and found that none of them shows similar jumps around that time; this also rules out instrumental problems. We also examined the light-curves of ULX-2 in all other Chandra observations, including those where eclipses or dips were found in the light-curve of ULX-1 (bottom panels of Figures 2, 4, and 5). We found no other unambiguous eclipses or deep dips.

ULX-2 does show significant intra-observational variability in ObsID 13815. Throughout the 67-ks observation, the source displays a much lower count rate than its average out of eclipse count rate, in both the soft and the hard band (Table 3 and Figure 6). The count rate further decreases during that Chandra epoch, until it becomes consistent with a non-detection at the end of the observation. The decrease is slow enough (compared with the eclipse in ObsID 13813, Figure 3) to rule out a stellar occultation. We do not have enough evidence or enough counts to test whether this flux decrease is due to intrinsic variability of ULX-2, or to an increased absorption by colder material in the outer disk. As usual, we checked the behaviour of ULX-1 and other bright sources in ObsID 13815 to ascertain that the lower count rate seen from ULX-2 is not an instrumental problem.

4.1.3 ULX-2 eclipse in the XMM-Newton data

We then searched for possible eclipses of either ULX-1 or ULX-2 during the XMM-Newton observations. Due to the poorer spatial resolution of EPIC relative to ACIS-S, ULX-1 and ULX-2 are not completely resolved by XMM-Newton; however, the point spread function in the combined EPIC MOS1+MOS2 images is clearly peanut-shaped, consistent with the position and relative intensity of the two Chandra sources, and the upper source (ULX-2) has significantly harder colors (Figure 7, top panel). Firstly, we extracted and examined background-subtracted EPIC light-curves for the combined emission of the two unresolved sources, for each XMM-Newton observation. Because the two sources have comparable count rates (Tables 2 and 3), an eclipse in either source would cause the observed count rate to drop by a factor of 2. This is the scenario we find in observation 0303420101: there is an apparent increase in the observed EPIC-MOS count rate by a factor of 2, some 22 ks from the start of the observation, which we tentatively interpret as the egress from an eclipse (Figure 8, top panel), superposed on short-term intrinsic variability. Unfortunately we cannot use EPIC-pn data for this crucial epoch, because the source falls onto a chip gap. To quantify the step change in the count rate between the first and second part of the observation (green and blue datapoints in Figure 8), we performed a Kolmogorov-Smirnov (KS) test on the two distributions of datapoints, to determine whether they are drawn from different populations. We find a KS statistic of and p-value of , suggesting that the two sections of the light-curve are indeed statistically different. The average MOS1+MOS2 net count rate in the “eclipse” part of the light-curve is ct s (90% confidence limit), while in the “non-eclipse” part it is ct s. Having ascertained from the X-ray light-curve that ObsID 0303420101 probably includes an eclipse, we extracted MOS1+MOS2 images from the low-rate and high-rate sections of that observation, and confirmed (Figure 7) that in the low-rate interval, the emission from ULX-2 is missing.

We extracted and inspected the light-curves of every other XMM-Newton observation. No eclipses of ULX-1 and no further eclipses of ULX-2 were detected; however, several of those observations are much shorter than the typical Chandra observations, and the background count rate is much higher in the EPIC cameras. Thus, ruling out the presence of an eclipse as opposed to intrinsic variability is no easy task in some of the XMM-Newton observations.

4.2. Constraints on the binary period of ULX-1

We noted (Table 2 and Section 4.1.1) that for ULX-1, two fractions of eclipses are seen 12 days apart, in ObsID 13813 and ObsID 13814. The egress from the eclipse in ObsID 13813 occurs at MJD 56180.21; the ingress into the eclipse in ObsID 13814 occurs at MJD 56191.64. This enables us to place some constraints on binary period, which must be,

 P≈(11.43+eclipse duration)n days, (1)

with . To refine this constraint, we take into account that the minimum duration of an eclipse is 90 ks (1.0 days), as observed in ObsID 13814. We also know that the maximum duration of an eclipse is 150 ks (1.7 days) as this is the time between the start of the eclipse in ObsID 13814 and the start of the next observation, ObsID 13815, which has no eclipse. Assuming the shortest possible duration of the eclipse implies a binary period of 12.5 days. If we use the maximum eclipse time, 1.7 days, the binary period is 13.1 days.

We tested a range of eclipse durations and binary periods, to determine which combination of parameters is consistent with the observed sequence of eclipses/non-eclipses in our Chandra observations. Based solely on the minimum duration of an uninterrupted non-eclipse phase (160 ks) and the minimum duration of the eclipse (90 ks), the minimum acceptable binary period, from Equations (1), is 230 ks 2.7 d (that is, ). However, if the binary period were 3 days, the eclipse found during ObsID 13813 implies that another eclipse should be detected in ObsID 13812. The start of ObsID 13812 is only 2.56 days after the end of the eclipse in ObsID 13813. We do not find an eclipse in ObsID 13812, and this rules out a period of 3 days. Moreover, an eclipse time 90 ks over a period of about 3 days would imply that ULX-1 should be in eclipse 30% of the time. A Roche-lobe filling donor star can eclipse a point-like X-ray source for such a long fraction of the orbit only for mass ratios a few 100 (Fig. 2 in Chanan et al., 1976), which is impossible for any combination of compact objects and normal donor stars. Next, we consider the possibility that in Equation (1), which corresponds to a period range between 4.08 and 4.38 days. In this case, too, we would have seen at least part and more likely all of an eclipse in ObsID 13812, which is not the case: this rules out the case, too. Therefore, the only two acceptable options for the binary period are (–6.5 days) or (–13 days).

We summarize the acceptable region of the period versus eclipse duration parameter space in Figure 9. We iterated over all possible eclipse durations (1.04–1.70 days, in iteration steps of 0.01 days) and for values of , and compared the predicted occurrences of eclipses with what is detected in the seven Chandra observations between 2012 September 9 to 2012 October 10. Along the line corresponding to each value of , some periods are consistent with the Chandra data (red intervals), others are ruled out (black intervals). In addition, for Roche-lobe-filling donors, eclipse durations of the binary period (dark shaded area in Figure 9) require a mass ratio at an inclination angle of 90, or at an inclination of 80 (Chanan et al., 1976). This is very implausible if the accretor is a BH, but it is acceptable for a neutron star accreting from an OB star. In the assumption that ULX-1 has a BH primary, the mass-ratio constraint further restricts the viable case to the narrow range –6.35 days, with an eclipse duration range of days. If we allow for a neutron star primary, the period can be as long as 6.55 days, corresponding to an eclipse fraction of 26%. Finally, for the case, the predicted fractional time in eclipse goes from 8% ( d, eclipse duration 1.0 d) to 13% ( d, eclipse duration 1.7 d), with mass ratios –1, more typical of a BH primary orbiting an OB star.

Based on the previous analysis, we compared the predicted eclipse fractions with the total fraction of time ULX-1 was observed in eclipse. Over all Chandra epochs, the system is seen in eclipse for a total of 135 ks out of 835 ks, equating to a total eclipsing fraction of . No eclipses of ULX-1 are significantly detected in 177 ks of XMM-Newton/EPIC observations; therefore, the combined eclipse fraction observed by Chandra plus XMM-Newton becomes 13.3%. This is slightly lower than the predicted time in eclipse in the case of (fractional eclipse duration : Figure 9). Conversely, for the case of , the observed time in eclipse is slightly larger than expected (between 8% and 13%). We do not regard such discrepancies as particularly significant, because of the limited and uneven sampling of the system; we may have been slightly lucky or slightly unlucky in catching ULX-1 during its eclipses.

We also note that dips in the X-ray flux can sometime provide phasing information in binary systems, if they are caused by bulging, denser material located where the accretion stream splashes onto the disk. For example, regular dips at phases 0.6–0.7 are sometimes seen in low-mass X-ray binaries (e.g. EXO 0748-676; Lubow, 1989; Homan et al., 2003), and other Roche-lobe overflow systems. We have already mentioned (Section 4.1.1) that ULX-1 shows a dip in ObsID 13812. A second possible dip can also be seen in the full light-curve (Figure 10) at the start of the final observation, ObsID 13816. Although only detected in a single 1000-s bin (the first 1000 s of the observation), this drop in flux appears to be intrinsic to ULX-1, as other nearby sources do not show this feature and there are no instrumental problems in those first 1000 s. Intriguingly, both dips appear to be at the same phase with respect to the preceding eclipses (i.e., 3.5 days after the eclipse), which strengthens our confidence that the second dip is also real. For a binary period 6 days, the dips would be at phase 0.6; for the alternative period range 12.5–13 days, the dips would be at phase 0.25–0.30.

Finally, two eclipses were found for ULX-2. Unfortunately, the large time interval between the two eclipses seen by Chandra in 2012 September and by XMM-Newton in 2006 May precludes any attempt to constrain the binary period. All we can say is that the total fraction of time spent in eclipse by ULX-2 in our Ms Chandra plus XMM-Newton dataset is 7%. Since the minimum eclipse duration is 48 ks (0.55 d), we expect the binary period to be 10 d.

4.3. Hardness ratios in and out of eclipses

Residual emission is detected at the position of ULX-1 during eclipses (Table 2). This is particularly evident in ObsID 1622, with a residual eclipse count rate 10% of the average out of eclipse count rate. It is also marginally significant in ObsIDs 13813 and 13814. The reason why the residual emission appears less significant in ObsID 13813 and 13814 than in ObsID 1622 is likely because of the decreased sensitivity of ACIS-S in the soft band between 2001 and 2012 (Plucinsky et al., 2004). We stacked the eclipse intervals from all three ULX-1 eclipses and display the resulting 135-ks ACIS-S X-ray-color image in Figure 11 (top panel). The residual emission of ULX-1 is centred at the same coordinates as the out of eclipse emission, and is unresolved, but appears softer (most photons below 1 keV). We also show (Figure 11, bottom panel) the 48-ks ACIS-S image corresponding to the only Chandra eclipse of ULX-2; the signal-to-noise ratio is lower, but there is significant residual emission for ULX-2 in eclipse, as well.

To quantify the colors and the color differences in and out of eclipse, we determined the hardness ratio between the net count rates in the 1.2–7 keV band and in the 0.3–1.2 keV band (i.e., the same bands used in our light-curve plots). It appears (particularly in ObsID 13813) that ULX-1 is softer in eclipse than out of eclipse (Figure 12 and Table 4).The difference becomes more significant when we compare the hardness ratio of the stacked eclipse data (Table 4) with that of the stacked out of eclipse ones. For ULX-2, we cannot identify significant color differences in and out of eclipse, because of the short duration of the lone detected Chandra eclipse. Our hardness ratio study also clearly shows (Table 4 and Figure 12) that ULX-1 is always softer than ULX-2, both in and out of eclipse.

Another difference between the two ULXs is their degree of hardness ratio variability from epoch to epoch in the Chandra series. For ULX-1, all the 2012 observations are consistent with the same hardness ratio (Table 4). The source appears softer in ObsIDs 354, 1622 and 3932; however, this is misleading because such observations were taken in Cycle 1, Cycle 2 and Cycle 4, respectively, when ACIS-S was more sensitive to soft photons (Plucinsky et al., 2004). A rough way to account for this effect is to assume simple power-law models and use the Chandra X-Ray Center online installation of PIMMS (Version 4.8) to convert the observed count rates into “equivalent” count rates that would have been observed in Cycle 13 (year 2012) when all the other observations took place. A more accurate conversion from observed count rates to Cycle 13-equivalent count rates requires proper spectral modelling in the various epochs. The corrected count rates listed in Table 4 and plotted in Figure 12, for both ULX-1 and ULX-2, were obtained with the latter method, after we carried out the spectral analysis discussed in Section 4.4; the best-fitting spectral models were convolved with response and auxiliary response functions of the detector at different epochs, to determine the predicted count rates. Inspection of the corrected count rates confirms that the hardness of ULX-1 is approximately constant; instead, that of ULX-2 is intrinsically variable from epoch to epoch.

Based on the observed hardness of the residual eclipse emission of ULX-1, we more plausibly interpret it as thermal-plasma emission, for the purpose of converting count rates into fluxes and luminosities. Assuming a temperature 0.5 keV, and using again the online PIMMS tool, we estimate a residual ULX-1 luminosity erg s in the 0.3–8 keV band. Again, this is only a simple, preliminary estimate. We will present a more accurate estimate of the residual emission based on spectral fitting, and we will discuss its physical origin, after we carry out a full spectral modelling of ULX-1 (Section 4.4.1). We do not have any constraints on plausible models for the residual eclipse emission of ULX-2; however, a selection of thermal-plasma and power-law models also give typical luminosities erg s. What is clear is that it is harder than the residual emission of ULX-1.

4.4. Spectral properties

The presence of eclipses implies that both system are viewed at high inclination. Therefore, these two ULXs, although not exceptionally luminous, can help us investigate the relationship between the spectral appearance of ULXs and their viewing angles. During out of eclipse intervals, both ULXs have sufficiently high count rates for multi-component spectral fitting. Here, we present the results of spectral fitting to the three longest Chandra observations: ObsIDs 13812 (158 ks), 13813 (179 ks) and 13814 (190 ks), taken between 2012 September 9–20. For each source, we fitted the three spectra simultaneously, keeping the intrinsic absorption column density and the parameters of any possible thermal-plasma components locked between them but leaving all other model parameters free. The reason for this choice is that we are assuming for simplicity that cold absorption and thermal-plasma emission vary on timescales longer than a few days, while the emission from the inner disk and corona may change rapidly. In addition, we assumed a line-of-sight absorption column cm (Dickey & Lockman, 1990; Kalberla et al., 2005).

4.4.1 Spectral models for ULX-1

The first obvious result of our modelling (Table 5) is that the spectrum of ULX-1 is intrinsically curved, not well fitted by a simple power-law () regardless of the value of . Therefore, we tried several other models, roughly belonging to two typical classes: disk-dominated models, in which the disk is responsible for most of the emission above 1 keV and for the high-energy spectral curvature; and Comptonization-dominated models, in which the disk (or other thermal component) provides the seed photon emission below 1 keV, and a cut-off power-law or Comptonization component provides the bulk of the emission above 1 keV. Much of the debate in the literature about the spectral classification of ULXs reduces to the choice between these two interpretations (e.g., Gladstone et al., 2009; Feng & Soria, 2011; Sutton et al., 2013). Finally, we tested whether the addition of a thermal-plasma emission component improves the fit: the justification for this component is that some ULXs (especially those seen at high viewing angles) may show emission features in the 1 keV region (Middleton et al., 2014, 2015b).

We started by fitting single-component disk models: TBabs TBabs diskbb for a standard disk (Mitsuda et al., 1984; Makishima et al., 1986), and TBabs TBabs diskpbb for a slim disk (Kubota et al., 2005). They fare relatively better ( and , respectively) than a power-law model, but they are still not good fits. They also require a surprisingly low peak color temperature, –0.7 keV; this is inconsistent with the disk temperatures expected near or just above the Eddington limit (–1.3 keV: e.g., Kubota & Makishima, 2004; Remillard & McClintock, 2006), and would require a heavy stellar-mass BH (as we shall discuss later).

Next, we tried adding a power-law component to the disk model: TBabs TBabs (diskpbb powerlaw). This is probably the most commonly used model in the literature for the classification of accretion states in stellar mass BHs (despite the interpretation problems caused by the unphysically high contribution of the power-law component at low energies). The quality of the fit improves slightly () but there are still significant systematic residuals. One source of fit residuals is the high-energy downturn. By using instead a TBabs TBabs (diskbb cutoffpl) model, we obtain a substantially better fit (), with an F-test statistical significance . The presence of a high-energy downturn is of course one of the main spectral features of ULXs (Stobbart et al., 2006), compared with stellar-mass BHs in sub-Eddington states. However, in this case the best-fitting cut-off energy keV, much lower than the typical 5-keV high-energy cutoff seen in other ULXs (Gladstone et al., 2009); this is quantitative evidence that the spectrum of ULX-1 is extremely soft compared with average ULX spectra. Fitting a cutoff power-law alone (without the disk component) gives a ; from this, we verify that an additional soft thermal component is significant to . There is still a third source of fit residuals, at energies around 1 keV, which we will discuss later.

The successful models discussed so far are phenomenological approximations of physical models; therefore, we fitted several alternative Comptonization models that produce a soft excess and a high-energy downturn: TBabs TBabs (diskbb comptt) (Titarchuk, 1994), TBabs TBabs (bb comptt), and TBabs TBabs diskir (Gierliński et al., 2009). They provide moderately good fits, with (Table 5). In this class of Comptonization models, the disk component is used as the source of seed photons, and the electron temperature sets the location of the high-energy cutoff. For ULX-1, typical seed photon temperatures are keV, and the range of electron temperatures in the Comptonizing region is –1.2 keV, with optical depths –9. The electron temperature of the Comptonization region is substantially cooler than in most other two-component ULXs (where –3 keV: Gladstone et al. 2009). This is the physical reason why ULX-1 appears as one of the softest sources in its class, with an unfolded spectrum peaking at 1 keV (c.f. the ULX classification of Sutton et al., 2013). The optically-thick thermal continuum component used as seed to the Comptonization models can be equally well modelled with a disk-blackbody or a simple blackbody, given its low temperature at the lowest edge of the ACIS-S sensitivity. Its direct flux contribution to the observed spectrum is small, although difficult to constrain precisely, because of the low number of counts at very soft energies. Individual fits to the three longest observations with a diskir model suggest a direct disk contribution an order of magnitude lower than the Comptonized component.

The main reason why none of the smooth continuum models described above are really good fits is the presence of strong residuals (F-test level of significance 99.99%) below and around 1 keV. A single-temperature mekal component is not sufficient to eliminate the residuals. Instead, two mekal components with temperatures keV and keV significantly improve the fits (Figure 13), providing for the Comptonization models and for the disk models. In the latter case, the temperature profile index is (“broadened disk”), favouring the slim disk over the standard disk model.

The presence of soft X-ray residuals and the strong continuum curvature are robust and independent of the choice of cold absorption model. We also tried combinations of neutral and ionized absorbers (tbabs varabs), but they do not reproduce the strong residual feature at energies 0.8–1.0 keV. Lower-energy residuals at 0.5–0.6 keV are relatively less constrained because of the degraded sensitivity of ACIS-S at low energies, rather than because of intrinsic absorption. For all our Comptonization-type and disk-type models, the intrinsic cold absorption is a few cm and in most cases, consistent with 0 within the 90% confidence level.

Disk models with additional thermal plasma emission produce equally good values as Comptonization models with thermal plasma emission (Table 5). Therefore, it is worth examining in more details whether disk models are physically self-consistent, and what their physical interpretation could be. Let us start with a standard disk, to account for the possibility that ULX-1 is a rather massive BH accreting at sub-Eddington rates. We may be tempted to discard this possibility straight away, because the best-fitting temperature profile index (rather than ) is generally considered the hallmark of a disk at the Eddington accretion rate, not of a standard disk; however, it was recently shown that broadened disks with may also occur in the sub-Eddington regime with accretion rates an order of magnitude lower (Sutton et al., 2016). Therefore, we will consider that case here. In standard accretion disk models, with the inner disk fixed at the innermost stable circular orbit, there are two approximate relations between the observable quantities and , and the non-directly-observable physical properties (Eddington-scaled mass accretion rate) and (BH mass):

 L ≈ 1.3×1039˙mM10  erg s−1 (2) kTin ≈ 1.3(˙m/M10)1/4  keV (3)

(Shakura & Sunyaev, 1973), where is the BH mass in units of 10 . Equation (2) is simply the luminosity as a fraction of Eddington, in the radiatively efficient regime; Equation (3) follows from the relation . For ULX-1, the best-fitting peak temperature is keV, and the luminosity erg s (as we shall discuss later). From Equations (2) and (3), this would correspond to a BH mass –50 at accretion rates –0.4. Even if we allow for the possibility of in a sub-Eddington disk, the presence of strong line residuals in the soft X-ray band is another, stronger piece of evidence against the standard disk model; it is instead indicative of Eddington accretion and associated outflows (Middleton et al., 2015b; Pinto et al., 2016a). For these reasons, we dis-favour the sub-Eddington standard disk model for ULX-1.

Next, we test the self-consistency of the slim disk model. In the super-Eddington regime, the disk luminosity is modified as

 L≈1.3×1039(1+0.6ln˙m)M10  erg s−1 (4)

(Shakura & Sunyaev, 1973; Poutanen et al., 2007). A slim disk is no longer truncated exactly at the innermost stable circular orbit: the fitted inner disk radius can be approximated by the empirical scaling (Watarai et al., 2000). As a result, the disk luminosity becomes . Substituting from Equation (4), and matching the normalization to that of the sub-Eddington case (Equation (3)), we obtain:

 kTin≈1.3(1+0.6ln˙m)1/2M−1/410  keV. (5)

A proper treatment of the observed properties of a slim disk requires additional parameters such as the viewing angle and the BH spin (Vierdayanti et al., 2008; Sa̧dowski, 2009; Vierdayanti et al., 2013); however, Equations (4) and (5) are already good enough as a first-order approximation to test the consistency of a slim disk model for ULX-1. From Equation (5), we find that a slim disk peak temperature keV requires a BH mass (confirmed also by the numerical results of Vierdayanti et al. 2008); but such BH would be far below Eddington at the observed luminosity erg s, contrary to our initial slim disk assumption. Therefore, the slim disk model cannot be self-consistently applied to the spectrum of ULX-1.

Based on those physical arguments, we conclude that the best fits in all three epochs are obtained with a Comptonization model, with the addition of multi-temperature optically-thin thermal-plasma emission. From the fits statistics alone, we cannot rule out a broadened disk model, with a heavy stellar-mass BH at sub-Eddington accretion rates; however, the observed presence of strong soft X-ray residuals points to an ultraluminous regime. We list the best-fitting parameters of two equivalent Comptonization models in Tables 6, 7; for comparison, we also list the best-fitting parameters of the broadened-disk model (Table 8).

4.4.2 Continuum and line luminosity of ULX-1

The unabsorbed X-ray luminosity is related to the absorption-corrected flux by the relation , where is the distance to the source and is the viewing angle, when the emission is from a (sub-Eddington) standard disk surface, and for a spherical or point-like emitter. We do not have direct information on the geometry of the emitting region in ULX-1; however, analytical models and numerical simulations of near-Eddington and super-Eddington sources predict mild geometrical beaming, that is, most of the X-ray flux is emitted along the direction perpendicular to the disk plane (Kawashima et al., 2012; Jiang et al., 2014; Sa̧dowski & Narayan, 2016), and the emission should appear fainter and down-scattered in a disk wind when a source is observed at high inclination (as in our case, given the presence of eclipses). Therefore, we choose to use the simplest angle-dependent expression for the luminosity . We also identify for simplicity (and in the absence of conflicting evidence) the viewing angle to the plane of the inner disk with the inclination angle of the binary system, which we have constrained to be high from the presence of eclipses; that is, we neglect the possibility of a warped, precessing disk. Instead, we estimate the unabsorbed X-ray luminosity of the thermal plasma components as , because we assume that the distribution of hot plasma is quasi-spherical above and beyond the disk plane, and its emission is approximately isotropic. With those caveats in mind, we estimate an emitted 0.3–8.0 keV luminosity of the two-temperature thermal-plasma component erg s (Tables 6, 7, 8), essentially all below 2 keV. Assuming , we then estimate the total (i.e., thermal plasma plus continuum) 0.3–8.0 keV luminosity as erg s, during the three longest Chandra observations, regardless of whether the continuum is fitted with a Comptonization model or with a slim disk.

To increase the signal-to-noise ratio of the thermal-plasma emission component, we extracted and combined (using CIAO’s specextract tool) the spectra and responses of all ten Chandra observations (Table 1), for a grand total of 700 ks out of eclipse. We fitted the resulting spectrum with the same smooth continuum models (Comptonization and slim-disk models) used for the three long spectra: regardless of the choice of continuum, strong systematic residuals appear at energies around 1 keV. As a representative case, we show the residuals corresponding to a comptt fit (Figure 14, top panel); for this continuum-only model, . When two mekal components are added to the continuum (as we did for the spectra of ObsIDs 13812, 13813 and 13814), the goodness-of-fit improves to . We tried introducing a third mekal component, and obtained a further improvement to the fit (significant to the 99.5% level), down to . We show the 3-mekal model fit and its residuals in Figure 14 (bottom panel), and list the best-fitting parameters in Table 9. The best-fitting mekal temperatures are keV, keV, and keV. (Instead, adding a third mekal component to the best-fitting model for ObsIDs 13812, 13813 and 13814 does not significantly improve that fit.) We estimate an unabsorbed 0.3–8.0 keV luminosity of the three-temperature thermal-plasma component erg s. This is moderately higher than the value we estimated for a two-temperature model, because now part of the emission at energies 2 keV is also attributed to optically-thin thermal plasma. The total (continuum plus thermal plasma) unabsorbed luminosity in the 0.3–8.0 keV band is erg s, consistent with the luminosity estimated in ObsIDs 13812, 13813 and 13814. Alternatively, we replaced the three mekal components with a cemekl, which is a multi-temperature thermal-plasma model with a power-law distribution of temperatures. The best-fitting cemekl diskbb comptt model has , maximum temperature keV, thermal-plasma luminosity erg s, and total luminosity erg s.

As noted earlier, the role of the disk component in this class of models is to provide the seed photons for the Comptonization component, as well as a soft excess due to the fraction of disk photons that reach us directly. A best-fitting seed temperature keV and inner-disk size km are consistent with the characteristic temperatures and sizes of the soft thermal components seen in other ULXs (e.g., Miller et al. 2004; Stobbart et al. 2006; Kajava & Poutanen 2009). The direct luminosity contribution of the disk in the 0.3–8 keV band is erg s.

Finally, we inspected the spectral emission in eclipse. We extracted a combined spectrum of the eclipse intervals in ObsIDs 1622, 13813 and 13814. Although we only have 60 counts, the energy distribution of the counts is similar (Figure 15) to the thermal-plasma emission component out of eclipse, rather than to the continuum emission. After rebinning the eclipse spectrum to 1 count per bin, we applied the Cash statistics (Cash 1979) to fit the normalization of the same two mekal components previously found in the out-of-eclipse spectra of ObsIDs 13812, 13813 and 13814 (temperatures fixed at keV and keV). We find a C-stat value of 49.5 over 54 degrees of freedom for the best-fitting model. The emitted luminosity erg s, consistent with our previous simpler estimate based on count rates (Section 4.2). We then let the temperatures free, but did not obtain any significant improvement to the quality of the fit (C-stat value of 49.2 over 52 degrees of freedom). Nor do we improve the fit by adding a third mekal component.

4.4.3 Spectral models and luminosity of ULX-2

As we did for ULX-1, we started by fitting the spectra of ULX-2 during Chandra ObsIDs 13812, 13813 and 13814 with a simple power-law model (Table 5). The fit is good (), but there are residuals consistent with a high-energy downturn. The best-fitting power-law index is ; however, this value may be an over-estimate if the high-energy steepening is not properly accounted for. Hence, we re-fitted the spectrum with a cutoff power-law (TBabs TBabs cutoffpl), and found that the fit is significantly improved: , with an F-test significance 99.99% with respect to the unbroken power-law. The power-law index below the cutoff is and the characteristic energy of the cutoff is () keV. This is evidence that the spectrum of ULX-2 is significantly curved. Therefore, as we did for ULX-1, we tried a series of models suitable to curved spectra: disk models and Comptonization model.

Among disk models, we find that a broadened disk is a significantly better fit (F-test significance 99.99%) than a standard disk; a TBabs TBabs diskpbb model provides (Table 5). The peak disk temperature 1.4–2.0 keV and , perfectly in line with the expected values for a mildly super-Eddington slim-disk model around a stellar-mass BH. The best-fit parameters can be found in Table 10; the model is illustrated in Figure 16. The diskpbb normalization, , translates into a characteristic inner disk radius

 Rin≈3.18K1/2d10kpc(cosθ)−1/2  km, (6)

using the conversion factors suitable for slim-disk models (Vierdayanti et al., 2008); is the distance in units of 10 kpc. A feature of super-critical slim disks is that is located slightly inside the innermost stable circular orbit (Watarai & Mineshige, 2003; Vierdayanti et al., 2008). When this correction is taken into account, the mass of a non-rotating BH can be estimated as . Characteristic radii 29–56 km are consistent with all the three long Chandra observations considered here (Table 10). For , this corresponds to characteristic masses –18, consistent with the observed mass distribution of Galactic BHs (Kreidberg et al., 2012). For a range of viewing angles , the corresponding BH mass range becomes –25. The emitted luminosity in the 0.3–8.0 keV band is erg s (assuming again a viewing angle ) and the bolometric disk luminosity is erg s 1–3 for the range of BH masses estimated earlier. In this model, ULX-2 would be classified as a broadened-disk ULX in the scheme of Sutton et al. (2013).

Although we favour the slim disk model because of its self-consistency, we cannot rule out the possibility that ULX-2 is fitted by a Comptonization model (Table 11): for example, TBabs TBabs (diskbb comptt) yields , statistically equivalent to the slim disk model (Table 5), with electron temperatures –1.5 keV and optical depth –13 (slightly hotter and more optically thick than the best-fitting comptt models in ULX-1). Similar values of and are also obtained from other Comptonization models such as diskir.

Regardless of the model, the unfolded spectrum peaks at 5 keV, similar to the sources classified as hard ultraluminous by Sutton et al. (2013). The original definition of the hard ultraluminous regime requires also the presence of a soft excess. In our spectra, it is difficult to constrain the significance of a direct soft emission component (in addition to the Comptonized component or the cutoff power-law) because of the low sensitivity of ACIS-S below 0.5 keV. When we fit the spectrum with a diskbb comptt model, we find that no more than 50% of the flux in the 0.3–1.0 keV band is in the direct diskbb component (90% upper limit), but the diskbb normalization is also consistent with 0 within the 90% confidence limit. Regardless of classification semantics, it is clear that ULX-2 has a hard spectrum in the Chandra band, with a high-energy curvature.

No significant residuals are found at 0.8–1 keV in the individual spectra from ObsIDs 13812, 13813 and 13814; however, in at least one observation (ObsID 13812), the spectrum shows two emission features with 90% significance at keV and keV. Similar lines are typically found in thermal plasma emission. They are usually interpreted as emission from a blend of Mg XI lines at 1.33–1.35 keV, and from a Si XIII line at 1.84 keV (with the likely additional contribution of slightly weaker Mg XII lines at 1.75 keV and 1.84 keV). To investigate these and possible other emission features, we extracted a combined Chandra spectrum of ULX-2 from all ten observations, as we did for ULX-1. We fitted the combined spectrum with a diskpbb model, and obtain an excellent fit, (Table 12). The significance of the two candidate emission features seen in ObsID 13812 fades to 90% in the combined spectrum (Figure 17). Adding thermal-plasma components to the combined spectrum does not produce any significant improvement. The characteristic disk temperature keV, radial temperature index and inner-disk radius km (Table 12) are consistent with those expected for a super-critical disk, and with the values obtained from the individual fits to ObsIDs 13812, 13813 and 13814. The corresponding range of BH masses is –20, for a non-rotating BH and a viewing angle .

Finally, we examined the spectrum of ULX-2 in eclipse (Figure 18). It appears different from what is seen in ULX-1: there is no evidence of a bimodal distribution of counts and it is not possible (from the few counts available) to determine whether the eclipse emission has the same origin as the out of eclipse continuum (e.g., a small fraction of the direct emission scattered into our line-of-sight by an extended corona), or comes from thermal-plasma at higher temperatures or from bremsstrahlung emission.

5. Discussion

5.1. Two eclipsing ULXs in one field: too unlikely?

Luminous stellar-mass BH X-ray binaries or ULXs with X-ray eclipses are very rare sources. SS 433 in the Milky Way shows eclipses of its X-ray emission caused by the donor star on a 13.1-d binary period (e.g., Stewart et al., 1987; Fabrika, 2004; Brinkmann et al., 2005; Kubota et al., 2010; Cherepashchuk et al., 2013; Marshall et al., 2013). Unlike M 51 ULX-1 and ULX-2, SS 433 does not appear as luminous as a ULX because the direct X-ray emission from the inner disk/corona region is already occulted from us. Its donor star periodically eclipses the thermal bremsstrahlung radiation ( erg s) from the base of the jet. The first unambiguous eclipsing behaviour in a candidate BH X-ray binary outside the Milky Way was found in IC 10 X-1, located in a Local Group dwarf galaxy, with a Wolf-Rayet donor star, a binary period of 1.45 days, and an X-ray luminosity erg s (Prestwich et al., 2007; Laycock et al., 2015a; Steiner et al., 2016). For IC 10 X-1, it is still disputed whether the accreting compact object is a BH or a neutron star (Laycock et al., 2015b). Outside the Local Group, NGC 300 X-1 ( erg s; binary period 33 hr) shows X-ray dips, consistent with occultation from geometrically thick structures in the outer disk, or absorption in the wind of the donor star, but not with true eclipses (Binder et al., 2015). A strong candidate for a true eclipse is the sharp dip in the Swift/X-Ray Telescope flux recorded once from the ULX P13 in NGC 7793, at an orbital phase consistent with the inferior conjunction of its supergiant donor star (Motch et al., 2014); however, there is no further confirmation of that single monitoring datapoint at subsequent epochs. Thus, we argue that the two M 51 ULXs discussed in this paper are the first unambiguous eclipsing sources observed at or near the Eddington regime.

It is rare enough to find two such bright sources projected close to each other in what is not a particularly active starburst region: it is obviously even stranger that both of them show eclipses. Therefore, we tried to assess the statistical significance of this finding. Firstly, we assume that any distance between ULX-1 and ULX-2 not in the plane of M 51 is negligible and thus only consider the separation. We want to discover the chances of finding two randomly distributed, luminous X-ray binaries in the same galaxy within of each other, both having inclination angles 80. Assuming for example 10 ULXs with erg s in the same spiral galaxy within a radius of (an over-estimate of the real number of ULXs detected in local-universe galaxies), we used a Monte-Carlo simulation, placing ULXs at random and recording the number of occurrences in which two ULXs were found within a radius of ; for 10 million trials, we find . The probability of finding two nearby ULXs then has to be multiplied by the probability that they both show eclipses. Assuming no preferential orientation angle, the likelihood of finding a ULX with a viewing angle, for example, is . However, if the orientation is too close to , the direct X-ray emission is likely blocked by the outer disk and the source would not appear as a ULX. The thickness of the disk in ULXs is unknown (likely a few degrees), and the minimum angle that produces eclipses is model-dependent, as a function of the ratio between stellar radius and binary separation , namely . For plausible distributions of such quantities, (Pooley & Rappaport, 2005). The final probability becomes a few : we only expect this to happen once every few hundred major galaxies.

5.2. Spectral properties: broadened disks, Comptonization and thermal plasma

Luminosity (or more precisely, mass accretion rate) and viewing angle are thought to be the main parameters that determine the observational appearance of ULXs in the X-ray band (e.g., Sutton et al., 2013; Middleton et al., 2015a; Urquhart & Soria, 2016); disentangling and quantifying their roles is still an unsolved problem. M 51 ULX-1 and ULX-2 have approximately the same luminosity and inclination angle (as they both show eclipses); however, they spectral appearance is substantially different. ULX-1 has soft colors in the Chandra band and is well modelled by a soft thermal component (blackbody or disk-blackbody) plus Comptonization; ULX-2 has hard colors and is well modelled by a slim disk with –2.0 keV (hotter than a standard disk). Also, ULX-1 has significant line residuals around 1 keV (consistent with thermal-plasma emission), which are not seen in ULX-2. Thus, we propose that there are other physical parameters that determine the spectral appearance of a ULX in addition to Eddington ratio and viewing angle. We also discovered that ULX-1 has strong radio and optical evidence of a jet (as we will discuss in a separate paper; Soria et al., in prep.) while ULX-2 does not; understanding the relation between outflow structure and spectral appearance remains a key unsolved problem.