CG X-1: an eclipsing Wolf-Rayet ULX in the Circinus galaxy

TL;DR

This paper studies the eclipsing Wolf-Rayet ultraluminous X-ray source CG X-1 in the Circinus galaxy, revealing its orbital period, spectral properties, and suggesting it as a key object for understanding super-critical accretion and binary evolution.

Contribution

It provides the first detailed analysis of CG X-1's timing, spectral variability, and wind interactions, establishing it as a unique Wolf-Rayet ULX candidate.

Findings

Orbital period of approximately 7.2 hours with a measurable period derivative.

Spectral analysis indicating complex absorption by neutral, ionized, and Compton-thick material.

Luminosity varies between 4×10^{39} and 3×10^{40} erg/s, consistent with super-Eddington accretion.

Abstract

We investigated the time-variability and spectral properties of the eclipsing X-ray source Circinus Galaxy X-1 (GG X-1), using Chandra, XMM-Newton and ROSAT. We phase-connected the lightcurves observed over 20 years, and obtained a best-fitting period $P = (25,970.0 \pm 0.1)$ s $\approx$7.2 hr, and a period derivative $\dot{P}/P = ( 10.2\pm4.6) \times 10^{-7}$ yr$^{-1}$. The X-ray lightcurve shows asymmetric eclipses, with sharp ingresses and slow, irregular egresses. The eclipse profile and duration vary substantially from cycle to cycle. We show that the X-ray spectra are consistent with a power-law-like component, absorbed by neutral and ionized Compton-thin material, and by a Compton-thick, partial-covering medium, responsible for the irregular dips. The high X-ray/optical flux ratio rules out the possibility that CG X-1 is a foreground Cataclysmic Variable; in agreement with previous studies, we conclude that it is the first example of a compact ultraluminous X-ray source fed by a Wolf-Rayet star or stripped Helium star. Its unocculted luminosity varies between $\approx$4 $\times 10^{39}$ erg s$^{-1}$ and $\approx$3 $\times 10^{40}$ erg s$^{-1}$. Both the donor star and the super-Eddington compact object drive powerful outflows: we suggest that the occulting clouds are produced in the wind-wind collision region and in the bow shock in front of the compact object. Among the rare sample of Wolf-Rayet X-ray binaries, CG X-1 is an exceptional target for studies of super-critical accretion and close binary evolution; it is also a likely progenitor of gravitational wave events.

Figures

Tables

Table 1. Table 1 Log of the X-ray observations considered in this work

ObsID	Instrument	Start date	Exp. Time (ks)	Off-Axis Angle (′)	counts
rh702058a02 (0)	ROSAT/HRI	1997-03-03 21:50:52	26.38	0.30	75
rh702058a03 (0)	ROSAT/HRI	1997-08-17 10:44:29	45.89	0.30	169
355	ACIS-S	2000-01-16 10:18:17	1.32	1.56	43
365∗	ACIS-S	2000-03-14 06:01:26	4.97	1.09	1639
356∗	ACIS-S	2000-03-14 07:46:18	24.72	1.09	3579
374	HETG	2000-06-15 22:01:09	7.12	1.06	85
62877	HETG	2000-06-16 00:38:28	60.22	1.06	923
2454	ACIS-S	2001-05-02 16:02:48	4.40	0.73	301
0111240101∗ (1)	MOS+PN	2001-08-06 08:54:51	109.85	0.27	15791
4770	HETG	2004-06-02 12:40:42	55.03	1.32	591
4771	HETG	2004-11-28 18:26:32	58.97	1.11	603
9140∗	ACIS-S	2008-10-26 10:24:46	48.76	4.28	1176
10226	HETG	2008-12-08 17:57:06	19.67	1.35	201
10223	HETG	2008-12-15 15:46:15	102.93	1.34	642
10832	HETG	2008-12-19 18:15:08	20.61	1.34	187
10833	HETG	2008-12-22 07:29:35	28.36	1.34	211
10224	HETG	2008-12-23 11:25:12	77.10	1.33	483
10844	HETG	2008-12-24 23:17:37	27.17	1.33	290
10225	HETG	2008-12-26 04:02:06	67.89	1.33	651
10842	HETG	2008-12-27 12:03:26	36.74	1.33	332
10843	HETG	2008-12-29 10:10:49	57.01	1.32	448
10873	HETG	2009-03-01 23:28:35	18.10	1.12	151
10850	HETG	2009-03-03 00:43:18	13.85	1.12	77
10872	HETG	2009-03-04 15:29:52	16.53	1.12	83
10937∗	ACIS-S	2009-12-28 21:10:27	18.31	2.96	464
12823∗	ACIS-S	2010-12-17 18:10:27	152.36	1.52	10464
12824∗	ACIS-S	2010-12-24 03:38:54	38.89	1.52	2368
0701981001∗ (2)	MOS+PN	2013-02-03 07:24:11	58.91	4.80	3177
0656580601∗ (3)	MOS+PN	2014-03-01 09:55:41	45.90	0.27	1833
0792382701∗ (4)	MOS+PN	2016-08-23 16:53:33	37.00	4.70	5939
0780950201∗ (5)	MOS+PN	2018-02-07 12:11:49	44.36	4.70	4302

Table 2. Table 2 Best-fitting parameters of phase-resolved spectra from combined Chandra  ObsIDs 9140, 10937, 12823, and 12824

Phase	$N_{H,\mathrm{intr}}$	$\Gamma$/ kT (keV)	p	$N_{H,\mathrm{absori}}$	$F_X$	$L_X$	${\chi }_{\nu }^2$ (dof)
	$({10}^{22}$ cm${}_{}{}^{-2})$			$({10}^{22}$ cm${}_{}{}^{-2})$	$({10}^{-13}$ erg cm-2 s-1)	$({10}^{38}$ erg s-1)
tbabs${}_1{}^{}\times$tbabs${}_2{}^{}\times$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		2.42${}_{-0.29}{}^{}{}_{}{}^{+0.38}$	-	-	0.57${}_{-0.07}{}^{}{}_{}{}^{+0.06}$	3.3${}_{-0.6}{}^{}{}_{}{}^{+1.7}$	0.76 (30)
tbabs${}_1{}^{}\times$tbabs${}_2{}^{}\times$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}$	$>24.3$	-	-	5.8${}_{-0.3}{}^{}{}_{}{}^{+0.3}$	17.4${}_{-0.9}{}^{}{}_{}{}^{+1.0}$	0.84 (100)
Faint		2.34${}_{-0.53}{}^{}{}_{}{}^{+0.45}$	-	-	0.52${}_{-0.06}{}^{}{}_{}{}^{+0.04}$	2.27${}_{-0.16}{}^{}{}_{}{}^{+0.17}$	0.91 (30)
tbabs${}_1{}^{}\times$tbabs${}_2{}^{}\times$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		2.50${}_{-0.03}{}^{}{}_{}{}^{+0.04}$	-	-	6.0${}_{-0.4}{}^{}{}_{}{}^{+0.4}$	15.6${}_{-0.8}{}^{}{}_{}{}^{+0.8}$	0.82 (100)
Faint		0.80${}_{-0.09}{}^{}{}_{}{}^{+0.10}$	-	-	0.46${}_{-0.05}{}^{}{}_{}{}^{+0.05}$	1.76${}_{-0.13}{}^{}{}_{}{}^{+0.13}$	1.36 (30)
tbabs${}_1{}^{}\times$tbabs${}_2{}^{}\times$diskpbb
Bright	0.49${}_{-0.05}{}^{}{}_{}{}^{+0.06}$	7.1${}_{-1.2}{}^{}{}_{}{}^{+\ast }$	0.54${}_{-0.01}{}^{}{}_{}{}^{+0.01}$	-	17.0${}_{-0.4}{}^{}{}_{}{}^{+0.4}$	60.1${}_{-3.0}{}^{}{}_{}{}^{+9.2}$	0.86 (266)
Intermediate		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		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)
tbabs${}_1{}^{}\times$absori$\times$tbabs${}_2{}^{}\times$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)
tbabs${}_1{}^{}\times$tbabs${}_2{}^{}\times$bbodyrad
Bright		0.87${}_{-0.01}{}^{}{}_{}{}^{+0.01}$	-	-	14.0${}_{-0.3}{}^{}{}_{}{}^{+0.3}$	35.7${}_{-0.7}{}^{}{}_{}{}^{+0.7}$	1.76 (267)
Intermediate		0.99${}_{-0.04}{}^{}{}_{}{}^{+0.04}$	-	-	5.0${}_{-0.3}{}^{}{}_{}{}^{+0.3}$	12.4${}_{-0.6}{}^{}{}_{}{}^{+0.7}$	1.46 (100)
Faint		0.49${}_{-0.05}{}^{}{}_{}{}^{+0.05}$	-	-	0.40${}_{-0.04}{}^{}{}_{}{}^{+0.05}$	1.4${}_{-0.1}{}^{}{}_{}{}^{+0.1}$	2.25 (30)

Table 3. Table 3 Best-fitting parameters for Chandra  and XMM-Newton  phase-resolved spectra, for our fiducial model.

Phase	$N_H$	$\xi$	f	$N_H$	$\Gamma$	Norm	$F_X$	$L_X$	${\chi }_{\nu }^2$ (dof)
	(${10}^{22}$ cm-2)	(${10}^2$)		(${10}^{22}$ 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}$	$\left[0\right]$	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}$	$\left[0\right]$	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.	$\left[0\right]$	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}$	$\left[0\right]$		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}$				9.1 ${}_{-1.5}{}^{}{}_{}{}^{+6.2}$
XMM-Newton obs. 0792382701 (4)
Brightb	3.21${}_{-0.88}{}^{}{}_{}{}^{+2.49}$	unconstr.	$\left[0\right]$	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}$				14.8${}_{-1.1}{}^{}{}_{}{}^{+1.6}$
XMM-Newton obs. 0780950201 (5)
Brightb	0.41${}_{-0.41}{}^{}{}_{}{}^{+2.08}$		$\left[0\right]$	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}$

Table 4. Table 4 Comparison of power-law and broken power-law models in the 2–6 keV band ( XMM-Newton  observation 1)

$\Gamma$	${\chi }^2$/dof	${\Gamma }_1$	$E_{\mathrm{break}}$	${\Gamma }_2$	${\chi }^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

Table 5. Table 5 Best-fitting parameters of the phase-averaged spectra from a sample of individual observations.

Obs.	$N_H$	$N_H$	$\Gamma$	Norm	$⟨F_X⟩$	$⟨L_X⟩$	$L_{X,\mathrm{peak}}$	${\chi }_{\nu }^2$ (dof)
	$({10}^{22}$ cm${}_{}{}^{-2})$	$({10}^{22}$ cm${}_{}{}^{-2})$		$\left({10}^{-4}\right)$	(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)

Equations

$$t_{\mathrm{c}}=\frac{{\sum }_{i=1}^N{\gamma }_i\left(I_0-I_i\right)t_i}{{\sum }_{i=1}^N{\gamma }_i\left(I_0-I_i\right)},$$

$$W={\left|\frac{{\sum }_{i=1}^N{\gamma }_i\left(I_0-I_i\right){\left(t_i-t_{\mathrm{c}}\right)}^2}{{\sum }_{i=1}^N{\gamma }_i\left(I_0-I_i\right)}\right|}^{1/2}.$$

$$O-C=(P_0-P_{\mathrm{best}})\times N,$$

$$O-C=c_1+c_{2}N+(1/2)P_{\mathrm{best}}\dot{P}\times N^2,$$

$$a$$

$$R_{\mathrm{L},2}=a\frac{0.49q^{2/3}}{0.6q^{2/3}+\ln (1+q^{1/3})}$$

$$\rho =0.2101\frac{q}{1+q}{\left(\frac{R_{\mathrm{L},2}}{a}\right)}^{-3}{\left(\frac{P}{7.2\mathrm{hr}}\right)}^{-2}\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}$$

$${\rho }_{w2}$$

$$t_{\mathrm{cool},2}=\frac{3(n_e+n{)}_{s2}k{T}_{s2}}{2n_{s2}^2\Lambda (T_{s2})}\approx \frac{3k{T}_{s2}}{n_{s2}\Lambda (T_{s2})}$$

$$v_{K1}={\left[\frac{G{M}_2^2}{a\left(M_1+M_2\right)}\right]}^{1/2}.$$

$$R_0={\left(\frac{{\dot{M}}_1{v}_{w1}}{4\pi {\rho }_{\mathrm{a}}{v}_{K1}^2}\right)}^{1/2}$$

$$R_0={\left(\frac{{\dot{M}}_1}{{\dot{M}}_2}\right)}^{1/2}{\left(\frac{v_{w1}{v}_{w2}}{v_{K1}^2}\right)}^{1/2}a.$$

Full text

1. Introduction

Ultraluminous X-ray sources (ULXs) are non-nuclear point-like sources with X-ray luminosity $\gtrsim 3\times$1039 erg s*-1* (Kaaret et al., 2017; Feng & Soria, 2011). The luminosity distribution of this population is now fairly well determined (Wang et al., 2016; Mineo et al., 2012; Walton et al., 2011; Swartz et al., 2011; Liu & Bregman, 2005), and is consistent with the high-luminosity end of the high-mass X-ray binary (XRB) distribution; we also know that the number of ULXs in a star-forming galaxy is roughly proportional to its star formation rate, and that the distribution may have a cut-off or downturn at an X-ray luminosity of $\approx$2 $\times$ 1040 erg s*-1* (Mineo et al., 2012; Swartz et al., 2011). However, more specific properties of these sources are still poorly constrained. It is not known what fraction of them are powered by a BH and what fraction by a neutron star (NS) (see e.g., Wiktorowicz2019, and references therein); so far, an identification of the compact object has been possible only for a handful of ULXs with X-ray pulsations, signature of a magnetized NS (Bachetti et al., 2014; Israel et al., 2017a, b; Fürst et al., 2016; Tsygankov et al., 2017; Carpano et al., 2018). The relative distribution of stellar types and ages for the donor stars is also poorly known. In a few cases, the donor is identified as a blue supergiant (Motch et al., 2014), or a red supergiant (Heida et al., 2016), or a low-mass star (Soria et al., 2012), but in many other cases, it is hard to tell whether the observed optical counterpart corresponds to the donor star or the irradiated accretion disk (Tao et al., 2011). Likewise, the mechanism of mass transfer (Roche lobe overflow (RLOF) or wind accretion), the duty cycle, the duration of the super-Eddington phases and the total amount of mass that may be accreted by the compact object during its lifetime remain generally unknown.

Population synthesis models predict super-Eddington mass transfer phases from several different types of donor stars at different ages, but our lack of empirical information on the fundamental system parameters for most ULXs makes it difficult to test such models. It also makes it hard to predict the future evolution of such binary systems, e.g., what fraction of ULXs will evolve into double compact binaries (BH-BH, BH-NS or NS-NS), potential progenitors of gravitational wave mergers (Marchant et al., 2017). Perhaps the most important piece of information that is usually missing is the binary period; without a period, also the mass ratio and the binary separation remain unconstrained. Only in a few cases have periodic variations (in either the X-ray or the optical lightcurve) been detected and interpreted as a binary period; those periods vary between a few days to a few months (Liu et al., 2013; Bachetti et al., 2014; Motch et al., 2014; Fürst et al., 2018; Urquhart & Soria, 2016; Wang2018; Vinokurov et al., 2018).

One rare ULX candidate with a well-determined period of 7.2 hr (as well as other intriguing X-ray properties) is the very bright X-ray source CG X-1 (Bauer et al., 2001; Weisskopf et al., 2004; Esposito et al., 2015). It is seen projected in the sky inside the inner region of the Circinus galaxy, $\approx$15*′′* to the east of its nucleus ($\approx$300 pc at the distance of Circinus), with the coordinates of $\alpha=14^{\rm h}13^{\rm m}12^{\rm s}.21$, $\delta=-65^{\circ}20^{\prime}13^{\prime\prime}.7$ (J2000). Circinus is a large Sb galaxy at a distance of 4.2 Mpc (Tully et al., 2009) with a Seyfert nucleus and intense star formation, at a rate of $\approx$3–8 $M_{\odot}$ yr*-1* (Freeman et al., 1977; For et al., 2012). If CG X-1 belongs to the Circinus galaxy, its X-ray luminosity is $\approx$1040 erg s*-1* (reaching $\approx$3 $\times 10^{40}$ erg s*-1* at some epochs), making it one of the most luminous ULXs in the local universe, near the potential break in the ULX luminosity function. By comparison, from the average X-ray luminosity function of Mineo et al. (2012), we expect $\approx$0.2–0.6 sources at or above an X-ray luminosity of 1040 erg s*-1* in a galaxy with the star formation rate of Circinus. Thus, the presence of 1–2 luminous ULX in Circinus is not unexpected. Besides CG X-1, there is another ULX (known as ULX5) in the outskirts of this galaxy, with $L_{\rm X}\approx 2\times 10^{40}$ erg s*-1* (Walton et al., 2013).

There has been some debate in the literature about whether CG X-1 really belongs to Circinus or is instead a foreground magnetic cataclysmic variable (mCV) in the Milky Way (as suggested by Weisskopf et al. 2004), projected by chance in front of Circinus. In the latter case, the 7.2-hr period would correspond to the spin period of the white dwarf rather than the binary period. There were several reasons behind the mCV suggestion. Firstly, Circinus is located almost behind the Galactic plane, in a field crowded with foreground stars (Figure 1); in fact, two other X-ray sources projected in the field of Circinus (but not as close to its starburst region as CG X-1) turned out to be foreground mCVs (Esposito et al., 2015). Secondly, when first discovered two decades ago, empirical knowledge of ULXs was still scant, and an intrinsic luminosity in excess of $10^{40}$ erg s*-1* for a non-nuclear source was still regarded as suspiciously unlikely. Thirdly, the X-ray lightcurve of CG X-1, with its asymmetric eclipses, is very unusual for a ULX or more generally for a luminous X-ray binary and is similar instead to the periodic eclipses of the accreting poles in an mCV. However, the mCV interpretation was substantially refuted by Esposito et al. (2015), based on probability arguments. They also showed that the short orbital period and eclipsing lightcurve are typical of systems with a Wolf-Rayet (WR) donor stars (as already proposed by Bauer et al. 2001); a few such systems have been discovered in recent years (Carpano et al., 2007; Esposito et al., 2013, 2015; Maccarone et al., 2014), although none of them has reached the super-Eddington luminosity of CG X-1. If CG X-1 survives the collapse of the WR donor and evolves into a double BH binary, the coalescence timescale via gravitational wave emission is only $\approx$50 Myr (Esposito et al., 2015).

Here, we present a detailed, multi-epoch analysis of CG X-1 to further probe the nature of this remarkable system. The main objectives of our paper are the following. In Section 2, we will describe the X-ray data we used. In Section 3, we will illustrate the ever-changing lightcurve profiles observed from individual orbital cycles, using data from two ROSAT, 24 Chandra and five XMM-Newton observations of Circinus, obtained over a period of 21 years; we will attempt to distinguish between periodic features and stochastic dipping, and we will obtain a more precise measurement of the period and of the period derivative. In Section 4, we will model the X-ray spectral properties of the systems during peak-luminosity phases, during dipping phases, and during occultation phases (which still show faint, residual emission). We corroborate the ULX interpretation of Esposito et al. (2015), and in this paper we will only add a short discussion about the X-ray to optical flux ratio of CG X-1, which again disfavours the mCV scenario (Section 5). In Section 6, we will propose a physical scenario that may explain the puzzling timing and spectral properties, in terms of partial covering by optically thick clumps located mostly in front of the compact object, as it moves through the dense wind of the donor star (bow shock scenario). Then, we will briefly discuss the possible origin and evolutionary state of this system in the context of population synthesis models. Finally (Section 7), we present a summary of our results.

2. Observations and data analysis

2.1. Chandra

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 centred 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 provided by the Chandra Interactive Analysis of Observations (ciao) software package, Version 4.9 (Fruscione et al., 2006). We corrected 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, to determine the coordinates of CG X-1 and to identify an ellipse around CG X-1 that contains $\approx$99% of the source counts (wavdetect parameter ellsigma $=3$). The typical semi-major axis of the elliptical source region is \approx 1\hbox{{}^{\prime\prime}}.3. 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 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 lightcurves 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 $\chi^{2}$ fitting.

For the HETG observations, we only extracted lightcurves (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 1*′′.3, 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).

2.2. XMM-Newton

XMM-Newton has observed CG X-1 on six occasions before October 2018, 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 lightcurves 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 count rate less than some certain criteria were remained to create the GTI intervals. The count rate criteria varies according to the observations and cameras, typically higher than 0.8 ct s*-1* for pn and 0.3 ct s*-1* for MOS, respectively. For the first four observations, we used the resulting GTI intervals for both lightcurve 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, 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 ($\approx$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 centred 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 lightcurve and the spectral extraction.

We used “#xmmea_em && (pattern $<=$ 12)” to filter MOS lightcurves and spectra; used “#xmmea_ep && (pattern $<=$ 4 && flag==0)” to filter the PN data. We extracted lightcurves from MOS and pn with the sas tasks evselect, followed by epiclccorr to correct for vignetting and bad pixels111http://www.cosmos.esa.int/web/xmm-newton/sas-threads. 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 $\chi^{2}$ fitting. Finally, as for the Chandra data, we used xspec for spectral modelling.

2.3. ROSAT

CG X-1 has been observed by the Roengten Satellite High Resolution Imager (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–11, and between 1997 August 17–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 lightcurve 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 lightcurves using the best-fitting period (determined from Chandra and XMM-Newton observations), for a better statistics. The epoch of this folded lightcurve was chosen at the orbital cycle, where the two observations accumulated half of the total observed counts, i.e., MJD 50681.437036 for phase $\phi$ = 1 (the phase will be defined in Section 3.1). The main purpose for using the folded ROSAT lightcurve is to study the variation of the periodicity of CG X-1 over a long-term duration of more than 20 years.

3. X-ray Timing Results

3.1. Binary period

A $(27\pm 0.7)$ ks X-ray periodicity was first discovered by Bauer et al. (2001) from Chandra/HETG zeroth-order data, and 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\pm 0.2)$ ks. To reduce the error even further, for this work we re-analyzed the whole archival set of twenty Chandra observations and five XMM-Newton observations longer than about half of the period ($\approx$13 ks): they cover a time span of 19 years (Table 1). (We did not use an additional four archival Chandra observations because their duration was only $\lesssim$0.25 times the period).

For each Chandra or XMM-Newton observation, we generated a 0.3–8 keV background-subtracted lightcurve, and rebinned each lightcurve to 500-s bins (Figure 2). We then applied the phase dispersion minimization (PDM) method to compute the best-fitting period (Stellingwerf, 1978; Schwarzenberg-Czerny, 1997). 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 archival observations (Figure 2), with a best-fitting period $P_{\rm best}=(25970.0\pm 0.1)$ s: an improvement in precision by a factor of 2000.

We computed six folded lightcurves for different sub-groups 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 lightcurves 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 lightcurves simultaneously. We re-defined phase $\phi\equiv 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 lightcurves: we estimate that the deepest points of the six folded eclipses have only a small scatter of $\delta\phi\approx 0.034$ around the global best-fitting phase $\phi=1$. (In Section 3.2, we will investigate this scatter further, and look for possible small changes in the period over the last 20 years.) This confirms that the X-ray eclipse is related to stable properties of the binary system: the most plausible explanation is that $\phi=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.

In summary, assuming a constant period $P_{\rm best}$, our best-fitting ephemeris is $\phi=1$ at MJD $50681.437036+0.300579\times N$ (d). We have chosen to set the reference time ($N=0$) at an epoch covered by the ROSAT observations222The ROSAT observations span several time intervals between 1997 March and September. We arbitrarily defined the cycle $N=0$ as the one in which we reached 50% of the counts in the stacked dataset.. 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.

3.2. Period derivative

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 $\phi=1$) is relatively easy to determine in all the lightcurves averaged over many cycles or several years (Figure 3), it is not obvious how to identify this point in any of the individual lightcurves (Figure 2), with eclipse durations and egress behaviours that differ markedly from cycle to cycle. Instead, we used the eclipse ingress as a phase marker for the individual lightcurves; the sharp flux drop leading to the eclipse is a feature that can be unequivocally identified in most of the individual lightcurves. 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, possible evidence of the period derivative $\dot{P}\neq 0$.

First, we determined the average mid-time of the ingress from the stacked Chandra and XMM-Newton lightcurves. We used the method developed by (Hu et al., 2008), which we adapted for our purposes. We defined the count rate $I_{0}$ of the un-eclipsed high state as the average count rate from $\phi=0.5$ to $\phi=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 lightcurve between two successive bins at times $t_{i}$ and $t_{i+1}=t_{i}+500$s. Finally, we defined a mid-time $t_{\rm c}$ of the ingress section of a lightcurve as

[TABLE]

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<\phi<1$), and $t_{i}$ and $I_{i}$ are the time and count rate of the $i$th bin, respectively. Eq.(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 weighted also by $\left(I_{0}-I_{i}\right)$ is because we want to reduce the uncertainty in the estimate of the time at which the lightcurve rolls over and the ingress starts; near the beginning of the ingress, $\left(I_{0}-I_{i}\right)\approx 0$. In any case, Eq.(1) provides a practical definition of a lightcurve feature that can be used to search for first-order period changes. Similarly, the physical width of the ingress is calculated as:

[TABLE]

The mid-ingress phase for the stacked Chandra and XMM-Newton observations occurs at phase $\phi=0.85\pm 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\times N$ (d) $\equiv T_{1}+P_{\rm best}\times N$. The empirically determined location of $t_{\rm c}$ in each of the individual lightcurves plotted in Figure 2 is marked with a magenta line.

In order to investigate whether the period is changing with time over a 20-year 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_{1}+P_{\rm best}\times N$), and are assumed to have no error. We calculated the O values from the data, using Eq.(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 only by two weeks), to increase the signal to noise ratio; in that case, we measured the time difference between the average ingress phase in the stacked lightcurve (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 datapoint was used for each observation, defined as the average of the individual time differences for all the ingresses covered during that observation. We estimate the 1-$\sigma$ uncertainty (standard deviation) for each individual ingress measurement as $\approx$700s; smaller errors are associated to datapoints that are the average of multiple ingress times.

In an $O-C$ diagram, if the period remains constant with time, datapoints are scattered along a line:

[TABLE]

where $P_{0}$ is the true period and $P_{\rm{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 follows a horizontal line around zero. Instead, if the orbital period has a small linear change, the datapoints follow a quadratic function:

[TABLE]

where $c_{1}$ and $c_{2}$ are normalization coefficients, and $\dot{P}$ is the period derivative.

The resulting $O-C$ diagram is shown in Figure 4; it visually suggests an upwards curvature. We fitted the $O-C$ datapoints with a quadratic function. The best fitting curve (dashed-dotted line in Figure 4) is $O-C=1320.6-0.29\times N+1.09\times 10^{-5}\times N^{2}$. The orbital period derivative is $\dot{P}=(8.4\pm 3.8)\times 10^{-10}$ s s*-1*, where the error range is the 90% confidence level (2.7$\sigma$); this corresponds to $\dot{P}/P=(10.2\pm 4.6)\times 10^{-7}$ yr*-1* at the 90% confidence level. The value of $\dot{P}/P$ is significantly $>$ 0 to approximately 10$\sigma$ (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\approx 1.1\times 10^{-3}$) confirms that the O-C datapoints 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 years.

3.3. Duration of eclipse and egress phases

Lightcurves folded over several cycles (Figure 3) show an eclipse lasting from $\phi\approx 0.9$ to $\phi\approx 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 $\delta\phi\approx$ 0.15. At other epochs, such as Chandra ObsID 9140 and ObsID 12823, the observed flux remains close to zero for a longer time, $\delta\phi\approx 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 behaviour 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 we shall discuss 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 indeed it 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 $\phi=0.8$–$1.8$. Then, we defined a duration of the occultations as the number of 500-s phase bins ($N_{\rm{occ}}$) in which the count rate is lower than 40% of the maximum count rate defined above. These bins typically include the candidate stellar eclipse around phase 1, as well as deep dips during the egress, 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 $\Delta\phi_{\rm{occ}}=N_{\rm{occ}}\times 500/P_{\rm{best}}\approx 0.0193N_{\rm{occ}}$. The uncertainty on the half-width was calculated as a Poisson error on the number of phase bins (i.e., $0.0193\sqrt{N_{\rm occ}}$).

Finally, in order to test whether the duration of the occultations is correlated with the peak luminosity, we used the online Portable, Interactive Multi-Mission Simulator333http://cxc.harvard.edu/toolkit/pimms.jsp. (pimms) version 4.9 to convert all 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 used for the conversion is a power-law with photon index $\Gamma=1.7$ and absorbing column density $N_{\rm H}=1\times 10^{22}\ cm^{-2}$.

We found a correlation (Figure 5) between the duration of the occultations and the peak luminosity, represented by the normalized maximum count rate $CR_{\rm max}$. At low luminosities ($L_{\rm X}\lesssim 10^{40}$ erg s*-1*), the duration of the occultations is linearly anticorrelated with the count rate, as $\Delta\phi_{\rm{occ}}=(-1.96\pm 0.97)\times CR_{\rm max}+(0.70\pm 0.13)$. To quantify the statistical significance of the correlation, we calculated the Spearman correlation coefficient $\rho$: we obtained $\rho\approx-0.73$, with a p value of 5$\times 10^{-6}$. At higher luminosities ($L_{\rm X}>10^{40}$ erg s*-1*), the duration of the occultations saturates around $\Delta\phi_{\rm{occ}}\approx 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.

4. X-ray Spectral Results

4.1. Outline of our spectral analysis

The next question we addressed is whether/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 three steps, as explained below.

First (Section 4.2), we split a selected sample of high signal-to-noise lightcurves (from Chandra/ACIS observations) into phenomenological sub-structures (ingress, eclipse, dips, egress, bright phase) and compared the cumulative energy distribution of the observed photons in the various sub-structures. 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 modelling 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 lightcurve 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; 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; 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 sub-intervals, 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 modelled spectra integrated over entire orbital cycles, from six Chandra observations and five XMM-Newton ones, between 2000 and 2018; the datasets used for this modelling are marked by asterisks in Table 1). We determined average and peak fluxes and luminosities during those observations.

4.2. Model-independent photon energy distribution

The four Chandra observations chosen for our model-independent study of the photon energy distribution are ObsIDs 9140, 10937, 12823 and 12824, covering a total of about 11 orbital cycles444We did not use Chandra ObsID 365, because its exposure time is too short to enable a meaningful definitions of phase substructures. Also, we did not use Chandra ObsID 356 because of its high pile-up fraction, $\approx$40%. However, we did model the individual spectra of those two observations among the others in Section 4.5.. The advantage of those particular observations is that it is relatively straightforward to identify five main sub-structures in the background-subtracted lightcurves: faint phase, bright phase, ingress, egress and dips; the five sub-structures are colour-coded in Figure 6. More specifically, for this part of the analysis we defined a maximum count rate as the average value of the ten 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; the bright phase as the time intervals when the count rates are higher than 70% of the maximum count rate; the ingress and egress phases 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 $\lesssim 2000$ s; they are colour-coded as red datapoints 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 took the photons falling within $\pm$200 s of that mid time, in the unbinned event file.

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\times 10^{-8}$. This suggests the presence of at least two emission components: a bright, harder one, and a faint, softer one that is seen as a residual component when the other component is occulted.

4.3. Chandra spectra

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 shall refer to them 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_{\rm H}$ fixed at $5.6\times 10^{21}$ cm*-2*, 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 one, especially for the faint phase, which has an unacceptable $\chi^{2}_{\rm\nu}=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 $\approx$1–5 keV range, but a lower intrinsic $N_{\rm H}$. Conversely, the faint spectrum (dominated by residual emission in the eclipse) is significantly softer (steeper).

Instead of a simple power law, we also tried Comptonization models (in particular, comptt, simpl $\times$ diskbb, and diskir). However, such 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 $\sim$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, from ionized gas. We do that 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_{\rm H}\sim$ a few $10^{22}$ cm*-2* and an ionization parameter $\sim$ 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^{2}_{\nu}=232.7/267$ to $\chi^{2}_{\nu}=218.7/265$: this is significant to $>$99.9% probability555 Model tbabs$\times$ absori $\times$ tbabs $\times$ po is equivalent to model tbabs$\times$ tbabs $\times$ po, when $N_{\rm H}$ from the absori component is equal to zero. The significance level is calculated using a $\chi^{2}$ distribution with two degrees of freedom, i.e. $\xi$ and $N_{\rm H}$ from the absori component. . For the intermediate spectrum, the improvement is from $\chi^{2}_{\nu}=85.0/100$ to $\chi^{2}_{\nu}=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$\sigma$ level ($N_{\rm H}=2.3^{+4.1}_{-2.3}\times 10^{22}$ cm*-2*).

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 ($\Gamma\approx 2.0$); the reason 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 ($\Gamma\approx 2.6$) than the other two (Table 2).

We can now start to identify some 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. There are changes in the fitted column densities of ionized and neutral components, but they only affect the shape of the spectrum below 2 keV. Instead, to 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 such evolution and 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 $\times$ tbabs $\times$ absori $\times$ pcfabs $\times$ 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 modelled with pcfabs has to be Compton thick, with $N_{\rm H}\gtrsim 2\times 10^{24}$ cm*-1*, 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^{2}_{\nu}=297.7/363=0.82)$. The photon index of the intrinsic emission component is $\Gamma\approx 2.1$ and the partial covering fraction of the intermediate spectrum is $f\approx 0.60$. The average intrinsic (de-absorbed) luminosity during the four Chandra observations used in our spectral modelling is $L_{\rm X}\approx 9\times 10^{39}$ erg s*-1*, in the 0.3–8 keV band. (Higher luminosities in excess of 1040 erg s*-1* 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. It cannot be explained with the same $\Gamma\approx 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\equiv 1$ in pcfabs) and there is instead a different emission component significantly detectable only in eclipse. We model this component also with a simple power-law, for want of higher signal-to-noise spectra; we find a photon index $\Gamma\approx 2.6$ and a luminosity L_{\rm X}$$\approx 5\times 10^{38} erg s*-1* in the 0.3–8.0 keV band (Table 3).

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 that, 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.

4.4. XMM-Newton spectra

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$\times$absori$\times$pcfabs$\times$tbabs$\times$po. The first tbabs component represents the line-of-sight absorption and is fixed at $N_{\rm H}=5.6\times 10^{21}$ cm*-2*; the second tbabs and the absori components represent Compton-thin neutral and ionized absorption around CG X-1, respectively; 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 $\sim$2–5 $\times 10^{21}$ cm*-2* and ionized-absorber column densities $\sim$2–4 $\times 10^{22}$ cm*-2* (Table 3), consistent with the properties of the Chandra spectra. The covering fraction for the intermediate spectra varies between $\sim$0.3–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 $\Gamma\approx 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 $\approx$5.5 $\times 10^{39}$ erg s*-1* to $\approx$2.6 $\times 10^{40}$ erg s*-1*, in the 0.3–8 keV band.

Finally, we used the bright-phase XMM-Newton spectra to test for the presence of a high-energy downturn, 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 so, 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 there. 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\pm 0.3)$ keV, and the continuum slope steepens from $\Gamma\approx 2.0$ to $\Gamma\approx 2.8$ (Figure 8) .

4.5. Long-term spectral and luminosity variations

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 so, we used six Chandra and five XMM-Newton observations from 2000 to 2018, as already mentioned in Section 4.1. This time, we extracted and fitted the average spectrum of each observation (i.e., no longer split into faint, intermediate and bright intervals); 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 $\times$ absori $\times$ tbabs $\times$ powerlaw. Here, we did not include a Compton-thick pcfabs term because we assume that it gives a grey 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). Then, we 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.

The highest peak luminosity was seen in the first XMM-Newton observation (2001 August 6), at $L_{\rm X,peak}\approx 2.9\times 10^{40}$ erg s*-1*; the lowest luminosity was recorded in Chandra ObsID 9140 (2008 October 26), with $L_{\rm X,peak}\approx 3.5\times 10^{39}$ erg s*-1* (Table 5). Thus, $L_{\rm X,peak}$ varies by a factor of $\approx$8 across our sub-sample of eleven observations (Figure 9), whereas the average luminosity of each observation varies by a factor of $\approx$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 $\Gamma$ is approximately the same for all observations. For 10 out of the 11 observations in our spectral analysis sub-sample, i.e., all except Chandra ObsID 365, the mean photon index is $\left\langle\Gamma\right\rangle\approx 2.11$ with a scatter $\sigma\approx 0.15$. Only the spectrum of Chandra ObsID 365 is significantly different from all the others, with $\Gamma=1.5\pm 0.2$.

5. The Nature of the donor star

5.1. A possible optical counterpart

A faint, point-like optical source was found at the X-ray position (Bauer et al., 2001; Weisskopf et al., 2004), in a combined 600s exposure from HST WFPC2, in the F606W filter (Figure 10). We re-processed the data to check or improve the astrometric alignment between X-rays and optical bands, and 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 $\approx$0.*′′*2 (white circle in Figure 10), improving on the previous results. The error circle does indeed include the optical source previously suggested as 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 zeropoints666http://www.stsci.edu/hst/acs/analysis/zeropoints, we confirm an apparent brightness m${}_{\rm F606W}\approx V\approx(24.3\pm 0.1)$ mag in the Vega system, in agreement with Weisskopf et al. (2004).

Taking into account the line-of-sight extinction $A_{F606W}\approx 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}\lesssim-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: Crowther 2007; Massey 2003). Given that CG X-1 is located in a star-forming region of Circinus, near dust lanes and filaments, we expect also 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$\alpha$ filter F656N and O III filter F502N, but no excess line emission is detected in those bands. In conclusion, with only one image in one optical filter, we do not have enough information to constrain the colour, spectral type and luminosity of the optical counterpart.

5.2. Ruling out foreground Galactic sources

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_{\rm X}/F_{\rm 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_{\rm X}$(0.1–4.0 keV)$/F_{\rm 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 relation777This relation comes from the standard definition of the apparent $V$ magnitude as $V=-2.5\log F_{\lambda}-21.100$ where $F_{\lambda}$ is the flux density at the top of the Earth’s atmosphere, in units of erg cm*-2* s*-1* Å*-1* (Bessell et al., 1998). Then, $F_{\rm opt}\equiv F_{\lambda}\times 1000$ Å. The same relation $m_{\lambda}=-2.5\log F_{\lambda}-21.100$ defines the “STmag” photometric system. $\log F_{\rm opt}=-0.4V-5.44$, we infer $F_{\rm X}$(2–10 keV)$/F_{\rm opt}\lesssim 7$ (Figure 3 in Mukai 2017), for both magnetic and non-magnetic CVs888The only exception to this rule is a small number of ultracompact white dwarf - white dwarf systems (Solheim, 2010), such as RX J1914$+$24, with orbital periods $\lesssim$10 min, which can reach $F_{\rm X}/F_{\rm opt}\sim 1000$: (Ramsay, 2008; Ramsay et al., 2005, 2002). However, such (very rare) sources have a super-soft, thermal X-ray spectrum with almost no emission above 1 keV, inconsistent with the spectrum of CG X-1. Their orbital period is also much shorter than the one measured in CG X-1, so there is no possibility of mis-identification.. A similar analysis by Revnivtsev et al. (2014) shows that $F_{\rm X}$(0.5–10 keV)$/F_{\rm opt}\lesssim 5$ for non-magnetic CV (for this relation, we have used the conversion $\nu F_{\nu}\approx 5F_{\lambda}\times 1000$ Å$\equiv F_{\rm opt}$ at $\lambda\approx 5000$ Å). All those 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_{\rm X}/F_{\rm opt}$ ratios, and compared them with the upper limits found in the CV surveys cited earlier. We have not removed the additional intrinsic absorption component derived from X-ray fitting, because that component is not removed in those CV surveys, either.

First, we compare the observed fluxes. From the net count rate in the WFPC2 F606W band, and the tabulated value of photflam999http://www.stsci.edu/hst/wfpc2/analysis/wfpc2_photflam.html, we obtain $F_{\lambda}\approx 5.2\times 10^{-19}$ erg cm*-2* s*-1* Å*-1*, and (by our definition, assuming a flat spectrum) $F_{\rm opt}\approx 5.2\times 10^{-16}$ erg cm*-2* s*-1*. 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_{\rm X}\approx 1.9\times 10^{-12}$ erg cm*-2* s*-1*, giving an implausibly high X-ray over optical flux ratio of $\approx 3650$. For a more physical result, we remove the line-of-sight extinction ($A_{V}\approx 4$ mag) from the optical flux, and the line-of-sight absorption (corresponding to $N_{\rm H}\approx 6\times 10^{21}$ cm*-2*) from the X-ray flux. We obtain $F_{\rm opt}\approx 2.06\times 10^{-14}$ erg cm*-2* s*-1* and $F_{\rm X}\approx 2.03\times 10^{-12}$ erg cm*-2* s*-1*, corresponding to $F_{\rm X}$(2–10 keV)$/F_{\rm opt}\approx 100$. These values are well outside the range observed in CVs, which are then definitively ruled out. Instead, $F_{\rm X}/F_{\rm opt}\sim 100$–1000 is what we expect and observe in typical ULXs (Ambrosi & Zampieri, 2018; Gladstone et al., 2013; Tao et al., 2011), because of an X-ray luminosity $\approx 10^{40}$ erg s*-1* (from the compact object) and an optical luminosity of $\approx 10^{38}$ erg s*-1* (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). In order to check whether there is 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.

5.3. Mass density of the secondary Roche lobe

The binary separation in CG X-1 is

[TABLE]

where G is the gravitational constant, $P$ the orbital period in unit of hr, $M_{1}$ and $M_{2}$ the masses of the compact object and donor star, respectively, in unit of M*⊙*. The radius of the secondary Roche lobe is

[TABLE]

(Eggleton, 1983), where $q\equiv M_{2}/M_{1}$ is the mass ratio. The average mass density inside the secondary Roche lobe (lower limit to the average density of the donor star) is

[TABLE]

(Eggleton, 1983).

In principle, the duration of the eclipse provides constraints on the mass ratio (coupled with the inclination angle); this 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 instead 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\lesssim q\lesssim 100$), the average density in the donor Roche lobe must be $0.55\lesssim\rho{\rm~{}g~{}cm}^{-3}\lesssim 2.2$ ($0.4\lesssim\rho/\rho_{\odot}\lesssim 1.6$). For a more likely range of $0.1\lesssim q\lesssim 10$, the density is $1\lesssim\rho~{}{\rm g~{}cm}^{-3}\lesssim 2.2$ ($0.7\lesssim\rho/\rho_{\odot}\lesssim 1.6$). Such 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.

5.4. Low-mass or high mass donor?

For a peak $L_{\rm X}$ $\approx 3\times 10^{40}$ erg s*-1*, the mass accretion rate must be a few 10*-6* $M_{\odot}$ yr*-1* even for radiatively efficient accretion. Considering also that in the super-critical 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*-5* M*⊙* $yr^{-1}$. 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 (van Haaften et al., 2019; Plotkin et al., 2014; David et al., 2005), and old globular clusters (Dage et al., 2018; Roberts et al., 2012; Shih et al., 2010; Maccarone et al., 2007); white dwarf donors in ultracompact systems are a plausible scenario for the ULX population in globular clusters (Dage et al., 2018; Steele et al., 2014). But ULXs in old stellar populations are much rarer than in young stellar environments (an order of magnitude rarer above $\approx$5 $\times 10^{39}$ erg s*-1*: Swartz et al. 2004; Plotkin et al. 2014). Moreover, we have already ruled out the possibility of a white dwarf donor, and 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: the donor may have started its life as an $\approx$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 $\approx$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, 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_{\rm X}\sim 10^{39}$–$10^{41}$ erg s*-1* (Wiktorowicz et al., 2015; Bogomazov, 2014), for timescales of a few 105 yr.

Second, WR stars have typical radii of only $\approx$1–2 $\times 10^{11}$ cm ($\approx$1.5–3 $R_{\odot}$) (Crowther, 2007), and can fit in a binary orbit with a period of a few hours. For a $M_{1}=10$ M*⊙* BH and a $M_{2}=20$ M*⊙* donor, and a period of 7.2 hr, $R_{L,2}\approx 2.6$ R*⊙* (large enough to fit a WR star), and $\rho\approx 1.6$ g cm${}^{-3}\approx 1.2~{}\rho_{\odot}$ (also consistent with a typical WR). Third, the short orbital period and eclipsing/dipping X-ray lightcurve 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; Carpano et al., 2007; Zdziarski et al., 2012; Laycock et al., 2015a; Ghosh et al., 2006; Maccarone et al., 2014; Esposito et al., 2013). The reason for this interpretation is that a WR provides a dense wind that strongly affects the X-ray lightcurve 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 lightcurve of CG X-1; we will revisit the comparison with other candidate WR X-ray binaries in Section 6.5.

6. DISCUSSION

6.1. What causes the asymmetric occultations?

We have shown (Section 3.2 and Figure 3) that the folded X-ray lightcurve (averaged over dozens of cycles) has an apparent eclipse (lasting for $\approx 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 showed (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’Aì et al., 2014; Díaz Trigo et al., 2006). In those 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, contrary to what we see in CG X-1.

Asymmetric eclipses are seen in some HMXBs (Falanga et al., 2015), such as the NS system Vela X-1. In that 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 such dips is partial covering of an extended X-ray emitting region by an opaque screen (Feng & Cui, 2002).

WR X-ray binaries are the subclass of HMXBs with the thicker stellar wind. Typical mass-loss rates are $\dot{M}\sim 10^{-4}$–$10^{-5}$ M*⊙* yr*-1*, with terminal velocities of $\sim$1000–3000 km $s^{-1}$ (Gräfener et al., 2017; Crowther, 2007). Thus, we expect that WR X-ray binaries are most likely to show the effect of variable wind absorption on X-ray lightcurves, 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; Crowther et al., 2010; Binder et al., 2015; Carpano et al., 2018) 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 lightcurve of NGC 300 X-1 also shows dips during a slow egress, consistent with variable partial covering of an extended X-ray emitting region (Binder et al., 2015); unlike CG X-1, the variable absorbing clumps are not Compton-thick (column densities of only $\approx$1023 cm*-2*), so 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 lightcurve 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 was suggested (Barnard et al., 2014) that in IC 10 X-1, the X-ray source is never completely occulted by the star, 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, in order to explain the asymmetric lightcurve 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 ($\approx$3 $\times 10^{40}$ erg s*-1*) is 300 times higher than the luminosity of IC 10 X-1 ($\approx$7 $\times 10^{37}$ erg s*-1*, Laycock et al. (2015a)). Therefore, the compact object itself in CG X-1 is expected to launch a much stronger wind (Ohsuga et al., 2005; King & Pounds, 2003; Pinto et al., 2016, 2017; Walton et al., 2016; Kosec et al., 2018) than in the other systems.

6.2. Colliding winds and bow shock

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, 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, 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., Pollock et al. 2005; Dougherty 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_{\rm{w1}}\sim v_{\rm{w2}}\sim 1000$ km s*-1*. 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 $\gtrsim$10${}^{-6}M_{\odot}$ yr*-1*. 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^{-6}M_{\odot}$ yr*-1*, in the super-Eddington disk outflow; this is consistent with our expectations that $\sim$10% of the mass lost by the WR is accreted by the BH101010A much lower fraction of stellar wind would be accreted by a NS, as the accretion cross section scales as $(M_{\rm NS}/M_{\rm BH})^{2}$.. The intersection of the contact discontinuity between the two colliding winds with the line of centres will be located at a distance $d_{1}$ from the compact object and $d_{2}$ from the donor star, with $d_{1}/d_{2}=(\dot{M}_{1}/\dot{M}_{2})^{1/2}\,(v_{\rm w1}/v_{\rm w2})^{1/2}$ and $d_{1}+d_{2}=a$ (Stevens et al., 1992). For example, for $\dot{M}_{1}\approx 0.1\dot{M}_{2}$, $d_{2}\approx(3/4)a$. The density $\rho_{\rm w2}$ of a spherically symmetric wind from the WR star at a distance $d_{2}$ from its centre is

[TABLE]

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_{\rm s2}\sim 10^{13}$ cm*-3*. The temperature of the shock-heated gas is $kT_{\rm s2}\approx(3/16)\,\mu\,m_{\rm H}\,v_{\rm w2}^{2}\approx 2$ keV for the assumed wind speed of $\approx$1000 km s*-1*, where $m_{\rm H}\approx 1.7\times 10^{-24}$ g is the mass of the hydrogen atom, and $\mu$ is the mean mass per particle of gas measured in unit of $m_{\rm 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_{\rm s1}\approx(v_{\rm w2}/v_{\rm w1})^{2}\,n_{\rm s2}\sim 10^{13}$ cm*-3*.

The cooling timescale for the shocked WR wind is

[TABLE]

(Dopita & Sutherland, 2003), where $\Lambda$ is the cooling function (Sutherland & Dopita, 1993). An analogous expression holds for the cooling timescale $t_{\rm cool,1}$ of the ULX wind. Thus, for a WR mass-loss rate of a few $10^{-5}M_{\odot}$ yr*-1*, we expect a cooling timescale of a few seconds for the shocked wind, and a few 10s of seconds for the shocked ULX wind, in a system with the 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_{\rm esc,1}=d_{1}/c_{s1}$ and $t_{\rm esc,2}=d_{2}/c_{s2}$, where $c_{\rm s1}$ and $c_{\rm s2}$ are the sound speeds in the two shocked layers. If $\chi\equiv t_{\rm cool}/t_{\rm esc}<1$, 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, $\chi_{2}\sim 10^{-3}$ for the shocked WR wind and $\chi_{1}\sim 10^{-2}$ for the shocked ULX wind. Thus, both shells are radiative and will collapse into geometrically thin, cold shells 111111As an aside, this is a crucial difference between the colliding winds in a compact system such as CG X-1 and those in WR-O star binaries. In the latter class of systems, because of their typically much larger binary separation, the wind density at the shock is lower, and both shells are usually adiabatic (Parkin & Pittard, 2008)..

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^{4}$ K, both for the WR and the ULX wind, and a temperature of $\approx$104 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 $\sim$10 km s*-1* and wind speeds $\sim$1000 km s*-1*, the compression factor is $\sim$104. Hence, we expect the two layers of radiatively cooled gas on both sides of the contact discontinuity to have densities of a few $10^{-8}$ g cm*-3*. From the solutions of Kenny & Taylor (2005), inserting typical parameters for this system, we obtain that each of the two shells has a thickness $\delta\sim 10^{-4}$ times the binary separation, that is a total thickness $\delta_{1}+\delta_{2}\sim 10^{8}$ cm. Thus, the surface density of the cold screen is $\sim$ a few g cm*-2*, corresponding to a Compton-thick column density of a few $10^{24}$ cm*-2*: 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 $\sim$108 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 108 cm, 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, we would not see sharp, individual dips in the lightcurve, but only smooth flux changes; conversely, if the clumps were much bigger than the X-ray emitting region, 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, b).

As we 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_{\rm K1}$ of the compact object in CG X-1 is

[TABLE]

For example, for $M_{1}=10M_{\odot}$ and $M_{2}=20M_{\odot}$, we obtain $v_{\rm K1}\approx 660$ km s*-1*; for $M_{1}=5M_{\odot}$ and $M_{2}=30M_{\odot}$, $v_{\rm K1}\approx 900$ km s*-1*; for $M_{1}=M_{2}=20M_{\odot}$, $v_{\rm K1}\approx 550$ km s*-1*. 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 is moving 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

[TABLE]

(Comeron & Kaper, 1998), where $\rho_{\rm a}$ is the density of the unshocked WR wind near the apex, at a distance of $\approx$ $a$ from the stellar centre. Using Equation (8), we can recast this distance as

[TABLE]

For our assumptions of $\dot{M}_{1}\sim 0.1\dot{M}_{2}$, and wind speeds of $\approx$1000 km s*-1*, this imples $R_{\rm 0}\approx 0.5a\approx 3R_{\odot}$. In terms of angular distance, the apex of the bow shock is $\approx$30*∘* in front of the BH.

As for 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 $\sigma$ of the cold gas can be approximated (Comeron & Kaper, 1998) as $\sigma\approx 1.25\rho_{\rm a}R_{\rm 0}\sim 1$ g cm*-2* for $\dot{M}_{2}\sim 10^{-5}M_{\odot}$ yr*-1*. This corresponds to an absorbing column density $\sim$1024 cm*-2*, 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 $\phi\approx 1.2$. The combined effect of true stellar eclipse, colliding winds and bow shock may keep the X-ray emission continuously occulted from phase $\phi\approx 0.90$ to $\phi\approx 1.1$, and partially covered from phase $\phi\approx 1.1$ to $\phi\approx 1.5$, as seen in some orbital cycles (Figure 2). Instead, we expect to see unocculted emission roughly from $\phi\approx 1.5$ to $\phi\approx 1.8$ when the bow shock is behind the X-ray source (Figures 3 and 11). As 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 considered for our model, ULXs also launch relativistic or semi-relativistic outflows at lower inclination angles from the polar axis (Walton et al., 2016; Pinto et al., 2016; Kosec et al., 2018; Pinto et al., 2017). Such outflows may create larger, parsec-scale shells of shocked circumstellar gas around the binary system. However, such shells will be Compton-thin and are not directly relevant to the phase-dependent occulting behaviour we are studying in this work; therefore, a study of the circumbinary ULX bubble is beyond the scope of this paper.

6.3. Anti-correlation between egress width and full luminosity

We showed (Figure 5) a trend of faster egress times in cycles with higher X-ray luminosities in the bright phase. We 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 corresponds also to stronger outflows from the compact object, and 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). 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 $\Delta\phi\approx 0.2$ at X-ray luminosities $\gtrsim$1040 erg s*-1*(Figure 5). We speculate that this minimum value is evidence of a true eclipse by the donor star, which 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.

6.4. Residual emission in the eclipse

One of the results of our spectral analysis of phase-resolved stacked spectra is that even during the true eclipse (or deepest occultation) at $\phi\approx 0.9$–1.1, we still detect residual flux from CG X-1, at a level of $\approx$1% of the bright-phase flux (L_{\rm X}$$\approx 2.5\times 10^{38} erg s*-1*). 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^{32}$ erg s*-1* (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 $\gtrsim$1000 km s*-1* for both the WR wind and the super-Eddington disk wind guarantee that the shocked gas reaches temperatures $\gtrsim$1 keV. In the simplest case of identical winds, each with mass loss rate $\dot{M}$ and speed $v_{\rm 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*-1* (Stevens et al., 1992), where the geometrical factor $f\approx 1/6$ parameterizes the fraction of wind kinetic luminosity perpendicular to the wind shock, $\dot{M}_{-5}$ is in units of $10^{-5}$ M*⊙* yr*-1*, and $v_{\rm w,3}$ is in units of 1000 km s*-1*. 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 w1}/v_{\rm w2})^{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^{38}$ erg s*-1*, values of $\dot{M}_{1}\gtrsim 10^{-6}$ M*⊙* yr*-1*, and $v_{\rm w1}\sim 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 $\sim$100 pc scale. On the other hand, before CG X-1, we had not identified any WR ULX, immersed in a circumstellar medium with densities $\sim$1012 cm*-3*. So, 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. Such thermal component has been invoked to explain the spectra of ultraluminous supersoft sources and the “soft excess” in ULX spectra (Tao et al., 2019; Soria & Kong, 2016; Urquhart & Soria, 2016; Middleton et al., 2015). It is generally consistent with a blackbody spectrum at temperatures $\approx$50–120 eV for supersoft sources, and $\approx$0.1–0.4 keV for ULXs, with a soft X-ray luminosity from $\sim$1038 erg s*-1* up to a few 1039 erg s*-1* (Zhou et al., 2019; Urquhart & Soria, 2016). If we fit the residual emission of CG X-1 with a blackbody model, we obtain a temperature of $\approx$0.5 keV (higher than seen in the soft excess of any other ULX) and an X-ray luminosity of $\approx$1038 erg s*-1* (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_{\rm{ph}}$) is related to the apparent blackbody radius ($r_{\rm{bb}}$) obtained from spectral fitting as $r_{\rm{ph}}\approx\sqrt{\tau_{\rm{es}}}\,r_{\rm{bb}}$ (Meier, 1982; Zhou et al., 2019), where $\tau_{\rm{es}}$ is the scattering optical depth at the photosphere ($\tau_{\rm{es}}\sim$ 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 $\tau_{\rm{es}}$ decreases to unity. The radius of this surface may be $\sim$100 times larger than the radius of the photosphere. In our case, $r_{\rm{bb}}\approx 130$ km, plausible values for $r_{\rm{ph}}$ are a few $10^{8}$ cm, and the radius of the last scattering sphere may be $\sim$1011 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^{2}_{\nu}>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 a fraction of X-ray photons emitted along the polar funnel and 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). It is usually a few percent of the direct luminosity out of eclipse, 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.

6.5. Binary evolution

To-date, we know of only seven compact (i.e., with a binary period $\lesssim 1$ d) 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 one (peak X-ray luminosity $L_{\rm X}\approx 3\times 10^{40}$ erg s*-1*), and the only ULX, clearly exceeding the Eddington limit for a $10$ M*⊙* BH.

We focus on compact WR binaries because they help us understand two key evolutionary features of binary systems, one in their past history and one in their future one. In the past, they must have come through a phase of dramatic shrinking of the binary separation from initial values of $\sim$100–1000 R*⊙* to $\lesssim$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 $\lesssim$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, 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). 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), provide crucial observational constraints to the predicted rate of LIGO/Virgo detections (Abbott et al., 2016a, b, 2017a, 2017b, 2017c). An analysis of the physical scenarios and possible steps between the X-ray binary stage and the final merger is beyond the scope of this work. Instead, we will briefly discuss the past and present evolutionary stages of these systems.

The common envelope phase required for the existence of systems such as CG X-1 is the one that occurs after the formation of the first compact object; this 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 star121212As the historical definition of a WR star is based on spectroscopic features rather than a precise range of physical parameters, the distinction between stripped-envelope helium stars and classical WR stars has always been somewhat blurred in the literature; helium stars more massive than $\approx$10 M*⊙* with WR-like spectra are usually simply called WR stars. See Götberg et al. (2018) for a comprehensive discussion of the relation between stripped helium-rich stars and WR stars across the mass range.. 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 $\sim$10*-5* M*⊙* yr*-1*, 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\approx$ 1.2–4 $\times 10^{-6}$ yr*-1*: Lommen et al. 2005). In this work (Section 3.2), we have shown a period increase at the 90% confidence limit also in CG X-1 ($\dot{P}/P=(10.2\pm 4.6)\times 10^{-7}$ yr*-1*).

Orbital widening inevitably leads to the donor star becoming detached from its Roche lobe, as 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 w2}^{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_{\rm 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 those two systems (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_{\rm X}\sim$ a few $10^{38}$ erg s*-1* (Koljonen et al., 2018, 2010; Hjalmarsdotter et al., 2008). 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_{1}^{2}$; 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 behaviour range from $\approx$2 M*⊙* (Zdziarski et al., 2013), 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 $\gtrsim$20 M*⊙* (Hjalmarsdotter et al., 2008).

This uncertainty leads us to the last issue we want 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 $\sim$1040 erg s*-1* (Bachetti et al., 2014; Israel et al., 2017a, b; Fürst et al., 2016), 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; Mushtukov2019; Pintore et al., 2017; Walton et al., 2018; Wiktorowicz et al., 2017; Wiktorowicz2019).

If CG X-1 is a wind-accreting rather than Roche lobe overflow system, one argument in favour 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}\sim$ a few $\times 10^{-5}$ M*⊙* yr*-1* (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 $\sim$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_{1}\sim M_{2}$) can this fraction approach $\sim$10%, perhaps with the aid of focused winds if the donor is closed to filling its Roche lobe. A WR + NS system is more likely to have $\dot{M}_{1}/-\dot{M}_{2}\lesssim 10^{-3}$, 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 super-critical 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 favour 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 Figs 1,2) that for a mass ratio $M_{2}/M_{1}\gtrsim 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_{2}\approx 40$–70 M*⊙); a 5-$\rm M_{\odot}$ 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 $\gtrsim$18 M⊙; more massive BHs ($M_{1}\gtrsim 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 $\gtrsim$1d (Pavlovskii et al., 2017; van den Heuvel et al., 2017; Bogomazov et al., 2018).

7. SUMMARY

We have studied the X-ray timing and spectral properties of the brightest X-ray source in the Circinus galaxy, confirming that it is a ULX with an X-ray luminosity $L_{\rm X}\sim$1040 erg s*-1*. The possibility of a foreground mCV is ruled out because the X-ray to optical flux ratio is $\gtrsim 100$, inconsistent with a CV.

We used a set of observations (two from ROSAT, 24 from Chandra, and five from XMM-Newton) that spans over 20 years. We phase-connected all the observations, thanks to the recurrent presence of an eclipse in each orbital cycle. This gave us a binary period $P=(25970.0\pm 0.1)$ s $\approx$7.2 hr (the most precise period measured in a ULX to date), and a period derivative $\dot{P}/P=(10.2\pm 4.6)\times 10^{-7}$ yr*-1* (error range at the 90% confidence level); $\dot{P}/P$ is larger than 0 at the 10$\sigma$ level. We showed that the X-ray profiles of each orbital cycle share the same general structure (fast ingress, eclipse, slow egress with dips), but each cycle has a unique shape, in terms of duration of the occultation phase, and dip structure. This is more complicated than a clean stellar eclipse, and suggests that the X-ray source is occulted by a complex structure of absorbers, located at approximately the same position with respect to the two binary components, but varying randomly over timescales shorter than a few hours.

As a first-order approximation, the X-ray spectrum is a power-law with photon index $\Gamma\approx 2$, seen through a cold local absorber (column density $\approx$ a few times $10^{21}$ cm*-2*). More detailed modelling of the spectra with highest signal to noise reveals two additional features. First, there is a spectral downturn above $\approx$ 4–5 keV. This curvature is one of the defining properties of ULXs, as opposed to sub-Eddington stellar-mass BHs. Second, there is a significant, additional absorption component from ionized gas, with column density $\sim$1022, more prominent during the egress phase than during the unocculted bright phase. This is typical of HMXBs, where we see the X-ray source through a thicker (and ionized) stellar wind when it passes behind the donor star.

The gradual transition from full occultation to fully bright phase does not correspond to a gradual decrease of the absorbing column density from Compton-thick to Compton-thin values. To a first approximation, we see an increase of the flux normalization during egress, without significant changes in the spectral shape; the increase is not monotonic, but is instead punctuated by irregular dips. This is consistent with a decrease of the covering fraction by an additional, Compton-thick partial-covering absorber. From the varying duration and structure of the occultation phases and dipping behaviour in each orbit, we suggest that the occulting material is made of opaque clouds with a characteristic size comparable to the size of the X-ray emitting region. From the orbital phase in which such occultations occur, we also infer that the absorbing material is mostly located between the two binary components and ahead of the accretor in its orbital trajectory. The duration of the occultation and dipping phases in each cycle is anti-correlated with the luminosity observed in the bright phase of that cycle. The peak luminosity varied between $\approx$4 $\times 10^{39}$ erg s*-1* and $\approx$3 $\times 10^{40}$ erg s*-1* over our 20-year coverage. Finally, we noticed residual X-ray emission in eclipse, at a luminosity $\sim$1038 erg s*-1*, with a softer spectrum than out of eclipse.

We proposed a simple model that can account for all these observational timing and spectral properties. We agree with the previous suggestions in the literature that CG X-1 is most likely a WR X-ray binary, as a WR star is the only type of massive star that can fit inside a Roche lobe radius with an inferred size $\lesssim$ 3 R*⊙, and at the same time provide enough mass transfer onto the compact object to generate X-ray luminosities $\sim$1040* erg s*-1*. In addition to the stellar wind, the super-Eddington accretor itself is expected to generate a powerful radiatively driven wind (this makes this system different from ordinary sub-Eddington HMXBs). Thus, we argued that the system should contain two regions of shocked wind: one between compact object and donor star (wind-wind collisions), and the other in front of the compact object along its orbital motion (bow shock into the dense medium created by the WR wind). Simple order-of-magnitude calculations show that the shocks are fully radiative; the shocked gas is so dense that it cools on timescales of a few seconds, collapses into a cold, Compton thick shell, and is expected to fragment into small clumps. Such clumps may be responsible for the partial covering and the irregular dips. The presence of a bow shock is the reason for the asymmetric location of the dips (mostly after rather than before the eclipse). The mass loss rate from the binary system is consistent with the observed increase in the period: $\dot{M}/M\sim\dot{P}/P\approx 10^{-6}$ yr*-1* (similar also to what is measured in the Galactic WR system Cyg X-3). The residual, softer emission in eclipse may come from the shocked wind, or, more likely, from X-ray photons scattered into our line of sight by the dense circumbinary medium.

We do not have direct observational evidence to determine whether the compact object is a BH or a NS. (BH identifications based only on super-Eddington luminosity arguments have been proven spectacularly unreliable in ULXs). However, we do favour the BH interpretation for two indirect arguments. First, the accretion cross section in wind-fed binaries scales as the square of the accretor mass, and a 1.5-$\rm M_{\odot}$ accretor is unlikely to intercept enough wind to generate steady luminosities $\gtrsim$1040 erg s*-1* for any plausible stellar wind model. Second, the short orbital period suggests that the system underwent a common envelope phase, which stripped the hydrogen envelope of the donor star and shrank the binary to a separation of a few solar radii; binary evolution models suggest that NSs are much less likely than stellar-mass BHs to survive common envelope.

The general significance of our study is that compact WR X-ray binaries (period $\lesssim$1 d) are a rare (only seven identified so far) and intriguing type of system, providing clues for crucial steps of binary evolution (see also the review of van den Heuvel 2019). They are also a step towards the formation of DCOs that can decay and merge via gravitational wave emission; constraining the formation rate of compact WR X-ray binaries helps predicting the LIGO/Virgo detection rate. On top of that, CG X-1 is the most luminous (in fact, the only persistent ULX) in this subclass, at least an order of magnitude above the others. Thus, its behaviour can be used for various tests of wind collisions, shocked shells, and super-critical accretion inflow/outflow models, in more extreme conditions than other HMXBs. Thirdly, CG X-1 is now the ULX with the most precise and accurate binary period, and even a tentative detection of binary period evolution. Further observational campaigns should be aimed at detecting or placing stronger upper limits on its point-like optical counterpart, ionized ULX bubble, and radio flux.

This work was supported by National Program on Key Research and Development Project (Grant No. 2016YFA0400800), and the National Natural Science Foundation of China (NSFC) through grants NSFC-11425313/11603035/11603038. We thank Ilaria Caiazzo, Paolo Esposito, JaeSub Hong, Chichuan Jin, Jiren Liu, Michela Mapelli, Manfred Pakull, Gavin Ramsay, Axel Schwobe, Lei Zhang for helpful discussions. Y.L. Q acknowledges the Harvard-Smithsonian Center for Astrophysics, and RS thanks Curtin University and The University of Sydney, for hospitality during part of this research. DJW acknowledges support from an STFC Ernest Rutherford fellowship. This work has made use of data obtained from the Chandra Data Archive, and software provided by the Chandra X-ray Center (CXC) in the application packages ciao. This work has made use of software obtained from the High Energy Astrophysics Science Archive Research Center (heasarc), a service of the Astrophysics Science Division at NASA/GSFC and of the Smithsonian Astrophysical Observatory’s High Energy Astrophysics Division. This work has also made use of the data from the XMM-Newton, an ESA science mission funded by ESA Member States and USA (NASA).