Quasi-periodic dipping in the ultraluminous X-ray source,
NGC 247 ULX-1

The PSD was estimated using the standard method of calculating the periodogram (e.g. Priestley 1981), with an ${\rm[rms/mean]}^{2}$ normalisation, then averaging over $M$ segments at each Fourier frequency. Non-overlapping neighbouring frequency bins are then averaged (see e.g. Vaughan et al. 2003). In order to pick out features over a broad range of frequencies we used a geometrically averaged frequency bins, each bin spanning a factor $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}1.3$ in frequency larger than the previous bin. Figure 2 shows the soft and hard band noise-subtracted PSDs in units of $\nu P(\nu)$ for both the dipping and non-dipping segments (the noise is estimated from standard formulae, e.g. Vaughan et al. 2003).

The dipping observations show clear peaks in power at $\sim 1\times 10^{-4}$ Hz and $\sim 2\times 10^{-4}$ Hz. A further broader peak can be seen at $\sim 8\times 10^{-4}$ Hz. At frequencies below $\sim 2\times 10^{-5}$ Hz the PSD is rolling over towards zero power. With the exception of the appearance of narrow peaks in power, the shape of the dip and no-dipping segments in the hard band looks similar. The non-dip PSD has a smooth power-law like frequency dependence, with a steeper slope in the hard band.

We study the correlated variations between the two bands using the Fourier coherence. This provides a measure of the linear correlation between two time series $x(t),y(t)$ as a function of Fourier frequency, in the range [0,1] (Vaughan & Nowak 1997). We estimated cross-spectrum, $C_{xy}(f)$, by first averaging the complex values over $M$ non-overlapping segments, then averaging over neighbouring frequency bins, in a similar manner to the PSD. This amplitude part of the cross-spectrum gives the coherence, $\gamma^{2}$, which is then corrected for Poisson noise, see e.g. Uttley et al. 2014 for more details.

The coherence is shown in Figure 2 for both the dipping and non-dipping light curves. A simple cubic spline has been fit to the data to illustrate the frequency dependence of the data. The coherence is high ($\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}0.9$) at frequencies where the soft and hard band PSDs peak, and drops at frequencies between these peaks in power. In contrast, a high coherence at low frequencies is observed in the non-dipping light curves, which falls off sharply to zero at $\sim 2\times 10^{-4}$ Hz, most likely caused by the onset of Poisson noise at high frequencies.

The phase part of the cross-spectrum provides an estimate of the frequency dependent delay between $x(t),y(t)$, which can be converted to a time delay in seconds (Vaughan & Nowak 1997). We estimated this for both the dipping and non-dipping segments and found the delay is consistent with zero lag at all Fourier frequencies.

To characterise the observed peaks in variability power in the dip segments we modelled the PSD using multiple Lorentzian functions;

The Akaike information criterion was used for model selection to compare models with additional complexity (Akaike 1974). The corrected form of this statistic takes account for the bias from the small sample size, which can be expressed as

We initially model the PSDs using a single broad Lorentzian component (model: Broadband), with all parameters free, giving $\Large\chi^{2}$/dof (degrees of freedom) $=1884.9/1594$. The $\Large\chi^{2}$ and AICc are given in Table 1 and the residuals of this model are shown in Figure 3. The next more complex model has the addition of a narrower Lorentzian ($l_{2}$) at $\sim 2\times 10^{-4}$ Hz (model: BB+1Lor). With ${\nu}_{\rm c}$ and $\sigma_{\rm l}$ tied between the two PSDs and $N_{\rm hard}$ and $N_{\rm soft}$ free, this improves the fit by $\Delta\Large\chi^{2}=22.7$ for 4 additional dof. Compared to the ‘Broadband’ model the ${\rm\Delta AICc}=14.6$, indicating a strong preference for the more complex model. Allowing any of the parameters to be untied between the two PSDs does not notably improve the model fit.

The addition of a Lorentzian ($l_{3}$) with ${\nu}_{\rm c}\sim 8\times 10^{-4}$ Hz (model BB+2Lor) reduces the $\Large\chi^{2}=14.9$ for 4 additional dof. Untying the parameters of the additional Lorentzian between the two PSDs does not significantly improve the fit. Compared to the previous model the ${\rm\Delta AICc}=10.9$, again suggesting the more complex model is strongly supported.

The most complex model we test involves the addition of a narrow Lorentzian ($l_{1}$) with ${\nu}_{\rm c}\sim 1\times 10^{-4}$ Hz (model: BB+3Lor). With all but the normalisation parameters tied across the two energy bands, the best fit model gives $\Large\chi^{2}$/dof $=1824.1/1584$, with 4 additional parameters compared to the previous model. The ${\rm\Delta AICc}=14.98$, again suggesting the more complex model is strongly supported. The low sampling resolution at low frequencies means the fit is largely insensitive to ${\nu}_{\rm c}$ of the broadband Lorentzian, so we fix this to zero (model: BB+3Lor($\nu_{\rm c_{\rm BB}}=0$)). This gives no notable change to the fit, but helps remove some parameter degeneracies between ${\nu}_{\rm c_{\rm BB}}$ and $\sigma_{\rm BB}$. Equally, ${\nu}_{\rm c,l_{1}}$ and $\sigma_{\rm l_{1}}$ (i.e. the Lorentzian at ${\nu}_{\rm c}\sim 1\times 10^{-4}$ Hz also display degeneracies, so we fix these at approximately the median MCMC values of $Q=\nu_{c}/\sigma_{l}=5$ (model: BB+3Lor($Q_{\rm l_{1}}=5$)). The parameters for this model are shown in Table 2. The addition of more Lorentzian components provides an improvement of ${\rm\Delta AICc}<2$, so we do not consider those models further.

For the no-dip segments, a single zero-centered Lorentzian provides an adequate description of the hard band PSD. The additional parameters are $\sigma_{\rm no-dip,hard}=4.55$${}_{-2.28}^{+1.25}$ $\times 10^{-5}$ Hz and $N_{\rm no-dip,hard}=0.29$${}_{-0.07}^{+0.06}$. Additional components do not significantly improve the AICc so we do not explore this further.

To check that no one particular observation is driving the peaks in the average PSD we estimate the unbinned periodogram. Figure 4 shows the periodogram for the four dipping observations. It is well known that the periodogram estimate at each frequency follows a ${\Large\chi}^{2}_{\nu=2}$ distribution so large scatter around the true underlying power spectrum is expected (e.g. Priestley 1981). We therefore also show a Gaussian smoothed with periodogram over $n=3$ neighbouring frequency bins. Peaks in power at $\sim 1\times 10^{-4}$ Hz and $\sim 2\times 10^{-4}$ Hz are found in each periodogram. Another peak or bend can be seen at $\sim 0.4\times 10^{-4}$ Hz in all but rev $3663$, which shows power at lower frequencies.

We refit these periodograms using the maximum likelihood methods (MLE) for unbinned periodograms outlined in Vaughan (2010) and Barret & Vaughan (2012). The power spectral modelling results are consistent with those found in Section 3.3, so we do not show those here. For rev $3676$ the Lorenztzian centroids $l_{1}$ and $l_{2}$ appear to shift to higher frequencies, however modelling reveals these to be consistent within error bars with the other revs. The only notable difference is the increase in power at low frequencies in rev $3663$, which requires a zero-centered Lorentzian with width $\sigma_{l_{\rm BB}}\sim 1\times 10^{-5}$ Hz. Removing rev $3663$ from the combined PSD modelling in Section 3.3 does not significantly alter the results, so we choose to leave it in to increase the number of estimates in each frequency bin.

In order to study what aspect of the dips is causing the peaks in the power spectra, we modelled the abrupt light curve changes in the dip segments using change point analysis. The R package bcp¹¹1https://cran.r-project.org/web/packages/bcp/ was used which implements the Bayesian change point method of Wang & Emerson 2015. This partitions the light curve into constant mean ‘blocks’ between change points. MCMC sampling is used to model the posterior probability at each change point, $p_{\rm ch}$, defined between [0,1]. Given the sharp changes in the light curves when dipping occurs, this method is particularly well suited to the data.

The Poisson noise starts to dominate at $\sim 2\times 10^{-3}$ Hz so a time bin of $dt=500$ s was used in order to smooth over any noise variations. A posterior probability $p_{\rm ch}>0.4$ was used to identify significant change point locations. From these we determine the dip onset time $t_{\rm 0}$, dip duration $t_{\rm d}$, and dip midpoint, $t_{m}=(t_{\rm 0}+t_{\rm d})/2$. The period $P$ is then taken as the time between successive dip midpoints.

Figure 5 shows the distribution of $t_{\rm d}$ and $P$ for all four dipping observations. A Gaussian kernel density estimate is overlaid on the observed distributions. The distribution of $t_{\rm d}$ peaks at $\sim 3$ ks and falls off approximately log-normally with increasing dip duration, i.e. there is a much larger probability for shorter dips. The distribution of $P$ is multi-modal with peaks at $\sim 5$ and $10$ ks, with a smaller probability of periods at $\sim 17$ ks. A similar distribution is observed if we represent the $P$ as the time between successive $t_{\rm 0}$, suggesting no relationship between the duration of the dips and the dip period.

We study the energy dependence of the variability at a given timescale using frequency resolved covariance spectra (Wilkinson & Uttley 2009). Using a high S/N reference band, the correlated variability is picked out in a given comparison band, thus improving the S/N of the variability spectrum compared to traditional rms spectra. We compute the covariance spectra in the Fourier domain (Uttley et al. 2011; Cassatella et al. 2012; Uttley et al. 2014). The $0.3-4.0$ keV band is used as the broad reference band, but the choice of this band does not affect the shape of the resultant covariance spectra.

We estimate covariance spectra in three frequency bands: $f_{1}=2.0-4.0\times 10^{-5}$ Hz, $f_{2}=0.9-1.1\times 10^{-4}$ Hz and $f_{3}=1.9-2.4\times 10^{-4}$ Hz. These timescales probe the long timescale spectral variability as well as the spectral variability associated with the two prominent dipping timescales. Figure 6 shows the resulting variability spectra. The top panels shows these data in absolute units (${\rm cts}~{\rm s}^{-1}~{\rm keV}^{-1}$) whilst the lower panels show the variability spectra in fractional ($\%$) units. In the dipping segments, the long timescale spectral variability is approximately flat, with a $10\%$ rms amplitude. The dipping timescales, $f_{2}$ and $f_{3}$, show an approximately linear dependence with energy, with larger variability amplitude at higher energies. This is true of all frequency bands up to higher frequencies, but we do not show those here.

For the dipping segments, as an additional check, we calculate the variability amplitude of each of the components in the preferred model (BB+3Lor $Q_{\rm l_{3}}=5$) by integrating the rms in each Lorentzian component. This shows the same trend in energy space as the frequency dependent covariance spectra. For the no-dip segments, the fractional variability amplitude is almost below $1$ keV, with a sharp increase in rms for energies above this.

The light curves of NGC 247 ULX-1 display deep and quasi-periodic dips imprinted on top of stochastic noise variations. Modelling the PSD with multiple Lorentzians reveals three preferred timescales for the dipping behavior. The peak in power of a Lorentzian occurs at frequency $\nu_{\rm max}=\sqrt{\nu_{c}^{2}+\sigma_{l}^{2}}$ (see e.g. Belloni et al. 2002). For the preferred model, this gives $\nu_{\rm max,BB}\sim 0.4\times 10^{-4}$ Hz, $\nu_{\rm max,l_{1}}\sim 0.9\times 10^{-4}$ Hz, $\nu_{\rm max,l_{2}}\sim 2\times 10^{-4}$ Hz and $\nu_{\rm max,l_{3}}\sim 9\times 10^{-4}$ Hz. The coherence between the hard and soft bands is high ($\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}0.9$) at the Lorentzian peaks and drops where the Lorentzians overlap, providing further evidence for multiple dipping timescales.

Time series change point analysis reveals the dip periods, $P$, cluster around $\sim 5$ and $10$ ks. This distribution of $P$ is consistent with the peak timescales observed in the power spectra (see Figure 3). The Lorentzian $Q$ values tell us that there is some scatter around a preferred periodicity, as is observed in the change point period analysis. The dip durations $t_{d}$ preferentially take on values between $1-5$ ks, which will add some additional broadening to the measured Lorentzian $Q$ values. There are several broad dips with $t_{d}>15$ ks, which clearly stand out from the other measured dip durations in Figure 5. This could be suggesting an alternative mechanism or timescale for these longer duration dips.

The change point analysis does not find any periods with $P<1000$ s suggesting this analysis is either insensitive to this dips, e.g. they could be of lower amplitude and therefore harder to detect, in agreement with the lower amplitude of the Lorentzian $l_{3}$. Alternatively, this high frequency component is related to an increase in the broadband noise continuum, which appears when the source enters the higher-flux / dipping epochs.

The time averaged and frequency resolved covariance spectra in Figure 6 show that the variations are larger at higher energies in both the dip and no-dip epochs. The exception being the $f_{1}=2.0-4.0\times 10^{-5}$ Hz frequency band in the dip epochs. Here, the fractional covariance is independent of energy. The coherence between the soft and hard bands (see Figure 2) also drops to $\sim 0.6$ at these frequencies. This suggests that the long timescale variability is different to the dip generating mechanism at shorter timescales.

Assuming a black hole with $M_{\rm BH}=10M_{\odot}$, the peak timescales $\nu_{\rm max}\sim 0.5-9\times 10^{-4}$ Hz correspond to an orbital radius (${}^{3}\sqrt{GMt^{2}/4\pi^{2}}$) of $10^{11}$ cm, or $2-10\times 10^{4}$ $R_{G}$ (where $R_{G}=GM/c^{2}$). A similar value is obtained for the case of a NS compact object. In both cases, this region is very distant from the compact object, suggesting an origin in the outer accretion flow. Two likely scenarios for the dips are warping of an accretion disc seen close to edge on, periodically obscuring the central source each orbit, or, the launching of a wind following changes in local $\dot{M}$. Alternatively, a combination of both processes could be causing the dips at different timescales.

The radius associated with Compton scattering through a medium around the spherisation radius can be expressed as $R_{\rm C}=\sqrt{ct_{\rm esc}/\sigma_{\rm T}n_{\rm e}}$ (where $c$ is the speed of light, $\sigma_{\rm T}$ is the Thomson scattering cross-section, $n_{\rm e}$ is the electron density, and $t_{\rm esc}$ is the timescale for photons to escape) would be of the order of $10^{9}$ cm, assuming a disc-like $n_{\rm e}$ density ($10^{20}$ cm^-3), while such radius would be comparable ($10^{11-12}$ cm) if we adopt a wind-like density ($10^{14}$ cm^-3). On the other hand, the wind radius, $R_{\rm wind}=L_{\rm ion}/n_{\rm e}\xi$ (where $L_{\rm ion}$ is the ionising luminosity and $\xi$ the ionisation parameter), that we estimate from the best-fit parameters of the XMM-Newton RGS-EPIC joint fit is also around $10^{11-13}$ cm (see Pinto et al. 2021). The lower density assumption is favoured as this agrees with both the outflow and dip timescale properties.

The mid plane of the accretion disc, i.e. the spherisation radius where the disc height is comparable to the distance from the accretor, is instead $R_{\rm sph}=5/3\dot{M}R_{\rm in}\sim 100$ $R_{G}$ (assuming a standard $R_{in}=6$ $R_{G}$ and $\dot{M}\sim 10$ $\dot{M}_{\rm Edd}$, which would be sufficient to achieve $L_{X}\sim 5\times 10^{39}\hbox{${\rm erg}~{\rm s}^{-1}$}$ for $M_{\rm BH}=10M_{\odot}$). This means that the timescales of the variability detected in NGC 247 ULX-1 are likely associated with the outer disc and, most likely, with the region where the wind becomes optically thin (e.g. $n_{\rm e}\lesssim 10^{14}$ cm^-3) and rather than with the inner accretion flow or the optically thick wind ($10^{20}$ cm^-3) which would be expected near the spherisation radius. For a comparison, the wind photosphere in the Galactic super-Eddington source SS 433 is $\sim 10^{12}$ cm (see, e.g. Fabrika 2004), which further indicates that NGC 247 ULX-1 quasi-periodic dips might be associated with some structure present in the outer disc where the photosphere becomes optically thin.

The presence of dips and their long timescales argue in favour of a moderate-to-high inclination in this source, which agrees with the scenario in which soft ULXs are primarily seen along the edge of the wind. The dips seen in both BH and NS X-ray Binaries (XRBs) are believed to be caused by periodic obscuration of the central X-ray source by structure located in the outer regions of a disc. These dips are observed in sources at high inclination ($i>60$ deg, e.g. Frank et al. 1987). For NGC 247 ULX-1, Yao & Feng (2019) derived an $i>67$ deg via modeling the UV/optical spectral energy distribution (SED) with an irradiation model, providing further evidence for a high inclination viewing angle in this source. Additionally, winds observed in BH or NS X-ray binaries (XRBs) have been mostly observed in sources in which the disc is inclined at a large angle to the line of sight (e.g. Ponti et al. 2012; Díaz Trigo & Boirin 2016).

Given the difficulty to access the inner regions of the accretion flow, mostly obscured by the scale height of the disc at the spherisation radius and the wind, we are not surprised that we did not detect the time lags between the soft and hard band that have been found in other, harder, ULXs. At these angles, the photons would undergo loads of Compton scattering, losing the original information.

If the dips are caused by a wind structure, what is not clear is whether the short-term variations in optical depth are just stochastic variability (at a given $\dot{M}$) due to the clumpy nature of the wind. Alternatively, each variation corresponds to a change in the underlying $\dot{M}$ variability in the local accretion rate where the wind is launched. This variability in the wind launching would thereby cause these quasi periodic dips, rather than being pre-existing structure with the appropriate scale height. This idea is supported by the onset of dipping behaviour in the higher flux observations. In fact, the bolometric luminosity of SSUL sources slightly increases when the blackbody temperature decreases and the blackbody radius increase. This behaviour is common to several other SSUL (e.g. Feng et al. 2016; Urquhart & Soria 2016a).

The broadband spectral model decomposition in D’Aì et al, in prep reveals two broad thermal components. These are interpreted as emission from an extended, optically thick, cold photosphere, and a hotter emission component, originating from regions closer to the compact object. In the ‘out of dip’ epochs (i.e. when a dip is not occurring), the short term variations are predominantly driven by changes in the flux of the hotter component. During the dips, both broad spectral components show significant decrease in their flux, particularly the hotter component. One scenario, where the dipping mechanism is caused by occultation further out than the typical radii where both components originate, would be consistent with what we infer from the timing properties in the present paper.

An alternative explanation for the dips is that of intrinsic luminosity variations due to the propeller effect switching on-and-off accretion in a magnetised NS (see e.g. Revnivtsev & Mereghetti 2015 for a review). This was suggested as the origin of large amplitude variations in the X-ray light curve of the ULX M 82 X-1 by Tsygankov et al. (2016), however these occur on timescales of $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}$ weeks. Interestingly, the ‘hiccuping’ X-ray light curve of the accreting millisecond X-ray pulsar, IGR J18245-2452, displays a remarkable resemblance to those in NGC 247 ULX-1 (Ferrigno et al. 2014), However, the dips in this source show complex spectral hardening, whereas in NGC 247 ULX-1 the dips are spectrally harder (see Pinto et al. 2021 Figure 2). The detection of pulsations in NGC 247 ULX-1 would provide more evidence for this scenario, which is scope for future work.

For NGC 55, Stobbart et al. (2004) report several dips in the XMM-Newton light curve. These are fewer in number, but have a similar dip duration and separation between consecutive dips. The authors argued the dips were likely due to material at $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}2\times 10^{11}$ cm from the central object (assumed to be a BH). At this radius the orbital timescale is $\approx 10^{4}$ s, consistent with our PSD peak timescales. This is in agreement with the idea that a disc warp or extended wind structure is causing the dips each orbit in NGC 247 ULX-1.

The two broad dips with $t_{d}\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}15$ ks in observations rev $3663$ and rev $3677$ are similar to the broad deep dips seen in M 51 ULX-1 (Urquhart & Soria 2016b). These authors argued that these dips are caused by eclipsing events and placed limits on the orbital period. CG X-1, a ULX in the galaxy Circinus, has also shown multiple eclipsing events (Qiu et al. 2019), with the authors finding an orbital period of $\sim 7$ hrs. It is possible that we are seeing a similar eclipses in NGC 247 ULX-1, meaning the separation of $\sim 19$ days for these deep dips may correspond to some integer multiple of the orbital period. If this is the case, then the dipping activity observed in NGC 247 ULX-1 may be triggered by the companion at some point in the orbital phase. This remains scope for future work.

A tentative $\sim 2$ hr ($\sim 1\times 10^{-4}$ Hz) quasi-periodic oscillation (QPO) was claimed in NGC 628 (Liu et al. 2005). The light curve of this source displays a similar dipping process to NGC 247 ULX-1. The authors scaled this frequency to QPOs typically seen in black hole X-ray binaries (XRBs; see e.g. Remillard & McClintock 2006) and inferred a black hole mass $M_{\rm BH}\sim 10^{3}M_{\odot}$. However, the light curve properties of the ULX dips are strikingly different to QPOs seen in XRBs and active galaxies (e.g. Ingram & Motta 2020, Alston et al. 2014; Alston et al. 2015, 2016). If the dips are caused by a similar process to NGC 247, then this would argue in favour of a stellar mass compact object in NGC 628.