CG X-1: An Eclipsing Wolf–Rayet ULX in the Circinus Galaxy

CG X-1 was observed 24 times by Chandra's Advanced CCD Imaging Spectrometer (ACIS) from 2000 to 2010 (Table 1). Eight of the observations were centered on the back-illuminated S3 chip with no gratings. The other 16 were taken with the High Energy Transmission Grating (HETG) in front of the ACIS-S detectors. We downloaded the data from the public archive and reprocessed them with the task chandra_repro that is provided by the Chandra Interactive Analysis of Observations (ciao) software package, Version 4.9 (Fruscione et al. 2006). We corrected the photon arrival times to the solar system barycenter using the ciao task axbary.

For the non-grating data, we used dmcopy to extract images in the 0.3–8 keV band from the event files of the individual observations. We then applied the source-finding task wavdetect to the images of individual observations, which allowed us to determine the coordinates of CG X-1 and to identify an ellipse around CG X-1 that contains ≈99% of the source counts (wavdetect parameter ellsigma = 3). The typical semimajor axis of the elliptical source region is ≈13. We used that ellipse as the source extraction region for our timing and spectral analysis. For the local background regions, we used an elliptical annulus that is centered at the position of the source The inner and outer sizes of the annulus were 2 and 4 times the size of the source ellipse, respectively. We extracted background-subtracted light curves of each observation in the 0.3–8 keV band, with the task dmextract. We used specextract to create spectra and associated background and response files, for each observation. We used the specextract option "correctpsf = yes" for aperture correction, and we grouped the spectra to a minimum of 15 counts per bin, for subsequent χ² fitting.

For the HETG observations, we only extracted light curves (also filtered to the 0.3–8 keV band), from the zeroth-order images. The source extraction region was a circle with a radius of 13, centered at the averaged position determined from the individual non-grating images. The background region is a circle with radius of 16'' in the northeast of CG X-1, free of point sources. Subsequent analysis of spectra was done with standard tool xspec from the heasoft package, Version 6.21 (Blackburn 1995).

XMM-Newton observed CG X-1 on six occasions before 2018 October, with the European Photon Imaging Camera (EPIC) as prime instrument. For this work, we used the first five observations (Table 1), which are already available for download from the public archive. We reduced the data with the XMM-Newton Science Analysis System (sas) version 16.7.0. In particular, we ran the standard tasks cifbuild, odfingest, epproc and emproc to obtain Processing Pipeline System (PPS) products from the Observation Data Files (ODFs), for EPIC-pn and EPIC-MOS. The sas task barycen was applied for barycentric corrections.

For each observation, we used xmmselect to create preliminary light curves of the whole field of view above 10 keV, for the purpose of flagging and removing time intervals affected by background flares. Time intervals with a count rate less than some certain criteria were kept to create the GTI intervals. The count-rate criteria varies according to the observations and cameras, typically higher than 0.8 ct s⁻¹ for pn and 0.3 ct s⁻¹ for MOS, respectively. For the first four observations, we used the resulting GTI intervals for both light curve and spectral extractions. Instead, the most recent observation in our sample (ObsID 0780950201) was affected by moderate background flares throughout most of the exposure; for this reason, and for this observation only, we decided not to remove the high background intervals, otherwise the GTI would have been too short for meaningful analysis.

Defining a source extraction region and a local background is much more complicated for XMM-Newton than for Chandra because of the larger point-spread function (PSF) of the former. In fact, the full-width at half-maximum of the pn and MOS PSFs is approximately the same as the distance between CG X-1 and the (stronger) nuclear X-ray source of Circinus (≈15''). Hence, we had to use a small extraction radius of only 5'' for CG X-1 to reduce the nuclear source contamination from the wings of its PSF. The local background was extracted from three circles with radius of 5'', located around the Circinus nucleus and centered at approximately the same distance of 15'' from the central source. This is the best way to ensure that the nuclear PSF-wing contamination to the source extraction region is as similar as possible to its contribution to the background regions and can be effectively subtracted out. We used the same source and background regions for both the light curve and the spectral extraction.

We used "#xmmea_em && (pattern < = 12)" to filter MOS light curves and spectra. We then used "#xmmea_ep && (pattern < = 4 && flag == 0)" to filter the PN data. We extracted light curves from MOS and pn with the sas tasks evselect, followed by epiclccorr to correct for vignetting and bad pixels.¹¹ The photon energy was filtered to 0.3–8 keV.

For each observation, we extracted individual PN, MOS1 and MOS2 spectra with evselect, in 0,3–8 keV band. Response and ancillary response files were generated with rmfgen and arfgen, respectively. We then combined the MOS and pn spectra of each observation with the sas task epicspeccombine, to improve the signal-to-noise ratio of possible narrow features. We grouped the merged spectra to at least 20 counts per bin for χ² fitting. Finally, as for the Chandra data, we used xspec for spectral modeling.

CG X-1 has been observed by the (ROSAT/HRI) five times. Two of those observations have exposure times longer than 25 ks (Table 1), while the other three are shorter than 5 ks. We only used the two observations with the longest exposure time; i.e., RH702058A02 and RH702058A03, taken at various intervals between 1997 March 3 and 11, and between 1997 August 17 and September 9, respectively. We used the heasoft tool rosbary to do the barycenter correction of the event lists. The source region used for extracting the light curve is a circle with a radius of 6''. The background region is an annulus with inner radius of 6'' and outer radius of 10'' around the source. Due to the low count rate and sparse snapshots in these two observations, we folded their light curves using the best-fitting period (determined from Chandra and XMM-Newton observations), for a better statistics. The epoch of this folded light curve was chosen at the orbital cycle, where the two observations accumulated half of the total observed counts; i.e., MJD 50681.437036 for phase ϕ = 1 (the phase will be defined in Section 3.1). The main purpose for using the folded ROSAT light curve is to study the variation of the periodicity of CG X-1 over a long-term duration of more than 20 yr.

A (27 ± 0.7) ks X-ray periodicity was first discovered by Bauer et al. (2001) from Chandra/HETG zeroth-order data. It was subsequently confirmed with BeppoSAX and Chandra/ACIS-S observations (Bianchi et al. 2002; Weisskopf et al. 2004). More recently, Esposito et al. (2015) analyzed two long Chandra/ACIS-S observations from 2010 (ObsID 12823 and ObsID 12824) and obtained a more precise value of (26.1 ± 0.2) ks. To reduce the error even further, for this work we reanalyzed the whole archival set of 20 Chandra observations and five XMM-Newton observations longer than about half of the period (≈13 ks): they cover a time span of 19 yr (Table 1). (We did not use an additional four archival Chandra observations because their duration was only ≲0.25 times the period.)

For each Chandra or XMM-Newton observation, we generated a 0.3–8 keV background-subtracted light curve, and then rebinned each light curve to 500 s bins (Figure 2). We then applied the phase dispersion minimization (PDM) method to compute the best-fitting period (Schwarzenberg-Czerny 1997; Stellingwerf 1978). The PDM method is more accurate if the amplitude of the modulation is constant. But the count rates measured by different instruments vary a lot because of their different instrumental responses (Figure 2). Thus, we normalized the count rate of each observation by dividing the original count rate by the maximum count rate in each observation. The value of the normalized count rate is in the range of 0–1. We managed to phase-connect all of the archival observations (Figure 2), with a best-fitting period P_best = (25970.0 ± 0.1) s: which is an improvement in precision by a factor of 2000.

Zoom In Zoom Out Reset image size

We computed six folded light curves for different subgroups of Chandra and XMM-Newton observations (Figure 3): all those taken in 2000–2001; in 2004; in 2006; in 2008; in 2010; and in 2013–2016. Each of the six folded light curves confirms an approximate description of the average profile as a sharp ingress, a short, deep eclipse, and a slow egress. Using a fourth-order Fourier model, we fitted the eclipse section (defined as the bins with a count rate less than 30% of the maximum count rate, in the phase interval between two consecutive peaks) of all the folded light curves simultaneously. We redefined phase ϕ ≡ 1 as the deepest point of the eclipse in the best-fitting model. As a first approximation, the phase of the eclipse is the same in each of the six folded light curves: we estimate that the deepest points of the six folded eclipses have only a small scatter of δϕ ≈ 0.034 around the global best-fitting phase ϕ = 1. (In Section 3.2, we will investigate this scatter further and will look for possible small changes in the period over the last 20 yr.) This confirms that the X-ray eclipse is related to stable properties of the binary system: the most plausible explanation is that ϕ = 1 corresponds to superior conjunction of the accreting compact object (i.e., when it passes behind the star). However, the profile of each individual cycle is much more irregular and variable from cycle to cycle, as we have shown (Figure 2). We will attempt to explain the irregular profiles in Section 6.1.

Zoom In Zoom Out Reset image size

In summary, assuming a constant period P_best, our best-fitting ephemeris is ϕ = 1 at MJD 50681.437036+0.300579 × N (d). We have chosen to set the reference time (N = 0) at an epoch covered by the ROSAT observations.¹² The cycle number N of subsequent observations is plotted as a label in the top left corner of each frame of Figure 2. Choosing a different zeropoint (for example defining N = 0 as the first Chandra observation) would obviously not change any physical interpretation.

The next step of our analysis was to search for possible small changes in the binary period over the two decades of observations. This presents a practical challenge. Although the time of mid-eclipse (used to define the reference phase ϕ = 1) is relatively easy to determine in all the light curves averaged over many cycles or several years (Figure 3), it is not obvious how to identify this point in any of the individual light curves (Figure 2), with eclipse durations and egress behaviors that differ markedly from cycle to cycle. Instead, we used the eclipse ingress as a phase marker for the individual light curves. The sharp flux drop leading to the eclipse is a feature that can be unequivocally identified in most of the individual light curves. It is already clear from a cursory inspection of Figure 2 that the phase of mid-ingress differs from cycle to cycle, what we want to determine is whether this is a random scatter around a mean value or represents a systematic drift, which is possible evidence of the period derivative $\dot{P}\ne 0$.

First, we determined the average mid-time of the ingress from the stacked Chandra and XMM-Newton light curves. We used the method developed by (Hu et al. 2008), which we adapted for our purposes. We defined the count rate I₀ of the un-eclipsed high state as the average count rate from ϕ = 0.5 to ϕ = 0.8. We then introduced a parameter ${\gamma }_{i}\equiv \left({I}_{i}-{I}_{i+1}\right)/(500\,{\rm{s}})$, to characterize the local slope of the light curve between two successive bins at times t_i and t_i+1 = t_i+500 s. Finally, we defined a mid-time t_c of the ingress section of a light curve as

$$\begin{eqnarray}&&{t}_{{\rm{c}}}=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{\gamma }_{i}\left({I}_{0}-{I}_{i}\right){t}_{i}}{{\displaystyle \sum }_{i=1}^{N}{\gamma }_{i}\left({I}_{0}-{I}_{i}\right)},\end{eqnarray} \tag{ 1 }$$

where i is an integer index running from the left boundary (i = 1) to the right boundary (i = N) of the ingress interval (approximately defined as the phase interval 0.7 < ϕ < 1), and t_i and I_i are the time and count rate of the ith bin, respectively. Equation (1) is essentially a weighted average of the ingress times; the weights are the steepness of the decline times the total drop at a certain point. The reason why the average is also weighted by $\left({I}_{0}-{I}_{i}\right)$ is because we want to reduce the uncertainty in the estimate of the time at which the light curve rolls over and the ingress starts; near the beginning of the ingress, $\left({I}_{0}-{I}_{i}\right)\approx 0$. In any case, Equation (1) provides a practical definition of a light-curve feature that can be used to search for first-order period changes. Similarly, the physical width of the ingress is calculated as:

$$\begin{eqnarray}&&W={\left|\displaystyle \frac{{\displaystyle \sum }_{i=1}^{N}{\gamma }_{i}\left({I}_{0}-{I}_{i}\right){\left({t}_{i}-{t}_{{\rm{c}}}\right)}^{2}}{{\displaystyle \sum }_{i=1}^{N}{\gamma }_{i}\left({I}_{0}-{I}_{i}\right)}\right|}^{1/2}.\end{eqnarray} \tag{ 2 }$$

The mid-ingress phase for the stacked Chandra and XMM-Newton observations occurs at phase ϕ = 0.85 ± 0.03, about 3896 s before mid-eclipse. Therefore, the predicted time of mid-ingress during an arbitrary cycle number N is MJD 50681.391949+0.300579 × N(d) ≡ T₁ + P_best × N. The empirically determined location of t_c in each of the individual light curves plotted in Figure 2 is marked with a magenta line.

To investigate whether the period changes with time over a 20 yr timescale, we computed an "Observed minus Calculated" (O − C) diagram. The predicted C values of the ingress times come from the ephemeris given above (C = T₁ + P_best × N), and are assumed to have no error. We calculated the O values from the data, using Equation (1). In particular, we used one ROSAT, four Chandra/ACIS, eleven Chandra/HETG, and five XMM-Newton/EPIC observations. For the Chandra/ACIS, XMM-Newton/EPIC, and three of the Chandra/HETG observations, we had sufficient counts to measure ingress times directly for individual cycles. We stacked the remaining eight Chandra/HETG observations (separated by only two weeks), to increase the signal to noise ratio. In this case, we measured the time difference between the average ingress phase in the stacked light curve (folded on the default ephemeris) and the predicted ingress phase. We used a similar method for those observations that covered multiple orbital cycles and therefore contained more than one ingress. In those cases, a single O − C data point was used for each observation, which is defined as the average of the individual time differences for all the ingresses covered during that observation. We estimate the 1σ uncertainty (standard deviation) for each individual ingress measurement as ≈700 s. Smaller errors are associated to data points that are the average of multiple ingress times.

In an O − C diagram, if the period remains constant with time, data points are scattered along a line:

$$\begin{eqnarray}&&O-C=({P}_{0}-{P}_{\mathrm{best}})\times N,\end{eqnarray} \tag{ 3 }$$

where P₀ is the true period and P_best is the best-fitting period used to determine the values of C. If the fitted period is exactly equal to the true period, the O − C diagrams follow a horizontal line around zero. Instead, if the orbital period has a small linear change, then the data points follow a quadratic function:

$$\begin{eqnarray}&&O-C={c}_{1}+{c}_{2}\,N+(1/2){P}_{\mathrm{best}}\dot{P}\times {N}^{2},\end{eqnarray} \tag{ 4 }$$

where c₁ and c₂ are normalization coefficients, and $\dot{P}$ is the period derivative.

The resulting O − C diagram is shown in Figure 4, which visually suggests an upwards curvature. We fitted the O − C data points with a quadratic function. The best-fitting curve (dashed–dotted line in Figure 4) is O − C = 1320.6 − 0.29 × N+1.09 × 10⁻⁵ × N². The orbital period derivative is $\dot{P}$ = (8.4 ± 3.8) × 10⁻¹⁰ s s⁻¹, where the error range is the 90% confidence level (2.7σ), which corresponds to $\dot{P}/P=(10.2\pm 4.6)\times {10}^{-7}$ yr⁻¹ at the 90% confidence level. The value of $\dot{P}/P$ is significantly >0 to approximately 10σ (99.9% confidence level). As a further check, we used the F-test to evaluate the improvement from a linear fit to a quadratic fit. The low F-test probability (p ≈ 1.1 × 10⁻³) confirms that the O − C data points follow a quadratic fit (curving upwards) better than a linear fit (constant period). This supports our conclusion that we are detecting a systematic increase of the binary period over 20 yr.

Zoom In Zoom Out Reset image size

Light curves folded over several cycles (Figure 3) show an eclipse lasting from ϕ ≈ 0.9 to ϕ ≈ 1.1, followed by a slow egress. In fact, as we have already mentioned, the profiles of the individual cycles tell a more complex story (Figure 2). In some cycles, such as those observed during Chandra ObsID 356 and the first three XMM-Newton observations, the faintest phases last δϕ ≈ 0.15. At other epochs, such as Chandra ObsID 9140 and ObsID 12823, the observed flux remains close to zero for a longer time, δϕ ≈ 0.4. In yet other cases, for example cycles N = 13,781 and N = 13,782 during Chandra ObsID 10224, we see a sequence of irregular dips and flares rather than a well-defined, single eclipse. This irregular behavior cannot be explained with a simple stellar occultation of a point-like X-ray source, other factors must (also) be contributing to the duration of the periodic occultation and the variable dipping profile; as will be discussed in Section 5.1.

We investigated whether there is a relation between the properties of the eclipsing phase in a cycle and the unobscured luminosity before the ingress. However, because of the irregular dipping, in many orbital cycles it is difficult to determine when the eclipse ends (if it indeed is a true stellar eclipse) and the egress begins. Instead, we defined a different, empirical quantity that parameterizes the total duration of the most occulted phases during each orbital cycle. First, we determined the maximum count rate of a cycle as the average of the rates in the five 500 s bins with the highest count rates, in the phase range ϕ = 0.8–1.8. We then defined a duration of the occultations as the number of 500 s phase bins (N_occ) in which the count rate is lower than 40% of the maximum count rate that we earlier defined. These bins typically include the candidate stellar eclipse around phase 1, and also deep dips during the egress, which are often seen around phase 1.1–1.5. We used only observations with a complete phase coverage between 0.85 and 1.65. In terms of orbital phase, the duration of the occultations is Δϕ_occ = N_occ × 500/P_best ≈ 0.0193N_occ. We calculated the uncertainty of the half-width as a Poisson error on the number of phase bins (i.e., $0.0193\sqrt{{N}_{\mathrm{occ}}}$).

Finally, to test whether the duration of the occultations is correlated with the peak luminosity, we used the online Portable, Interactive Multi-Mission Simulator¹³ (pimms) version 4.9 to convert all of the count rates from different epochs and instruments to the effective count rate and unabsorbed luminosity for Cycle 3 Chandra ACIS-S3 (no grating) in 0.3–8 keV band. The spectral model that we used for the conversion is a power law with photon index Γ = 1.7 and absorbing column density N_H = 1 × 10²² cm⁻².

We found a correlation (Figure 5) between the duration of the occultations and the peak luminosity, which is represented by the normalized maximum count rate CR_max. At low luminosities (L_X ≲ 10⁴⁰ erg s⁻¹), the duration of the occultations is linearly anticorrelated with the count rate, as Δϕ_occ = (−1.96 ± 0.97) × CR_max + (0.70 ± 0.13). To quantify the statistical significance of the correlation, we calculated the Spearman correlation coefficient ρ: we obtained ρ ≈ −0.73, with a p value of 5 × 10⁻⁶. At higher luminosities (L_X > 10⁴⁰ erg s⁻¹), the duration of the occultations saturates around Δϕ_occ ≈ 0.2–0.3, independent of luminosity. The precise slope of the correlation depends on how the peak count rate and the occultation bins are defined, but the existence of a significant trend is a robust result.

Zoom In Zoom Out Reset image size

The next question that we addressed is whether and how the spectrum of the observed photons changes between fainter and brighter sections of an orbital cycle, and between orbital cycles over the years. We did this in the following three steps.

First (Section 4.2), we split a selected sample of high signal-to-noise light curves (from Chandra/ACIS observations) into phenomenological substructures (ingress, eclipse, dips, egress, bright phase) and compared the cumulative energy distribution of the observed photons in the various substructures. The objective of this part of our analysis is to detect spectral changes in a model-independent way.

Second, we did a phase-resolved spectral modeling of the four highest-quality Chandra/ACIS-S3 observations (Section 4.3) and all five XMM-Newton/EPIC observations (Section 4.4). The objective of this part of our analysis is to model how the fit parameters change between the brighter and fainter sections of an orbital cycle. Thus, we split each of the selected observations into three phase groups: a bright phase, an intermediate phase, and a faint phase, drawing on the results of our light-curve analysis. We defined the bright phase as the time intervals when the count rates are higher than 70% of the maximum count rate for that orbital cycle; we defined the faint phase as the intervals when the count rates are less than 15% of the maximum count rate for XMM-Newton, and 10% for Chandra; and we defined the intermediate phase as the time bins in between the faint and the bright phase. To increase the number of counts in each of the three subintervals, we combined their spectra from the four Chandra datasets. Instead, we had enough photons to analyze the five XMM-Newton observations individually. For XMM-Newton, we combined pn and MOS spectra together with the sas task epicspeccombine to increase the signal-to-noise ratio.

Finally, in the third step of our analysis (Section 4.5), we extracted and modeled spectra integrated over entire orbital cycles, from six Chandra observations and five XMM-Newton ones, between 2000 and 2018 (the datasets that we used for this modeling are marked by asterisks in Table 1). We determined average and peak fluxes and luminosities during those observations.

The four Chandra observations that we chose for our model-independent study of the photon energy distribution are ObsIDs 9140, 10937, 12823 and 12824, covering a total of about 11 orbital cycles.¹⁴ The advantage of these particular observations is that it is relatively straightforward to identify the following five main substructures in the background-subtracted light curves: faint phase, bright phase, ingress, egress and dip (these five substructures are color-coded in Figure 6). More specifically, for this part of the analysis we defined a maximum count rate as the average value of the 10 500 s bins with the highest count rate during each observation. We then defined the faint phase as the time intervals when the count rates are lower than 10% of the maximum count rate. Next, we defined the bright phase as the time intervals when the count rates are higher than 70% of the maximum count rate. Finally, the ingress and egress phases are defined as the intervals of decreasing and increasing count rates, respectively, between the bright and faint phases. The definition of a dip is somewhat more arbitrary but corresponds to an approximate flux drop of at least a factor of 2 followed by an immediate recovery within ≲2000 s, which are color-coded as red data points in Figure 6. To extract the observed net counts from the dips, we considered the time of the local minimum count rate as the mid-time of a dip and we took the photons falling within ±200 s of that mid-time in the unbinned event file.

Zoom In Zoom Out Reset image size

We performed the Kolmogorov–Smirnov (K-S) statistical test for the null hypothesis that the cumulative photon energy distributions of faint, bright, ingress, egress and dip phases are the same. We find (Figure 6, bottom panel) that photons from the last four (bright, ingress, egress and dips) of those five structures do indeed follow the same energy distribution. Instead, photons from the faint intervals are significantly softer. The K-S test rejects the null hypothesis that the faint intervals follow the same distribution as the other intervals, with a p value of 7 × 10⁻⁸. This suggests the presence of at least two emission components: a bright and harder component, and a faint and softer component that is seen as a residual component when the other component is occulted.

We used the same four Chandra observations selected in Section 4.2 (ObsIDs 9140, 10937, 12823 and 12824). We extracted three phase-resolved spectra, averaged over the four observations but distinguished by count-rate brackets; we refer to these as the bright-phase, intermediate-phase and faint-phase spectra (or, more simply, the bright, intermediate and faint spectra). First, we fitted the three phase-resolved spectra independently, in the 0.5–7 keV range, using standard one-component models suitable to X-ray binaries: power-law, bremsstrahlung, diskbb, diskpbb, and bbodyrad. The emission components were convolved with two neutral absorbers (tbabs model): one for the Galactic line of sight absorption (column density N_H fixed at 5.6 × 10²¹ cm⁻², from Kalberla et al. 2005) and one left free for the intrinsic absorption inside the Circinus galaxy and the binary system. Three of those models (power-law, bremsstrahlung and diskpbb with p < 0.75) give good fits for all three phases (Table 2): the diskbb model is significantly less good, both for the bright phase and for the faint phase; the bbodyrad model is the worst, especially for the faint phase, which has an unacceptable ${\chi }_{\nu }^{2}$ = 2.25. The phenomenological interpretation is that both the bright and the intermediate spectra have a low degree of curvature in the Chandra bandpass, so they are best represented either by a power law or by the (relatively flat) low-frequency section of a thermal continuum component (e.g., bremsstrahlung with a temperature >8 keV, or p-free disk with peak temperature >2.5 keV). The intermediate-phase spectrum has a harder (flatter) slope than the bright spectrum, at least in the ≈1–5 keV range, but a lower intrinsic N_H. Conversely, the faint spectrum (dominated by residual emission in the eclipse) is significantly softer (steeper).

Table 2.  Best-fitting Parameters of Phase-resolved Spectra from Combined Chandra ObsIDs 9140, 10937, 12823, and 12824

Phase	NH,intr	Γ/kT (keV)	p	NH,absori	FX	LX	${\chi }_{\nu }^{2}$ (dof)
	(1022 cm−2)			(1022 cm−2)	(10−13 erg cm−2 s−1)	(1038 erg s−1)
tbabs1 × tbabs2 × po
Bright	${0.50}_{-0.05}^{+0.06}$	${1.76}_{-0.06}^{+0.06}$	⋯	⋯	${17.1}_{-0.4}^{+0.4}$	${66.5}_{-2.7}^{+3.1}$	0.87 (267)
Intermediate	${0.14}_{-0.10}^{+0.11}$	${1.27}_{-0.13}^{+0.14}$	⋯	⋯	${6.4}_{-0.4}^{+0.4}$	${18.2}_{-1.0}^{+1.1}$	0.85 (100)
Faint	<0.25	${2.42}_{-0.29}^{+0.38}$	⋯	⋯	${0.57}_{-0.07}^{+0.06}$	${3.3}_{-0.6}^{+1.7}$	0.76 (30)
tbabs1 × tbabs2 × bremss
Bright	${0.38}_{-0.04}^{+0.04}$	${8.39}_{-1.09}^{+1.36}$	⋯	⋯	${16.7}_{-0.4}^{+0.4}$	${56.5}_{-0.3}^{+1.4}$	0.85 (267)
Intermediate	${0.13}_{-0.04}^{+0.09}$	$\gt 24.3$	⋯	⋯	${5.8}_{-0.3}^{+0.3}$	${17.4}_{-0.9}^{+1.0}$	0.84 (100)
Faint	<0.82	${2.34}_{-0.53}^{+0.45}$	⋯	⋯	${0.52}_{-0.06}^{+0.04}$	${2.27}_{-0.16}^{+0.17}$	0.91 (30)
tbabs1 × tbabs2 × diskbb
Bright	${0.19}_{-0.03}^{+0.03}$	${1.75}_{-0.08}^{+0.08}$	⋯	⋯	${15.8}_{-0.4}^{+0.4}$	${45.9}_{-1.0}^{+1.0}$	0.95 (267)
Intermediate	<0.05	${2.50}_{-0.03}^{+0.04}$	⋯	⋯	${6.0}_{-0.4}^{+0.4}$	${15.6}_{-0.8}^{+0.8}$	0.82 (100)
Faint	<0.82	${0.80}_{-0.09}^{+0.10}$	⋯	⋯	${0.46}_{-0.05}^{+0.05}$	${1.76}_{-0.13}^{+0.13}$	1.36 (30)
tbabs1 × tbabs2 × diskpbb
Bright	${0.49}_{-0.05}^{+0.06}$	${7.1}_{-1.2}^{+* }$	${0.54}_{-0.01}^{+0.01}$	⋯	${17.0}_{-0.4}^{+0.4}$	${60.1}_{-3.0}^{+9.2}$	0.86 (266)
Intermediate	<0.17	${2.68}_{-0.68}^{+6.77}$	${0.73}_{-0.12}^{+0.06}$	⋯	${6.0}_{-0.5}^{+0.5}$	${15.8}_{-1.2}^{+2.8}$	0.83 (99)
Faint	<0.06	${1.47}_{-0.37}^{+0.69}$	${0.50}_{-0.05}^{+0.03}$	⋯	${0.53}_{-0.7}^{+0.08}$	${2.53}_{-0.28}^{+0.19}$	0.86 (29)
tbabs1 × absori × tbabs2 × po
Bright	${0.54}_{-0.08}^{+0.08}$	${2.03}_{-0.14}^{+0.15}$	⋯	${1.75}_{-0.89}^{+1.03}$	${16.8}_{-0.4}^{+0.4}$	${90.1}_{-13.1}^{+17.8}$	0.83 (265)
Intermediate	${0.31}_{-0.18}^{+0.21}$	${1.89}_{-0.43}^{+0.41}$	⋯	${3.8}_{-2.7}^{+3.1}$	${6.1}_{-0.4}^{+0.4}$	${30.4}_{-9.7}^{+19.2}$	0.80 (98)
Faint	${0.12}_{-0.12}^{+0.09}$	${2.63}_{-0.31}^{+0.25}$	⋯	${2.3}_{-2.3}^{+4.1}$	${0.58}_{-0.06}^{+0.07}$	${5.3}_{-2.4}^{+1.1}$	0.73 (28)
tbabs1 × tbabs2 × bbodyrad
Bright	<0.01	${0.87}_{-0.01}^{+0.01}$	⋯	⋯	${14.0}_{-0.3}^{+0.3}$	${35.7}_{-0.7}^{+0.7}$	1.76 (267)
Intermediate	<0.01	${0.99}_{-0.04}^{+0.04}$	⋯	⋯	${5.0}_{-0.3}^{+0.3}$	${12.4}_{-0.6}^{+0.7}$	1.46 (100)
Faint	<0.01	${0.49}_{-0.05}^{+0.05}$	⋯	⋯	${0.40}_{-0.04}^{+0.05}$	${1.4}_{-0.1}^{+0.1}$	2.25 (30)

Note. tbabs₁ is the line of sight absorption in the direction of Circinus, fixed at N_H = 0.56 × 10²² cm⁻²; tbabs₂ is the intrinsic neutral absorption; absori is the ionized absorption. F_X is the absorbed flux in the 0.3–8 keV band; L_X ≡ 4πd²F_X is the emitted 0.3–8 keV luminosity. See Section 4.1 for the definition of "bright," "intermediate," and "faint" phase intervals.

Download table as:  ASCIITypeset image

Instead of a simple power law, we also tried Comptonization models (in particular, comptt, simpl × diskbb, and diskir). However, these models add a layer of complexity and additional free parameters without providing any improvement to the fit. Below 1 keV, the thermal seed component of the Comptonization models is unconstrained because of moderately high absorption; the possible high-energy downturn above ∼5 keV (typical of ULXs, Stobbart et al. 2006; Gladstone et al. 2009; Sutton et al. 2013; Walton et al. 2018) is also unconstrained because of low signal-to-noise in the Chandra spectra. Thus, we stick to the simple power-law model for the Chandra analysis.

Before we can attempt a more physical interpretation of the spectra, consistent with the proposed super-Eddington HMXB scenario, we need to account for one additional source of absorption, which comes from ionized gas. We do this by adding an absori component. Here, we take the power-law model as the benchmark to determine whether the addition of an ionized absorber improves the fit. We find that an ionized absorber with column density N_H ∼ a few 10²² cm⁻² and an ionization parameter ∼a few 100 does provide a significant improvement for both the bright and the intermediate spectra (Table 2). For the bright spectrum, the goodness-of-fit improves from ${\chi }_{\nu }^{2}\,=232.7/267$ to ${\chi }_{\nu }^{2}=218.7/265$: this is significant to >99.9% probability.¹⁵ For the intermediate spectrum, the improvement is from ${\chi }_{\nu }^{2}$ = 85.0/100 to ${\chi }_{\nu }^{2}=78.2/98$, significant at the >95% level. A similar amount of ionized absorption is also consistent with the faint spectrum; however, in that case, because of the lower signal-to-noise level, the ionized absorber only improves the fit at the 1σ level (${N}_{{\rm{H}}}={2.3}_{-2.3}^{+4.1}\,\times {10}^{22}$ cm⁻²).

Moreover, when the ionized absorber is included, we find that the intrinsic power-law slope of the bright and intermediate spectra is consistent with being the same (Γ ≈ 2.0); the reason why the intermediate spectrum looks flatter is because of a factor-of-two higher column density of the ionized absorber. The faint spectrum remains significantly steeper (Γ ≈ 2.6) than the other two (Table 2).

We can now start to identify some of the physical properties of the spectral evolution. The main reason for the increase in observed count rates from the intermediate to the bright spectrum is neither a dramatic decrease in absorption column density nor a state transition in the intrinsic emission properties. While there are changes in the fitted column densities of ionized and neutral components, they only affect the shape of the spectrum below 2 keV. Instead, in a first approximation, the main difference between bright and intermediate spectra is consistent with a reduced normalization of the broadband emission component during the egress and dipping phases. Considering that this evolution happens regularly during each orbital cycle, we consider it unlikely that it is due to intrinsic changes in the source emission. A much less contrived explanation (also, by analogy, with other dipping X-ray binaries) is that the flux changes are due to variable partial covering of the emission region by clumps of optically thick material.

We tested this physical interpretation with a new spectral model. We fitted the bright and intermediate spectra simultaneously, with a tbabs × tbabs × absori × pcfabs × power-law model, keeping both the slope and the normalization of the power-law component locked for the two spectra, while allowing column densities and ionization parameters to vary independently. The partial-covering absorber modeled with pcfabs has to be Compton thick, with N_H ≳ 2 × 10²⁴ cm⁻¹, so that it blocks at least 99% of the flux below 7 keV. We also assumed a covering fraction f = 0 for the spectrum in the bright phase. We find that this model provides the best fit (${\chi }_{\nu }^{2}$ = 297.7/363 = 0.82). The photon index of the intrinsic emission component is Γ ≈ 2.1 and the partial-covering fraction of the intermediate spectrum is f ≈ 0.60. The average intrinsic (de-absorbed) luminosity during the four Chandra observations that we used in our spectral modeling is L_X ≈ 9 × 10³⁹ erg s⁻¹, in the 0.3–8 keV band. (Higher luminosities in excess of 10⁴⁰ erg s⁻¹ were found during some of the XMM-Newton observations.)

We also recovered the result that the faint spectrum (residual emission in eclipse) is softer and less absorbed than the emission out of eclipse, which cannot be explained with the same Γ ≈ 2 model used for the intermediate and bright spectra. The simplest way to model this spectrum is to assume that the dominant emission component seen in the bright and intermediate spectra is completely blocked by the opaque screen (f ≡ 1 in pcfabs) and there is instead a different emission component that is only significantly detectable in eclipse. We also model this component with a simple power law, for want of higher signal-to-noise spectra. We find a photon index Γ ≈ 2.6 and a luminosity L_X ≈ 5 × 10³⁸ erg s⁻¹ in the 0.3–8.0 keV band (Table 3).

Table 3.  Best-fitting Parameters for Chandra and XMM-Newton Phase-resolved Spectra, for Our Fiducial Model

Phase	NH	ξ	f	NH	Γ	Norm	FX	LX	${\chi }_{\nu }^{2}$ (dof)
	(1022 cm−2)	(102)		(1022 cm−2)		(10−4)	(10−13 CGS)	(1038 CGS)
	absori	absori	pcfabs	tbabs2	po	po	cflux
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
Chandra combined ObsIDs 9140, 10937, 12823 and 12824
Fainta	${2.34}_{-2.34}^{+4.13}$	unconstr.	⋯	${0.12}_{-0.12}^{+0.09}$	${2.63}_{-0.31}^{+0.25}$	${0.50}_{-0.03}^{+0.07}$	${0.58}_{-0.06}^{+0.07}$	${5.3}_{-2.4}^{+1.1}$	0.73 (28)
Brightb	${1.96}_{-0.91}^{+1.01}$	${4.9}_{-2.4}^{+3.8}$	$[0]$	${0.57}_{-0.8}^{+0.8}$	${2.10}_{-0.16}^{+0.16}$	${9.04}_{-1.94}^{+2.48}$	${16.5}_{-0.5}^{+0.5}$	${89}_{-12}^{+16}$	0.82 (363)
Interm.b	${5.11}_{-1.68}^{+2.02}$	${6.0}_{-2.5}^{+4.3}$	${0.60}_{-0.04}^{+0.04}$	${0.40}_{-0.11}^{+0.11}$			${5.9}_{-0.3}^{+0.3}$
XMM-Newton obs. 0111240101 (1)
Brightb	${2.04}_{-0.59}^{+0.73}$	${1.7}_{-0.9}^{+1.5}$	[0]	${0.38}_{-0.16}^{+0.13}$	${1.97}_{-0.10}^{+0.11}$	${22.72}_{-3.45}^{+4.35}$	${48.5}_{-1.0}^{+1.0}$	${255}_{-27}^{+36}$	1.00 (192)
Interm.b	${3.38}_{-0.72}^{+0.84}$	${1.4}_{-0.4}^{+1.4}$	${0.42}_{-0.03}^{+0.03}$	${0.16}_{-0.13}^{+0.28}$			${26.5}_{-0.8}^{+0.8}$
XMM-Newton obs. 070198100 (2)
Brightb	${4.31}_{-2.15}^{+1.41}$	unconstr.	[0]	${0.71}_{-0.27}^{+0.11}$	${2.44}_{-0.31}^{+0.16}$	${14.03}_{-2.15}^{+1.41}$	${15.7}_{-2.4}^{+1.1}$	${131}_{-45}^{+38}$	1.17 (134)
Interm.b	${2.39}_{-1.42}^{+0.92}$	${8.4}_{-0.8}^{+30.6}$	${0.38}_{-0.11}^{+0.07}$	${0.47}_{-0.25}^{+0.15}$			${10.0}_{-2.3}^{+2.5}$
XMM-Newton obs. 0656580601 (3)
Brightb	${1.65}_{-1.54}^{+2.49}$	${2.4}_{-2.4}^{+5.5}$	[0]	<0.59	${2.04}_{-0.36}^{+0.37}$	${7.60}_{-3.52}^{+5.89}$	17.2 ${}_{-2.2}^{+2.4}$	${71.1}_{-10.6}^{+34.4}$	0.95 (113)
Interm.b	${2.07}_{-2.07}^{+1.38}$	unconstr.	${0.43}_{-0.11}^{+0.09}$	<1.0			9.1 ${}_{-1.5}^{+6.2}$
XMM-Newton obs. 0792382701 (4)
Brightb	${3.21}_{-0.88}^{+2.49}$	unconstr.	[0]	${0.59}_{-0.14}^{+0.15}$	${2.10}_{-0.18}^{+0.12}$	${10.81}_{-1.62}^{+3.35}$	21.7 ${}_{-1.2}^{+0.6}$	${112}_{-12}^{+27}$	1.12 (98)
Interm.b	${1.29}_{-1.29}^{+3.04}$	${24.6}_{-2.1}^{+1.4}$	${0.25}_{-0.15}^{+0.13}$	<0.98			${14.8}_{-1.1}^{+1.6}$
XMM-Newton obs. 0780950201 (5)
Brightb	${0.41}_{-0.41}^{+2.08}$	<9.2	[0]	${0.29}_{-0.24}^{+0.16}$	${1.99}_{-0.22}^{+0.28}$	${4.85}_{-0.09}^{+1.83}$	${11.9}_{-0.9}^{+0.9}$	${54.8}_{-8.0}^{+24.4}$	0.93 (110)
Interm.b	${12.38}_{-12.37}^{+11.41}$	unconstr.	${0.45}_{-0.22}^{+0.22}$	${0.34}_{-0.22}^{+0.22}$			${5.59}_{-0.58}^{+0.76}$

Notes. Columns: (1) Identification of the orbital phase corresponding to the best-fitting spectral parameter listed on each row (see Section 4 for the definition of Faint, Intermediate and Bright phases); (2) equivalent hydrogen column density of the ionized absorber (modeled with absori); (3) absori ionization parameter ξ = L/(N_eR²), where L is the integrated luminosity between 5 eV and 300 keV, N_e is the electron number density, and R is the distance from the source to the ionized material (Done et al. 1992); (4) dimensionless covering fraction of the Compton-thick absorber modeled with pcfabs (0 ≤ f ≤ 1); (5) Equivalent hydrogen column density of the intrinsic neutral absorber (in the host galaxy and around CG X-1); (6) photon index of the power-law component, assumed identical for Bright and Intermediate phases; (7) power-law normalization, in units of 10⁴ photons keV⁻¹ cm⁻² s⁻¹ at 1 keV; (8) absorbed flux in the 0.3–8 keV band; (9) unabsorbed luminosity in 0.3–8 keV band (L_X ≡ 4πd²F_X, with d = 4.2 Mpc); (10) reduced χ² and degrees of freedom. The models used for the phase-resolved spectra are as follows:

^atbabs₁ × absori × tbabs₂ × po. ^btbabs₁ × absori × pcfabs × tbabs₂ × po. Other parameters not shown in this table are fixed at default or assumed values. Specifically: (i) for absori, the temperature, Fe abundance and redshift are fixed at 3 × 10⁴ K, 1.00, and 0, respectively; (ii) for pcfabs, the column density of the Compton-thick absorber is fixed at N_H = 2 × 10²⁴ cm⁻², but identical fitting results are obtained for any other higher values; (iii) the line of sight absorption was modeled with an additional tbabs₁ component, with column density fixed at N_H = 0.56 × 10²² cm⁻².

Download table as:  ASCIITypeset image

It is plausible (in the framework of our interpretation of the system) that part or all of the softer component seen in the eclipse spectra is also present in the intermediate and bright spectra. To check this, we refitted all three spectra simultaneously, including the eclipse-phase component as an additional (fixed) component in the intermediate and bright spectra. However, this does not change their best-fitting parameters because the residual component is lost in the noise of the much brighter out-of-eclipse spectra. Therefore, for simplicity, we ignore this additional term when fitting the bright and intermediate spectra, for Chandra and, in the next section, for XMM-Newton.

For the XMM-Newton data, we only fitted the spectra of the bright and intermediate phases because the faint-phase emission is highly contaminated by the diffuse emission in the inner region of the Circinus galaxy and by the PSF wings from the active galactic nucleus. We also ignored all data above 6 keV because the strong and inhomogeneous background emission (including Fe lines from the nuclear source) always dominates over the CG X-1 emission in that band.

We used the same combination of models discussed for the Chandra data: tbabs × absori × pcfabs × tbabs × po. The first tbabs component represents the line of sight absorption and is fixed at N_H = 5.6 × 10²¹ cm⁻²; the second tbabs and the absori components represent Compton-thin neutral and ionized absorption around CG X-1, respectively; and the pcfabs component is a partial-covering, Compton-thick absorber responsible for the dips and occultations. As before, we assumed that in the bright phase, the covering fraction f of pcfabs is f = 0. We also assumed that the photon index and normalization of the power-law emission is the same in the bright and intermediate spectra, so that the difference between the two phases is entirely due to changes in the column densities of the Compton-thin absorbers and to an increased covering fraction of the Compton-thick medium. Although the best-fitting values change from observation to observation without a clear trend, most spectra are well fitted with intrinsic cold-absorber column densities ∼2–5 × 10²¹ cm⁻² and ionized-absorber column densities ∼2–4 × 10²² cm⁻² (Table 3), which is consistent with the properties of the Chandra spectra. The covering fraction for the intermediate spectra varies between ∼0.3 and 0.6. This is not surprising because the intermediate spectra were extracted (by definition) from time bins with count rates between 15% and 70% of the maximum count rate.

In summary, as for the Chandra spectra, the main difference between the bright and intermediate spectra from each individual observation is the energy-independent covering fraction of the Compton-thick screen (Figure 7). The photon index of the power-law emission is consistent with Γ ≈ 2.1 within the 90% confidence limit of every observation (which is also consistent with the Chandra result). The intrinsic luminosity of the bright phase is itself variable from orbit to orbit, ranging from ≈5.5 × 10³⁹ erg s⁻¹ to ≈2.6 × 10⁴⁰ erg s⁻¹, in the 0.3–8 keV band.

Zoom In Zoom Out Reset image size

Finally, we used the bright-phase XMM-Newton spectra to test for the presence of a high-energy downturn, which is one of the defining properties of ULXs (Stobbart et al. 2006; Gladstone et al. 2009; Sutton et al. 2013; Walton et al. 2018). To do this, we fitted the XMM-Newton spectra above 2 keV with a simple power law and a broken power law (bknpow in xspec), independently. We restricted our fit to energies >2 keV because the effect of neutral and ionized absorption is statistically negligible in this region. We find that the broken power-law model is statistically preferred over the simple power law, at the 98% significance level, according to the F-test (Table 4). The best-fitting break energy is (4.4 ± 0.3) keV, and the continuum slope steepens from Γ ≈ 2.0 to Γ ≈ 2.8 (Figure 8).

Zoom In Zoom Out Reset image size

Table 4.  Comparison of Power-law and Broken Power-law Models in the 2–6 keV Band (XMM-Newton Observation 1)

Γ	χ2/dof	Γ1	Ebreak	Γ2	χ2/dof	1 − P
po	po	bknpo	bknpo	bknpo	bknpo
(1)	(2)	(3)	(4)	(5)	(6)	(7)
${2.36}_{-0.16}^{+0.17}$	73.82/72	${1.97}_{-0.33}^{+0.29}$	${4.44}_{-0.53}^{+0.27}$	${2.78}_{-0.37}^{+0.43}$	66.41/70	97.5

Note. Columns: (1) Photon index of the single power-law model; (2) goodness-of-fit for the single power-law model; (3) photon index below the break for the bknpo model; (4) break energy (keV); (5) photon index above the break; (6) goodness-of-fit for the bknpo model; and (7) percentage probability that the broken power-law model is a better fit than the single power-law model, from the F-test.

Download table as:  ASCIITypeset image

After analyzing the differences between brighter and fainter phases within individual orbital cycles (Sections 4.3 and 4.4), we set out to estimate the range of luminosity variability for the bright phase from cycle to cycle, over the past two decades. We have already shown a range of luminosities based on the maximum observed count rates converted through pimms (Section 3.3 and Figure 5), now we want to estimate peak luminosities more accurately from spectral analysis. To do this, we used six Chandra and five XMM-Newton observations from 2000 to 2018, as already mentioned in Section 4.1. However, this time we extracted and fitted the average spectrum of each observation (i.e., no longer split into faint, intermediate and bright intervals), while MOS and pn data were fitted simultaneously for each XMM-Newton observation. We applied our fiducial model with intrinsic ionized and neutral absorption, and foreground neutral absorption: tbabs × absori × tbabs × powerlaw. We did not include a Compton-thick pcfabs term because we assume that it gives a gray correction to the spectrum, equivalent to a rescaling of its normalization. For each observation, we determined the average de-absorbed luminosity in the 0.3–8 keV band (Table 5). We also measured the average count rate and the maximum count rate (see Section 4.2 for the definition of maximum count rate). We then multiplied the average luminosity of each observation by the ratio of maximum over average count rate for that same observation. This gave us the peak luminosity of each observation.

Table 5.  Best-fitting Parameters of the Phase-averaged Spectra from a Sample of Individual Observations

Obs.	NH	NH	Γ	Norm	$\left\langle {F}_{{\rm{X}}}\right\rangle$	$\left\langle {L}_{{\rm{X}}}\right\rangle$	LX,peak	${\chi }_{\nu }^{2}$ (dof)
	(1022 cm−2)	(1022 cm−2)		(10−4)	(10−13 erg cm−2 s−1)	(1038 erg s−1)	(1038 erg s−1)
	absori	tbabs	po	po
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
365	${0.48}_{-0.48}^{+3.16}$	${0.16}_{-0.16}^{+0.51}$	${1.51}_{-0.20}^{+0.22}$	${10.32}_{-2.49}^{+3.79}$	${46.2}_{-3.0}^{+3.0}$	${158}_{-16}^{+25}$	${158}_{-16}^{+25}$	0.89 (85)
356a	${1.94}_{-1.10}^{+1.55}$	${0.56}_{-0.30}^{+0.19}$	${2.09}_{-0.52}^{+0.28}$	${14.98}_{-6.20}^{+9.03}$	${29.0}_{-2.3}^{+1.8}$	${160}_{-38}^{+76}$	${205}_{-49}^{+91}$	0.87 (166)
1	${2.43}_{-0.90}^{+1.36}$	${0.71}_{-0.16}^{+0.16}$	${1.99}_{-0.12}^{+0.13}$	${16.15}_{-2.83}^{+4.02}$	${33.2}_{-0.3}^{+0.6}$	${179}_{-22}^{+34}$	${288}_{-35}^{+55}$	1.00 (752)
9140	${6.47}_{-4.59}^{+2.13}$	${0.21}_{-0.17}^{+0.09}$	${1.79}_{-0.15}^{+0.21}$	${1.12}_{-0.47}^{+0.69}$	${2.98}_{-0.16}^{+0.32}$	${13.9}_{-2.6}^{+5.7}$	${35.2}_{-9.5}^{+14.7}$	1.02 (63)
10937	${2.12}_{-2.12}^{+7.04}$	${0.20}_{-0.20}^{+0.50}$	${2.25}_{-0.70}^{+0.66}$	${2.62}_{-1.62}^{+3.84}$	${4.2}_{-0.5}^{+0.6}$	${26.8}_{-12.8}^{+39.5}$	${53.2}_{-25.3}^{+78.2}$	0.88 (22)
12823	${1.57}_{-0.89}^{+1.05}$	${0.49}_{-0.07}^{+0.07}$	${2.02}_{-0.14}^{+0.14}$	${4.18}_{-0.83}^{+1.09}$	${8.82}_{-0.21}^{+0.22}$	${46.1}_{-6.7}^{+9.1}$	${107}_{-16}^{+21}$	0.97 (303)
12824	${2.25}_{-1.70}^{+2.01}$	${0.50}_{-0.16}^{+0.15}$	${2.16}_{-0.25}^{+0.29}$	${4.19}_{-1.55}^{+2.03}$	${7.06}_{-0.36}^{+0.37}$	${43.8}_{-12.6}^{+18.3}$	${87.8}_{-25.2}^{+36.6}$	0.87 (120)
2	${1.90}_{-1.72}^{+2.28}$	${0.61}_{-0.54}^{+0.07}$	${2.31}_{-0.31}^{+0.17}$	${6.07}_{-2.33}^{+4.21}$	${9.60}_{-0.47}^{+0.47}$	${61.8}_{-20.1}^{+36.4}$	${127}_{-41}^{+75}$	0.98 (300)
3	${1.97}_{-1.74}^{+1.93}$	${0.50}_{-0.50}^{+0.21}$	${2.15}_{-0.34}^{+0.27}$	${5.22}_{-2.01}^{+2.56}$	${10.01}_{-0.57}^{+0.60}$	${54.9}_{-16.9}^{+21.8}$	${84.0}_{-25.9}^{+33.3}$	0.94 (212)
4	${1.33}_{-1.33}^{+1.47}$	${0.59}_{-0.59}^{+0.22}$	${2.15}_{-0.31}^{+0.26}$	${8.42}_{-3.19}^{+5.34}$	${15.33}_{-0.79}^{+0.95}$	${88.6}_{-25.4}^{+45.5}$	${132}_{-38}^{+68}$	0.88 (182)
5	${2.60}_{-2.60}^{+4.03}$	${0.58}_{-0.22}^{+0.42}$	${2.28}_{-0.34}^{+0.25}$	${3.21}_{-0.34}^{+0.25}$	${5.34}_{-0.43}^{+0.46}$	${29.6}_{-7.6}^{+16.4}$	${56.6}_{-14.5}^{+31.5}$	1.01 (152)

Note.

^aChandra ObsID 356 suffers from pileup at the 40% level (ACIS-S3 count rate ≈0.2 ct s⁻¹). We fitted its spectrum with pileup × tbabs × absori × tbabs × po. The model used for all the other observations is tbabs × absori × tbabs × po. $\left\langle {F}_{{\rm{X}}}\right\rangle$ is the average observed flux of each observation and $\left\langle {L}_{{\rm{X}}}\right\rangle$ the average unabsorbed luminosity ($\left\langle {L}_{{\rm{X}}}\right\rangle \equiv 4\pi {d}^{2}\left\langle {F}_{{\rm{X}}}\right\rangle$), which are both in the 0.3–8 keV band. L_X,peak is the inferred luminosity of the unocculted portions of each observation (see Section 4.5).

Download table as:  ASCIITypeset image

The highest peak luminosity was seen in the first XMM-Newton observation (2001 August 6), at L_X,peak ≈ 2.9 × 10⁴⁰ erg s⁻¹. The lowest luminosity was recorded in Chandra ObsID 9140 (2008 October 26), with L_X,peak ≈ 3.5 × 10³⁹ erg s⁻¹ (Table 5). Thus, L_X,peak varies by a factor of ≈8 across our subsample of 11 observations (Figure 9), whereas the average luminosity of each observation varies by a factor of ≈13. The variability range of the average luminosities is higher because observations with lower peaks also have longer occultation phases (Figure 5) and therefore even lower average luminosities. There is no trend in the hardness-luminosity plane; that is, the best-fitting photon index Γ is approximately the same for all observations. For 10 out of the 11 observations in our spectral analysis subsample, i.e., all except Chandra ObsID 365, the mean photon index is $\left\langle {\rm{\Gamma }}\right\rangle \approx 2.11$ with a scatter σ ≈ 0.15. Only the spectrum of Chandra ObsID 365 is significantly different from the others, with Γ = 1.5 ± 0.2.

Zoom In Zoom Out Reset image size

A faint, point-like optical source was found at the X-ray position (Bauer et al. 2001; Weisskopf et al. 2004), in a combined 600 s exposure from HST-WFPC2 with the F606W filter (Figure 10). We reprocessed the data to check or improve the astrometric alignment between X-rays and optical bands, and to verify whether that source is indeed the most likely counterpart of CG X-1. We aligned the Chandra and HST positions using five bright, point-like sources detected in both images. The residual random scatter between the best-fitting positions of the reference sources in the two bands gives us an estimate of the error radius for the position of CG X-1. We are able to narrow down the error circle to ≈02 (white circle in Figure 10), which improves on the previous results. The error circle does indeed include the optical source that was previously suggested to be the most likely counterpart. Applying standard techniques of HST/WFPC2 aperture photometry (Sirianni et al. 2005) to the source, and using the most updated values of the WFC2 zeropoints,¹⁶ we confirm an apparent brightness m_F606W ≈ V ≈ (24.3 ± 0.1) mag in the Vega system, which is in agreement with Weisskopf et al. (2004).

Zoom In Zoom Out Reset image size

Taking into account the line of sight extinction A_F606W ≈ 3.5 mag (Schlafly & Finkbeiner 2011), and a distance modulus of 28.1 mag for the Circinus galaxy, we estimate the intrinsic brightness of the source as M_V ≲ −7.3 mag (depending on the amount of additional local extinction). This is already at the high end of the luminosity range for a WR star (M_V ranging approximately from −7 to −2.5 mag: Massey 2003; Crowther 2007). Given that CG X-1 is located in a star-forming region of Circinus, near dust lanes and filaments, we also expect a non-negligible local extinction, which is likely very high and uncertain. Thus, it is more likely that the optical counterpart is an unresolved young star cluster. The field of CG X-1 was also observed with HST/WFPC2 in the narrow-band Hα filter F656N and O III filter F502N, but no excess line emission is detected in these bands. In conclusion, with only one image in one optical filter, we do not have enough information to constrain the color, spectral type, or luminosity of the optical counterpart.

We have already mentioned (Section 1) that the identification of CG X-1 as a possible foreground Galactic source (in particular, an mCV projected by chance in front of Circinus) has been a point of contention in previous studies (Weisskopf et al. 2004), but was subsequently rejected by Esposito et al. (2015). Here, we confirm those arguments and suggest additional ones against the mCV interpretation, based on our X-ray and optical results.

The X-ray over optical flux ratio F_X/F_opt is a useful criterion to identify and classify CVs because it removes the uncertainty on the source distances. Although the precise definition of the energy bands used for the X-ray and optical fluxes may change from author to author, the general conclusions remain the same (within a factor of two). The ratio F_X(0.1–4.0 keV)/F_opt(5000–6000 Å) shows an upper boundary at around 4 (Patterson & Raymond 1985). Using the standard V band, and converting it to a flux with the relation¹⁷ log F_opt = −0.4V − 5.44, we infer F_X(2–10 keV)/F_opt ≲ 7 (Figure 3 in Mukai 2017), for both magnetic and non-magnetic CVs.¹⁸ A similar analysis by Revnivtsev et al. (2014) shows that F_X(0.5–10 keV)/F_opt ≲5 for non-magnetic CV (for this relation, we have used the conversion νF_ν ≈ 5F_λ × 1000 Å ≡ F_opt at λ ≈ 5000 Å). All of these empirical relations are based on observed (rather than de-absorbed) fluxes, but most of the sources are within a few 100 pc and their line of sight absorbing column density is negligible. For CG X-1, even if it were a CV in the Milky Way, its foreground line of sight optical extinction and X-ray absorption would be significant because of its location in the direction of the Galactic plane. Given the uncertainty on its true distance, for CG X-1 we have calculated both the observed and the de-absorbed F_X/F_opt ratios, and compared them with the upper limits that we found in the CV surveys that we cited earlier. We have not removed the additional intrinsic absorption component derived from X-ray fitting because that component is also not removed in those CV surveys.

First, we compare the observed fluxes. From the net count rate in the WFPC2 F606W band, and the tabulated value of photflam,¹⁹ we obtain F_λ ≈ 5.2 × 10⁻¹⁹ erg cm⁻² s⁻¹ Å⁻¹, and (by our definition, assuming a flat spectrum) F_opt ≈ 5.2 × 10⁻¹⁶ erg cm⁻² s⁻¹. For the X-ray flux, we take the unocculted phases of Chandra/ACIS ObsID 12823 (the longest in our series, Table 1 and Figure 6). The observed 2–10 keV flux is F_X ≈ 1.9 × 10⁻¹² erg cm⁻² s⁻¹, giving an implausibly high X-ray over optical flux ratio of ≈3650. For a more physical result, we remove the line of sight extinction (A_V ≈ 4 mag) from the optical flux and the line of sight absorption (corresponding to N_H ≈ 6 × 10²¹ cm⁻²) from the X-ray flux. We obtain F_opt ≈ 2.06 × 10⁻¹⁴ erg cm⁻² s⁻¹ and F_X ≈ 2.03 × 10⁻¹² erg cm⁻² s⁻¹, corresponding to F_X(2–10 keV)/F_opt ≈ 100. These values are well outside the range observed in CVs, which are then definitively ruled out. Instead, F_X/F_opt ∼ 100–1000 is what we expect and observe in typical ULXs (Tao et al. 2011; Gladstone et al. 2013; Ambrosi & Zampieri 2018) because of an X-ray luminosity ≈10⁴⁰ erg s⁻¹ (from the compact object) and an optical luminosity of ≈10³⁸ erg s⁻¹ (combination of the contributions from the irradiated disk and the massive donor star).

Finally, as an additional test of the foreground source scenario, we point out that the 6.4 keV Fe line is usually strong in mCVs, and their absence is very rare (Butters et al. 2011). To check if there is a significant 6.4 keV line emission in CG X-1, we extracted Chandra/ACIS X-ray images in the 6.4–6.7 keV band and then combined all the images from different observations. We found no concentration of photons in this band at the position of CG X-1. The few detected photons are clearly distributed as diffuse background. Also, no significant 6.4–6.7 keV lines were detected in stacked background-subtracted X-ray spectra (Section 4.4).

In conclusion, we suggest that our new arguments about the flux ratio and the lack of Fe lines, together with those of Esposito et al. (2015), put the final nail in the coffin of the Galactic mCV scenario.

The binary separation in CG X-1 is

$$\begin{eqnarray}\begin{array}{rcl}a & = & {\left[\displaystyle \frac{{{GP}}^{2}({M}_{1}+{M}_{2})}{4{\pi }^{2}}\right]}^{1/3}\\ & \simeq & \ 5.8\,{\left(\displaystyle \frac{P}{7.2\mathrm{hr}}\right)}^{2/3}\,{\left(\displaystyle \frac{{M}_{1}+{M}_{2}}{30{M}_{\odot }}\right)}^{1/3}\,{R}_{\odot },\end{array}\end{eqnarray} \tag{ 5 }$$

where G is the gravitational constant, P the orbital period in unit of hr, M₁ and M₂ the masses of the compact object and donor star, respectively, in unit of M_⊙. The radius of the secondary Roche lobe is

$$\begin{eqnarray}&&{R}_{{\rm{L}},2}=a\,\displaystyle \frac{0.49{q}^{2/3}}{0.6{q}^{2/3}+\mathrm{ln}(1+{q}^{1/3})}\end{eqnarray} \tag{ 6 }$$

(Eggleton 1983), where q ≡ M₂/M₁ is the mass ratio. The average mass density inside the secondary Roche lobe (lower limit to the average density of the donor star) is

$$\begin{eqnarray}&&\rho =0.2101\displaystyle \frac{q}{1+q}{\left(\displaystyle \frac{{R}_{{\rm{L}},2}}{a}\right)}^{-3}{\left(\displaystyle \frac{P}{7.2\mathrm{hr}}\right)}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-3}\end{eqnarray} \tag{ 7 }$$

(Eggleton 1983).

In principle, the duration of the eclipse provides constraints on the mass ratio (coupled with the inclination angle), which was already well discussed by Weisskopf et al. (2004). However, the problem for CG X-1 is that we are no longer sure what the duration of the true eclipse is (given the large cycle-to-cycle variability) and what part of the occultation is, due to other thick structures between donor star and accretor. So, at this point we want to use only the most general constraint from Equation (7). For any plausible value of q suitable to stellar-mass binaries (0.01 ≲ q ≲ 100), the average density in the donor Roche lobe must be 0.55 ≲ ρ g cm⁻³ ≲ 2.2 (0.4 ≲ ρ/ρ_⊙ ≲ 1.6). For a more likely range of 0.1 ≲ q ≲ 10, the density is 1 ≲ ρ g cm⁻³ ≲ 2.2(0.7 ≲ ρ/ρ_⊙ ≲ 1.6). This density range rules out massive early-type main-sequence stars, supergiants, red giants, asymptotic giants, and white dwarfs as the donor star. Instead, low-mass main-sequence stars, slightly evolved stars (such as low-mass helium stars) and massive WR stars are consistent with this tight constraint.

For a peak L_X ≈ 3 × 10⁴⁰ erg s⁻¹, the mass accretion rate must be a few 10⁻⁶ M_⊙ yr⁻¹, even for radiatively efficient accretion. Considering also that in the supercritical regime, radiative efficiency decreases because of advection and outflows (e.g., Poutanen et al. 2007), we suggest that the mass transfer rate into the Roche lobe of the accretor approaches 10⁻⁵ M_⊙ yr⁻¹. However, can a low-mass donor (with the constraints mentioned in Section 5.3) provide such high mass transfer rate? A small number of ULXs have been found in elliptical galaxies (David et al. 2005; Plotkin et al. 2014; van Haaften et al. 2019) and old globular clusters (Maccarone et al. 2007; Shih et al. 2010; Roberts et al. 2012; Dage et al. 2018). Meanwhile, white dwarf donors in ultracompact systems are a plausible scenario for the ULX population in globular clusters (Steele et al. 2014; Dage et al. 2018). However, ULXs in old stellar populations are much rarer than in young stellar environments (an order of magnitude rarer above ≈5 × 10³⁹ erg s⁻¹: Swartz et al. 2004; Plotkin et al. 2014). Moreover, we have already ruled out the possibility of a white dwarf donor and we have also ruled out a globular cluster as the optical counterpart of CG X-1. Other ULXs with a low- or intermediate-mass donor have been seen in intermediate-age environments of star-forming galaxies (Soria et al. 2012). Population synthesis models show the possibility of ULXs with a NS accretor and a low-mass companion stars, where the donor may have started its life as an ≈6 M_⊙ (Wiktorowicz et al. 2015). Those models show that after a common-envelope phase, the companion star is stripped of its hydrogen envelope, the binary separation is reduced, and the system goes through a ULX phase with a helium-star donor mass of ≈1–2 M_⊙ (either a Helium Hertzsprung gap or a Helium giant branch star).

While we cannot dismiss the possibility of a stripped, low-mass helium-star, we suggest that the most likely type of donor is a WR star, which is in agreement with Esposito et al. (2015). We briefly summarize a few arguments that are consistent with this scenario. First, a WR donor can provide a sufficiently high mass transfer rate onto a BH accretor, producing persistent X-ray luminosities L_X ∼ 10³⁹–10⁴¹ erg s⁻¹ (Bogomazov 2014; Wiktorowicz et al. 2015), for timescales of a few 10⁵ yr. Second, WR stars have typical radii of only ≈1–2 × 10¹¹ cm (≈1.5–3 R_⊙) (Crowther 2007) and can fit in a binary orbit with a period of a few hours. For a M₁ = 10 M_⊙ BH and a M₂ = 20 M_⊙ donor, and a period of 7.2 hr, R_L,2 ≈ 2.6 R_⊙ (large enough to fit a WR star), and ρ ≈ 1.6 g cm⁻³ ≈ 1.2 ρ_⊙ (also consistent with a typical WR). Third, the short orbital period and eclipsing/dipping X-ray light-curve profile of CG X-1 resemble the observed properties of other X-ray binaries that have been interpreted as WR systems (e.g., Bauer & Brandt 2004; Ghosh et al. 2006; Carpano et al. 2007; Zdziarski et al. 2012; Esposito et al. 2013; Maccarone et al. 2014; Laycock et al. 2015a). The reason for this interpretation is that a WR provides a dense wind that strongly affects the X-ray light curve via photoelectric absorption, especially when the compact object transits behind the donor star. Stronger stellar winds and smaller binary separations make this effect more important in WR HMXBs than, for example, in those with a supergiant donor. In Sections 6.1 and 6.2, we will discuss in more detail how the dense WR wind may explain the peculiar X-ray light curve of CG X-1. We will revisit the comparison with other candidate WR X-ray binaries in Section 6.5.

We have shown (Section 3.2 and Figure 3) that the folded X-ray light curve (averaged over dozens of cycles) has an apparent eclipse (lasting for ≈1/4 of the period), with a sharp ingress and a slow egress. We have also shown (Figure 2) that this simple picture is complicated by irregular dips in each single cycle, and that the X-ray profile and the duration of the full occultation differ substantially from cycle to cycle. From spectral analysis, we have shown (Sections 4.2 and 4.3) that the transition from lower to higher fluxes is energy-independent; hence, it is best explained as a decreasing level of partial covering by a totally opaque medium—rather than, for example, by a gradual decrease of the column density of the photoelectric absorber. Any model of the system need to explain these unusual X-ray properties.

The sharp ingress and regular occurrence of the occultation around phase 1 suggest that there is a proper eclipse by the companion star, but the irregular duration of this phase also suggests a contribution from other optically thick material, such as clouds or clumps of dense gas (See Figure 11 for a cartoon of the binary geometry). Irregular dips are well-known in several LMXBs seen at moderately high inclination (White & Swank 1982; Frank et al. 1987; Díaz Trigo et al. 2006; D'Aì et al. 2014). In these systems, the occulting material is either the accretion stream or the geometrically thick bulge where the stream impacts the outer rim of the accretion disk. However, because of conservation of angular momentum along its ballistic trajectory, the stream is always trailing the compact object along the orbit. Therefore, dips caused by the accretion stream and impact bulge happen before eclipse ingress, which is contrary to what we see in CG X-1.

Zoom In Zoom Out Reset image size

Asymmetric eclipses are seen in some HMXBs (Falanga et al. 2015), such as the NS system Vela X-1. In this case, the asymmetric absorption has been attributed to a stream-like region of slower and denser wind trailing the NS (Doroshenko et al. 2013). Photoelectric absorption gradually increases during the slow ingress and reaches the maximum at superior conjunction. Again, this type of enhanced wind absorption is not consistent with the energy-independent dips of CG X-1 and with their preferential location after the eclipse. However, energy-independent occultations were detected ("Type B dips": Feng & Cui 2002) in another well-known HMXB, the BH system Cyg X-1. A possible explanation for these dips is partial covering of an extended X-ray emitting region by an opaque screen (Feng & Cui 2002).

WR X-ray binaries are a subclass of HMXBs with a thicker stellar wind. Typical mass-loss rates are $\dot{M}$ ∼ 10⁻⁴–10⁻⁵ M_⊙ yr⁻¹, with terminal velocities of ∼1000–3000 km s⁻¹ (Crowther 2007; Gräfener et al. 2017). Thus, we expect that WR X-ray binaries are most likely to show the effect of variable wind absorption on X-ray light curves, especially just before and just after superior conjunction of the compact object. Indeed, this is what is seen in the BH–WR systems NGC 300 X-1 (Carpano et al. 2007, 2018; Crowther et al. 2010; Binder et al. 2015) and IC 10 X-1 (Silverman & Filippenko 2008; Barnard et al. 2014; Laycock et al. 2015b; Steiner et al. 2016). In particular, the X-ray light curve of NGC 300 X-1 also shows dips during a slow egress, which is consistent with variable partial covering of an extended X-ray emitting region (Binder et al. 2015). However, unlike CG X-1, the variable absorbing clumps are not Compton-thick (column densities of only ≈10²³ cm⁻²), thus they have an energy-dependent effect on the 0.3–10 keV spectrum. A Compton-thick absorber was found in IC 10 X-1, from NuSTAR observations (Steiner et al. 2016), when the BH passes behind the WR star. The X-ray light curve of IC 10 X-1 also has the same type of asymmetric profile as CG X-1, with a steep ingress and a slow egress; however, it has been suggested that in IC10 X-1, the X-ray source is never completely occulted by the star (Barnard et al. 2014), and even its eclipse (more exactly, a dip in flux by a factor of 10) is caused by the thick wind.

Based on the analogy with the wind absorption in WR HMXBs, in the next section we will try to assess where the Compton-thick clumps could be located in CG X-1 to explain the asymmetric light curve and dips. We also keep in mind another crucial piece of evidence: CG X-1 is a super-Eddington source, while NGC 300 X-1 and IC 10 X-1 are sub-Eddington BHs. For example, the peak luminosity of CG X-1 (≈3 × 10⁴⁰ erg s⁻¹) is 300 times higher than the luminosity of IC 10 X-1 (≈7 × 10³⁷ erg s⁻¹, Laycock et al. 2015a). Therefore, the compact object in CG X-1 is expected to launch a much stronger wind (King & Pounds 2003; Ohsuga et al. 2005; Pinto et al. 2016, 2017; Walton et al. 2016; Kosec et al. 2018) than in the other systems.

There are two obvious regions where we expect a density enhancement (over the already high value inside a WR wind) with possible Compton-thick clump formation. The first location is along the line between WR and compact accretor, which is caused by wind–wind collision (WR wind and super-Eddington wind from the accretion disk). The second location is the bow shock in front of the compact object, which is caused by its fast orbital motion inside a dense medium. For simplicity, we shall also assume that the compact object is a stellar-mass BH (Figure 11).

Let us start from the wind collision. Shocks produced by colliding winds have been extensively studied (especially in the context of WR-O star binaries), both theoretically (e.g., Stevens et al. 1992; Usov 1992; Kenny & Taylor 2005; Parkin & Pittard 2008; Lamberts et al. 2012) and observationally (e.g., Dougherty et al. 2005; Pollock et al. 2005; Zhekov 2012; Hill et al. 2018; Nazé et al. 2018). Revisiting the detailed properties of the shocked gas is well beyond the scope of this paper: here we just want to check (with order-of-magnitude arguments) whether the shocked gas layer can act as a Compton-thick absorber. Let us assume for simplicity that both the ULX disk wind (component 1) and the WR wind (component 2) have a similar projected speed in the orbital plane, v_w1 ∼ v_w2 ∼ 1000 km s⁻¹. For the ULX outer disk wind, this is comparable to the escape velocity from the outermost disk annulus. From the observed X-ray luminosity, we inferred an average mass accretion rate ≳10⁻⁶ M_⊙ yr⁻¹. At such super-Eddington rates, the amount of mass lost in disk outflows is predicted to be larger than the mass accreted through the BH event horizon. Therefore, we estimate a mass-loss rate of a few times 10⁻⁶ M_⊙ yr⁻¹, in the super-Eddington disk outflow; this is consistent with our expectations that ∼10% of the mass lost by the WR is accreted by the BH.²⁰ The intersection of the contact discontinuity between the two colliding winds with the line of centers will be located at a distance d₁ from the compact object and d₂ from the donor star, with d₁/d₂ = (${\dot{M}}_{1}/{\dot{M}}_{2}{)}^{1/2}$ (v_w1/v_w2)^1/2 and d₁ + d₂ = a (Stevens et al. 1992). For example, for ${\dot{M}}_{1}\approx 0.1{\dot{M}}_{2}$, d₂ ≈ (3/4)a. The density ρ_w2 of a spherically symmetric wind from the WR star at a distance d₂ from its center is

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{{\rm{w}}2} & = & \displaystyle \frac{{\dot{M}}_{2}}{4\pi {d}_{2}^{2}\,{v}_{{\rm{w}}2}}\approx 6.4\times {10}^{-12}\left(\displaystyle \frac{\dot{{M}_{2}}}{{10}^{-5}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)\\ & & \times \ {\left(\displaystyle \frac{{d}_{2}}{4{R}_{\odot }}\right)}^{-2}{\left(\displaystyle \frac{{v}_{{\rm{w}}2}}{{10}^{3}\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-1}{\rm{g}}\,{\mathrm{cm}}^{-3}.\end{array}\end{eqnarray} \tag{ 8 }$$

From standard shock theory, the density of the shock-heated WR wind is a factor of 4 higher, which corresponds to a number density n_s2 ∼ 10¹³ cm⁻³. The temperature of the shock-heated gas is kT_s2 ≈ (3/16) μ m_H v_w2² ≈ 2 keV for the assumed wind speed of ≈1000 km s⁻¹, where m_H ≈ 1.7 × 10⁻²⁴ g is the mass of the hydrogen atom, and μ is the mean mass per particle of gas measured in unit of m_H. On the other side of the contact discontinuity, there will be another layer of shocked gas from the disk wind, with a similar temperature and density n_s1 ≈ (v_w2/v_w1)² n_s2 ∼ 10¹³ cm⁻³.

The cooling timescale for the shocked WR wind is

$$\begin{eqnarray}&&{t}_{\mathrm{cool},2}=\displaystyle \frac{3{\left({n}_{e}+n\right)}_{{\rm{s}}2}\,{{kT}}_{{\rm{s}}2}}{2{n}_{{\rm{s}}2}^{2}\,{\rm{\Lambda }}({T}_{{\rm{s}}2})}\approx \displaystyle \frac{3\,{{kT}}_{{\rm{s}}2}}{{n}_{{\rm{s}}2}\,{\rm{\Lambda }}({T}_{{\rm{s}}2})}\end{eqnarray} \tag{ 9 }$$

(Dopita & Sutherland 2003), where Λ is the cooling function (Sutherland & Dopita 1993). An analogous expression holds for the cooling timescale t_cool,1 of the ULX wind. Thus, for a WR mass-loss rate of a few 10⁻⁵ M_⊙ yr⁻¹, we expect a cooling timescale of a few seconds for the shocked wind, and a few 10 s of seconds for the shocked ULX wind, in a system with a characteristic size of CG X-1. The cooling timescales can be compared with the timescales for the hot shocked gas on either side of the contact discontinuity to leave the interaction region; that is, the escape timescales t_esc,1 = d₁/c_s1 and t_esc,2 = d₂/c_s2, where c_s1 and c_s2 are the sound speeds in the two shocked layers. If χ ≡ t_cool/t_esc < 1, then a shocked shell is radiative; otherwise, it is adiabatic (Stevens et al. 1992). For the physical parameters assumed for our toy model of CG X-1, χ₂ ∼ 10⁻³ for the shocked WR wind and χ₁ ∼ 10⁻² for the shocked ULX wind. Thus, both shells are radiative and will collapse into geometrically thin, cold shells.²¹

Next, we want to estimate the surface density of the radiative shells, to estimate whether they can provide the column density required to absorb essentially all X-ray emission below 10 keV. For this purpose, we adapt the analytical solution of Kenny & Taylor (2005). In addition to our previous assumptions on wind speed and mass-loss rate, we assumed initial wind temperatures of a few 10⁴ K, both for the WR and the ULX wind, and a temperature of ≈10⁴ K for the radiatively cooled shocked shells. The density in a radiative shell is enhanced by a compression factor approximately equal to ${{ \mathcal M }}^{2}$, where ${ \mathcal M }$ is the Mach number (Weaver et al. 1977). For typical sound speeds ∼10 km s⁻¹ and wind speeds ∼1000 km s⁻¹, the compression factor is ∼10⁴. Hence, we expect the two layers of radiatively cooled gas on both sides of the contact discontinuity to have densities of a few 10⁻⁸ g cm⁻³. From the solutions of Kenny & Taylor (2005), inserting typical parameters for this system, we obtain that each of the two shells has a thickness δ ∼ 10⁻⁴ times the binary separation, that is a total thickness δ₁+δ₂ ∼ 10⁸ cm. Thus, the surface density of the cold screen is ∼a few g cm⁻², which corresponds to a Compton-thick column density of a few 10²⁴ cm⁻²: this is enough to block all X-ray photons in the Chandra and XMM-Newton bands.

Crucially, the opaque screen is subject to a dynamical process of fragmentation, known as thin-shell instability (Stevens et al. 1992). We expect the cold gas to get shredded into clumps and filaments, with characteristic sizes ∼10⁸ cm. The X-ray emitting region in a ULX is approximately the region of the geometrically thick accretion flow inside the spherization radius, which is a function of the mass accretion rate. A plausible value of the spherization radius for the observed X-ray luminosity of CG X-1 is a few 1000 km = a few 10⁸ cm, which is comparable with the size of the occulting clumps. Thus, we conclude that the shocked wind between the two system components is dense enough to produce Compton-thick absorption and it can fragment into clumps with the right size to produce partial covering of the X-ray emission. If the obscuring clumps were much smaller in size, then we would not see sharp, individual dips in the light curve, but only smooth flux changes. Conversely, if the clumps were much bigger than the X-ray emitting region, then we would either see total eclipses of the full flux, rather than variable partial covering. An analogous scenario of partial covering of the X-ray source by Compton-thick clouds was used to explain occultations in the Seyfert galaxy NGC 1365 (Risaliti et al. 2009a, 2009b).

As mentioned earlier, the wind collision region is not the only place where a thin shell of shocked gas can be formed. Under the plausible approximation of a circular orbit, the Keplerian orbital velocity v_K1 of the compact object in CG X-1 is

$$\begin{eqnarray}&&{v}_{{\rm{K}}1}={\left[\displaystyle \frac{{{GM}}_{2}^{2}}{a\left({M}_{1}+{M}_{2}\right)}\right]}^{1/2}.\end{eqnarray} \tag{ 10 }$$

For example, for M₁ = 10 M_⊙ and M₂ = 20 M_⊙, we obtain v_K1 ≈ 660 km s⁻¹; for M₁ = 5 M_⊙ and M₂ = 30 M_⊙, v_K1 ≈ 900 km s⁻¹; for M₁ = M₂ = 20 M_⊙, v_K1 ≈ 550 km s⁻¹. An object with a strong wind, moving into dense ambient medium at such high speed, will necessarily produce a strong bow shock. The physics of bow shocks in front of fast-moving OB stars is also a well-studied problem (e.g., Wilkin 1996; Comeron & Kaper 1998; Meyer et al. 2017), based on the same physical elements as the colliding wind problem (forward shock into the WR wind; contact discontinuity; reverse shock into the ULX wind). Near the apex of the bow shock, the WR wind moves perpendicularly to the shock but the forward shock is driven by the supersonic orbital motion of the BH, which, in the case of CG X-1, is only slightly lower than the WR wind speed. The distance between the BH and the contact discontinuity at the apex is

$$\begin{eqnarray}&&{R}_{0}={\left(\displaystyle \frac{{\dot{M}}_{1}{v}_{{\rm{w}}1}}{4\pi {\rho }_{{\rm{a}}}{v}_{{\rm{K}}1}^{2}}\right)}^{1/2}\end{eqnarray} \tag{ 11 }$$

(Comeron & Kaper 1998), where ρ_a is the density of the unshocked WR wind near the apex, at a distance of ≈a from the stellar center. Using Equation (8), we can recast this distance as

$$\begin{eqnarray}&&{R}_{0}={\left(\displaystyle \frac{{\dot{M}}_{1}}{{\dot{M}}_{2}}\right)}^{1/2}{\left(\displaystyle \frac{{v}_{{\rm{w}}1}{v}_{{\rm{w}}2}}{{v}_{{\rm{K}}1}^{2}}\right)}^{1/2}\,a.\end{eqnarray} \tag{ 12 }$$

For our assumptions of ${\dot{M}}_{1}\sim 0.1{\dot{M}}_{2}$, and wind speeds of ≈1000 km s⁻¹, this implies R₀ ≈ 0.5a ≈ 3 R_⊙. In terms of angular distance, the apex of the bow shock is ≈30° in front of the BH.

In the case of colliding winds, it is easy to verify that both the shocked ambient gas and the shocked ULX wind are fully radiative and collapse to two thin, cold shells on either side of the contact discontinuity. The surface density σ of the cold gas can be approximated (Comeron & Kaper 1998) as σ ≈ 1.25ρ_aR₀ ∼ 1 g cm⁻² for ${\dot{M}}_{2}\sim {10}^{-5}{M}_{\odot }$ yr⁻¹. This corresponds to an absorbing column density ∼10²⁴ cm⁻², which is again sufficient to block most of the X-ray emission below 10 keV. The bow shock will be seen directly in front of the BH at the orbital phase ϕ ≈ 1.2. The combined effect of true stellar eclipse, colliding winds and bow shock may keep the X-ray emission continuously occulted from phase ϕ ≈ 0.90 to ϕ ≈ 1.1, and partially covered from phase ϕ ≈ 1.1 to ϕ ≈ 1.5, as can be seen in some orbital cycles (Figure 2). Instead, we expect to see unocculted emission roughly from ϕ ≈ 1.5 to ϕ ≈ 1.8 when the bow shock is behind the X-ray source (Figures 3 and 11). For the colliding wind shell, the cold bow-shock shell is also subject to instabilities and fragmentation (Blondin & Koerwer 1998), which cause irregular dips during each cycle, variable partial covering, and profile variability from cycle to cycle (including cycles where the X-ray emission is already close to its maximum value around phase 1.2).

In addition to the slower equatorial disk winds that are considered in our model, ULXs also launch relativistic or semi-relativistic outflows at lower inclination angles from the polar axis (Pinto et al. 2016, 2017; Walton et al. 2016; Kosec et al. 2018). These outflows may create larger, parsec-scale shells of shocked circumstellar gas around the binary system. However, these shells will be Compton-thin and they are not directly relevant to the phase-dependent occulting behavior that we are studying in this work; therefore, a study of the circumbinary ULX bubble is beyond the scope of this paper.

We have shown (Figure 5) a trend of faster egress times in cycles with higher X-ray luminosities in the bright phase. We have also argued from spectral analysis (Section 4) that the transition between bright and occulted phases of a cycle is due to an increase of partial covering by an opaque screen. Thus, we suggest that the egress width is correlated with the area and location of the opaque screen with respect to the X-ray emission region: higher X-ray luminosities seem to reduce the phase interval affected by partial covering.

One simple explanation for the anti-correlation is that at higher (intrinsic) luminosities, larger parts of the opaque screen are ionized by the strong X-ray flux and become optically thin to X-rays; this leads to shorter occultations and faster egress times. An alternative explanation is that the higher luminosity also corresponds to stronger outflows from the compact object. This may push the opaque screen and the bow shock to larger distances from the X-ray source (as discussed in Section 6.2), or even blow it away, thus reducing the fraction of the X-ray emission region covered along our line of sight. A third, more speculative, scenario is that changes in the observed X-ray luminosity are due to disk precession: when the X-ray-emitting inner disk is seen more face-on, we expect to see a higher X-ray luminosity in the unocculted phases, and a smaller degree of occultation by opaque clouds in the orbital plane (Hu et al. 2013). However, further investigations on these possibilities are beyond the scope of this paper.

It is significant that the duration of the full or partial occultation phase (as defined in Section 3.3) saturates at a minimum value Δϕ ≈ 0.2 at X-ray luminosities ≳10⁴⁰ erg s⁻¹ (Figure 5). We speculate that this minimum value is evidence of a true eclipse by the donor star that can never be eliminated, even when the partial-covering clouds are ionized or blown away. However, even during the full occultation/eclipse phase, there is still residual emission, which we will discuss in the next section.

One of the results of our spectral analysis of phase-resolved stacked spectra is that even during the true eclipse (or deepest occultation) at ϕ ≈ 0.9–1.1, we still detect residual flux from CG X-1, at a level of ≈1% of the bright-phase flux (L_X ≈ 2.5 × 10³⁸ erg s⁻¹). The residual flux is softer and less absorbed than the flux in the bright phase (Section 4). Chandra's image resolution shows that the residual emission is associated with the point-like source, and is not, for example, due to incomplete background subtraction of diffuse X-ray emission in the inner region of the Circinus galaxy. Here, we assess possible origins for this residual emission.

The residual emission cannot be from the WR star itself: the X-ray luminosity of an isolated WR star does not exceed a few times 10³² erg s⁻¹ (Skinner et al. 2012). Wind–wind collisions might offer a more promising explanation, analogous to the X-ray emission seen from WR + O star binaries (e.g., Guerrero & Chu 2008). Expected wind speeds ≳1000 km s⁻¹ for both the WR wind and the super-Eddington disk wind guarantee that the shocked gas reaches temperatures ≳1 keV. In the simplest case of identical winds, each with mass-loss rate $\dot{M}$ and speed v_w, in the fully radiative regime, the predicted luminosity is ${L}_{{\rm{X}}}\approx f\dot{M}{v}_{{\rm{w}}}^{2}\approx {10}^{36}{\dot{M}}_{-5}\,{v}_{{\rm{w}},3}^{2}$ erg s⁻¹ (Stevens et al. 1992), where the geometrical factor f ≈ 1/6 parameterizes the fraction of wind kinetic luminosity perpendicular to the wind shock, ${\dot{M}}_{-5}$ is in units of 10⁻⁵ M_⊙ yr⁻¹, and v_w,3 is in units of 1000 km s⁻¹. For unequal winds, in the isothermal regime, the ratio of the luminosities from the two shocked winds is ${L}_{1}/{L}_{2}\approx {({\dot{M}}_{1}/{\dot{M}}_{2})}^{1/2}{({v}_{{\rm{w}}1}/{v}_{{\rm{w}}2})}^{3/2}$ (Stevens et al. 1992). In the case of CG X-1, we can plausibly assume ${\dot{M}}_{1}\lesssim 0.1{\dot{M}}_{2}$, as already discussed in Section 6.2. For the X-ray luminosity of the shocked winds to exceed 10³⁸ erg s⁻¹, values of ${\dot{M}}_{1}\gtrsim {10}^{-6}$ M_⊙ yr⁻¹, and v_w1 ∼ 0.1c would be required. The luminosity would be completely dominated by the disk wind component. We cannot rule out this scenario for a ULX because we know that the kinetic power carried by the outflows is comparable to the radiative power (Siwek et al. 2017). So far, we have not identified any other ULXs where such kinetic power is dissipated as X-rays in a strong shock at small distances from the compact object, as opposed to large ionized bubbles (ULX bubbles: Pakull & Mirioni 2002; Pakull et al. 2006) at ∼100 pc scale. Furthermore, before CG X-1, we have not identified any WR ULX, immersed in a circumstellar medium with densities ∼10¹² cm⁻³. Consequently, the theoretical properties of a compact WR ULX bubble are an intriguing topic for further work.

Another scenario is that the residual soft X-ray emission comes from the photosphere of the optically thick outflows driven by the super-Eddington source. This thermal component has been invoked to explain the spectra of ultraluminous supersoft sources and the "soft excess" in ULX spectra (Middleton et al. 2015; Soria & Kong 2016; Urquhart & Soria 2016; Tao et al. 2019). This is generally consistent with a blackbody spectrum at temperatures ≈50–120 eV for supersoft sources, and ≈0.1–0.4 keV for ULXs, with a soft X-ray luminosity from ∼10³⁸ erg s⁻¹ up to a few 10³⁹ erg s⁻¹ (Urquhart & Soria 2016; Zhou et al. 2019). If we fit the residual emission of CG X-1 with a blackbody model, we obtain a temperature of ≈0.5 keV (higher than seen in the soft excess of any other ULX) and an X-ray luminosity of ≈10³⁸ erg s⁻¹ (Table 2). In the outflow model of Zhou et al. (2019), this solution corresponds to a NS accreting at almost 100 times Eddington. The physical size of the outflow photosphere (r_ph) is related to the apparent blackbody radius (r_bb) obtained from spectral fitting as ${r}_{\mathrm{ph}}\approx \sqrt{{\tau }_{\mathrm{es}}}\,{r}_{\mathrm{bb}}$ (Meier 1982; Zhou et al. 2019), where τ_es is the scattering optical depth at the photosphere (τ_es ∼ a few 100). The photons emitted from the photosphere will be scattered multiple times before they can propagate freely to us from the last scattering sphere, where τ_es decreases to unity. The radius of this surface may be ∼100 times larger than the radius of the photosphere. In our case, r_bb ≈ 130 km, plausible values for r_ph are a few 10⁸ cm, and the radius of the last scattering sphere may be ∼10¹¹ cm, which is comparable to the size of the occulting star and the bow shock. We note that a blackbody model is a poor fit to our eclipse spectrum, with ${\chi }_{\nu }^{2}$ > 2 (Table 2). However, the spectrum emerging from the last scattering sphere may be significantly different from the input blackbody emitted at the photosphere of the optically thick outflows.

Alternatively, we suggest that the most likely source for the residual eclipse emission is that a fraction of the X-ray photons are emitted along the polar funnel and are then downscattered by the wind into our line of sight, at scales larger than the companion star. Residual emission in eclipse is seen in various other Galactic HMXBs, and is usually interpreted as scattered light (Sako et al. 1999; Schulz et al. 2002; Wojdowski et al. 2003; Lopez et al. 2006). This is usually a few percent of the direct luminosity out of eclipse, which is consistent with our flux ratio. A WR system such as CG X-1, with strong outflows from both binary components, is precisely the type of system where we expect that some radiation emitted perpendicular to the binary plane will be scattered into our line of sight, even when the direct emission is obscured.

To-date, we know of only seven compact (i.e., with a binary period ≲1 day) HMXBs with a candidate WR donor: in addition to CG X-1, they are: NGC 4214 X-1 (P = 3.6 hr, Ghosh et al. 2006); Cyg X-3 in the Milky Way (P = 4.8 hr, Zdziarski et al. 2012); NGC 4490 CXOU J123030.3+413853 (X-1) (P = 6.4 hr, Esposito et al. 2013); NGC 253 CXOU J004732.0-251722.1 (P = 14.5 hr Maccarone et al. 2014); NGC 300 X-1 (P = 32.8 hr, Carpano et al. 2007); IC 10 X-1 (P = 34.8 hr, Laycock et al. 2015a). One more candidate WR X-ray binary is M 101 ULX-1 (Liu et al. 2013), but it has a tentative period of 8.2 d, which suggests a different evolutionary history than the compact class. Within this group, CG X-1 is the most luminous (peak X-ray luminosity L_X ≈ 3 × 10⁴⁰ erg s⁻¹) and is the only ULX, clearly exceeding the Eddington limit for a 10M_⊙ BH.

We focus on compact WR binaries because they help us to understand two key evolutionary features of binary systems, one in their past history and one in their future. In the past, they must have come through a phase of dramatic shrinking of the binary separation from initial values of ∼100–1000 R_⊙ to ≲10 R_⊙. The most likely mechanism for such shrinking is a common-envelope phase after the formation of the compact object (e.g., Dominik et al. 2012; Bogomazov 2014; Belczynski et al. 2016; van den Heuvel et al. 2017; Bogomazov et al. 2018; Giacobbo et al. 2018). In some cases, theoretical calculations suggest (Pavlovskii et al. 2017; van den Heuvel et al. 2017; Bogomazov et al. 2018) that a substantial but gradual shrinking of a compact binary system (a "spiral-in" during stable Roche-lobe overflow mass transfer) may occur without a common-envelope phase, if the donor star has a radiative envelope and the mass ratio is ≲3; however, a common-envelope phase seems unavoidable (van den Heuvel et al. 2017) if we want to obtain binary periods as short as a few hours (like in CG X-1 and Cyg X-3). In the future, if they remain bound and their orbit does not widen too much after core collapse of the WR star, then systems like CG X-1 and Cyg X-3 will form double compact objects (DCOs) that can spiral in and merge via gravitational wave emission on timescales much shorter than the Hubble time (e.g., Bulik et al. 2011; Belczynski et al. 2013; Esposito et al. 2015; Belczynski et al. 2016). However, other types of DCOs originating for example from a compact object orbiting a supergiant donor star (i.e., evolved from supergiant HMXBs) are too widely separated to do so. Thus, population studies of short-period WR X-ray binaries, and especially ULX WR systems (detectable with X-ray telescopes in a larger volume of space), will provide crucial observational constraints to the predicted rate of LIGO/Virgo detections (Abbott et al. 2016a, 2016b, 2017a, 2017b, 2017c). An analysis of the physical scenarios and possible steps between the X-ray binary stage and the final merger is, however, beyond the scope of this work. Instead, we will briefly discuss the past and present evolutionary stages of these systems.

The common-envelope phase that is required for the existence of systems such as CG X-1 is the one that occurs after the formation of the first compact object, which results in the loss of the hydrogen-rich envelope of the companion star and the dramatic shrinking of the binary separation. After the end of common envelope, the donor star is helium-rich and can be classified as a WR star.²² From this point, further evolution of the system (which is now an X-ray binary) depends on the mass ratio, on how mass and angular momentum are transferred from the donor to the compact object, and on how they are lost from the system (Huang 1963; Lommen et al. 2005; van den Heuvel et al. 2017). For example, mass transfer from a more massive WR donor to a less massive compact component and outflows from the accretion disk wind are two factors that lead to orbital shrinking, while mass loss in the wind from the more massive component results in orbital widening. In particular, for WR donors with wind mass-loss rates of ∼10⁻⁵ M_⊙ yr⁻¹, we expect a secular increase of binary separation and orbital period (Tutukov et al. 2013). This is already observed in Cyg X-3 ($\dot{P}/P$ ≈ 1.2–4 × 10⁻⁶ yr⁻¹: Lommen et al. 2005). In this work (Section 3.2), we have shown a period increase at the 90% confidence limit in CG X-1 ($\dot{P}/P=(10.2\pm 4.6)\times {10}^{-7}$ yr⁻¹).

Orbital widening inevitably leads to the donor star becoming detached from its Roche lobe because WR stars do not expand in their final stages of evolution (unlike supergiants). For a fixed mass-loss rate ${\dot{M}}_{2}$ in the WR wind, the mass transfer rate ${\dot{M}}_{1}$ into the Roche lobe of the compact object scales as ${\dot{M}}_{1}/-{\dot{M}}_{2}\approx {G}^{2}{M}_{1}^{2}/({a}^{2}\,{v}_{{\rm{w}}2}^{4})$ (Frank et al. 2002). As the WR donor becomes more detached, we expect a decrease in the mass accretion rate and luminosity, both because of the increase in the binary separation a, and because the WR wind reaches a higher velocity v_w2 before being intercepted by the compact object; the latter effect makes it harder to form an accretion disk around the accretor (Ergma & Yungelson 1998).

We have already noted that CG X-1 is the only ULX in the known sample of seven close WR X-ray binaries. The higher luminosity of CG X-1 compared with NGC 300 X-1 and IC 10 X-1 can be explained by a larger orbital separation in these two systems (i.e., a period longer than a day). However, Cyg X-3 has an even shorter period than CG X-1, and yet its X-ray luminosity is almost two orders of magnitude lower, L_X ∼ a few 10³⁸ erg s⁻¹ (Hjalmarsdotter et al. 2008; Koljonen et al. 2010, 2018). Clearly, the binary period alone cannot be used to predict the X-ray luminosity of compact WR binaries; the mass, evolutionary stage, metal abundance of the WR star certainly plays an equally important role. Moreover, the accretion cross-section scales as M₁²; the mass of the compact object is not known, either in Cyg X-3 or in CG X-1. For example, for Cyg X-3, mass estimates for the compact object based on accretion-state behavior range from ≈2 M_⊙ (Zdziarski et al. 2013), which is consistent with either a BH or a NS, to a run-of-the-mill stellar-mass BH around 10 M_⊙ (Shrader et al. 2010), to a heavy BH with a mass ≳20 M_⊙ (Hjalmarsdotter et al. 2008).

This uncertainty leads us to the last issue that we want to briefly to mention in this work: whether the compact object in CG X-1 is more likely a BH or a NS. Based on X-ray luminosity alone, until a few years ago it would have been straightforward to classify it as a BH in the ultraluminous state. However, we now know from observational studies that NS ULXs do exist and can also reach apparent luminosities of ∼10⁴⁰ erg s⁻¹ (Bachetti et al. 2014; Fürst et al. 2016; Israel et al. 2017a, 2017b), although so far no NS ULX has been found with a period as short as a few hours. Theoretical models have also caught up with the observations to explain the existence and abundance of NS ULXs (Koliopanos et al. 2017; Pintore et al. 2017; Wiktorowicz et al. 2017, 2019; Walton et al. 2018; Mushtukov et al. 2019).

If CG X-1 is a wind-accreting rather than Roche-lobe overflow system, one argument in favor of the BH scenario is that a BH can intercept a larger fraction of stellar wind, and reach higher luminosities. For a typical WR wind with $-{\dot{M}}_{2}$ ∼ a few × 10⁻⁵ M_⊙ yr⁻¹ (e.g., Hillier 2003; Vink & de Koter 2005; Belczynski et al. 2010; Smith 2014), the fraction intercepted by the compact object must be at least ∼10% of that rate to explain the luminosity of CG X-1. The intercepted wind fraction can be estimated as ${\dot{M}}_{1}/-{\dot{M}}_{2}\approx {(1/4)({M}_{1}/{M}_{2})}^{2}\,{({R}_{2}/a)}^{2}$ (Frank et al. 2002). Only with a BH accretor (M₁ ∼ M₂) can this fraction approach ∼10%, perhaps with the aid of focused winds if the donor is close to filling its Roche lobe. A WR + NS system is more likely to have ${\dot{M}}_{1}/-{\dot{M}}_{2}\lesssim {10}^{-3}$, which is insufficient to explain CG X-1. In addition to the smaller accretion cross-section, there is another (possibly even more important) reason why a NS orbiting in the wind of a WR donor is not expected to become a ULX (Lipunov 1982). Most of the WR wind is expected to be deflected by the magnetosphere of the NS (i.e., the magnetospheric radius can be larger than the Bondi capture radius). Moreover, any phase of supercritical accretion will quickly spin-up the NS above the threshold for the propeller or ejector regimes (Lipunov 1982; Ergma & Yungelson 1998). The characteristic spin-down timescale is then longer than the WR lifetime, so that the NS would not have time to become a ULX again.

Another argument in favor of the BH scenario in short-period WR X-ray binaries (van den Heuvel et al. 2017) is that a NS would be unlikely to survive the common-envelope stage and would likely merge with the WR star. In particular, van den Heuvel et al. (2017) showed (their Figures 1, 2) that for a mass ratio M₂/M₁ ≳ 3.5, a 1.5 M_⊙ NS may survive common envelope only for a very limited range of initial binary periods and donor masses (M₂ ≈ 40–70 M_⊙); a 5 M_⊙ BH has the highest chance to survive and produce a short-period WR system such as Cyg X-3 and CG X-1, for initial donor masses ≳18 M_⊙. Meanwhile, more massive BHs (M₁ ≳ 15 M_⊙) may avoid common envelope altogether for any mass of the donor star, and go through a slower spiral-in phase during stable Roche-lobe overflow, with periods always ≳1 day (Pavlovskii et al. 2017; van den Heuvel et al. 2017; Bogomazov et al. 2018).