# Fast spectroscopic imaging using extreme ultraviolet interferometry

Hannah C. Strauch  1st Institute of Physics, University of Göttingen, Germany
Fengling Zhang  Advanced Research Center for Nanolithography ARCNL, Amsterdam,
the Netherlands  Stefan Mathias  1st Institute of Physics, University of
Göttingen, Germany  Thorsten Hohage  Institute for Numerical and Applied
Mathematics, University of Göttingen, Germany  Stefan Witte  Advanced Research
Center for Nanolithography ARCNL, Amsterdam, the Netherlands  G. S. Matthijs
Jansen  1st Institute of Physics, University of Göttingen, Germany
<EMAIL_ADDRESS>
###### Abstract

Extreme ultraviolet pulses as generated by high harmonic generation (HHG) are
a powerful tool for both time-resolved spectroscopy and coherent diffractive
imaging. However, the integration of spectroscopy and microscopy to harness
the unique broadband spectra provided by HHG is hardly explored due to the
challenge to decouple spectroscopic and microscopic information. Here, we
present an interferometric approach to this problem that combines Fourier
transform spectroscopy (FTS) with Fourier transform holography (FTH). This is
made possible by the generation of phase-locked pulses using a pair of HHG
sources. Crucially, in our geometry the number of interferometric measurements
required is at most equal to the number of high-harmonics in the illumination,
and can be further reduced by incorporating prior knowledge about the
structure of the FTH sample. Compared to conventional FTS, this approach
achieves over an order of magnitude increase in acquisition speed for full
spectro-microscopic data, and furthermore allows diffraction-limited
computational imaging.

††journal: opticajournal

##  1 Introduction

Coherent extreme ultraviolet light produced by high-harmonic generation [1]
provides a unique opportunity for the study of ultrafast dynamics in condensed
matter at the nanoscale. On the one hand, the short wavelength yields a
favourable Abbe resolution limit, and extreme ultraviolet (EUV) coherent
diffractive imaging nowadays enables imaging at nanometer-scale resolution [2,
3, 4]. On the other hand, the wide range of photon energies and ultrashort,
attosecond (as) to femtosecond (fs) pulse durations of EUV light pulses from
high-harmonic generation (HHG) provide access to a wide range of elemental
absorption edges. Here, time-resolved spectroscopy yields a sensitive element-
resolved probe of dynamics of the electrons, the spins and the lattice in
optically excited matter [5, 6, 7, 8, 9, 10, 11]. The combination of both
these approaches, namely nanometer-scale microscopic time-resolved EUV
spectroscopy, is highly appealing, as it would enable the study of complex
dynamics in nanoscale structures ranging from naturally-inhomogeneous quantum
materials to fabricated heterostructures. However, such multi-dimensional
measurements remain out of reach for HHG light sources as is, and can so far
only be performed at significant experimental cost at accelerator-based EUV
light sources [12, 13].

Holography plays an important role in the development of advanced EUV
microscopy [14, 15, 16, 17, 18, 19]. For example, time-resolved EUV microscopy
is commonly performed by combining Fourier-transform holography (FTH) with
numerical phase retrieval [20, 13]. In FTH, the diffracted wave from the
sample is interfered with a reference wave that is generated by a point-like
structure (usually a pinhole in transmission geometry). The resulting fringe
pattern in the far-field diffraction pattern allows for direct image
reconstruction by a single Fourier transform. The single-shot nature of the
FTH measurement is critical for the implementation of time-resolved studies,
as it implies that only the pump-probe delay needs to be scanned during an
experiment. In contrast, ptychography is a powerful imaging method that
enables diffraction-limited resolution and can handle experimental challenges
such as partial coherence, multi-spectral illumination and structured
illumination profiles [21, 22, 2]. To achieve this, ptychography relies on
spatial scanning of the sample in overlapping steps such that the phase
retrieval problem is sufficiently constrained. The need for overlap is
significantly increased for multi-spectral illumination [22, 23, 2]. A full
time- and spectrum-resolved ptychography measurement would therefore place
significant requirements on the brightness and stability of the EUV light
source. Thus, there is a clear need for a fast spectromicroscopy technique
that reduces the data collection time for scanning in the spatial and spectral
domains.

In this Article, we demonstrate Fourier-transform spectroscopic holography
(FTSH) (Fig. 1), which combines FTH with an interferometric Fourier-transform
spectroscopy measurement. For an HHG spectrum with $N_{\lambda}$ harmonics,
FTSH requires at most $N_{\lambda}$ measurements and allows for further
reduction by incorporating prior knowledge of the sample structure. In the
experiment, we achieve flux-efficient FTSH by focusing the output of two
phase-locked HHG sources [24] to a conventional FTH sample, such that the
sample and reference structures are illuminated by separately controlled EUV
pulses. In this geometry, we exploit the intrinsic coupling between spectral
and spatial information in the diffraction pattern, and show that this enables
a strong reduction of the FTS sampling requirements. For our HHG spectrum
(covering up to 11 harmonics from 80 to 29 nm, see Fig. 2b) and sample
geometry, we find that only 6 phase steps are necessary to recover all
spectral components completely.

Figure 1: Overview of the Fourier transform spectroscopic holography (FTSH)
workflow. (a) Two phase-locked polychromatic EUV beams from a pair of HHG
sources illuminate the sample and reference structures. (b) The diffraction
pattern is detected by a camera in the far-field. (c) As in conventional FTH,
a single Fourier transform is used to analyze the interference pattern,
however in this case the holograms of the individual photon energies are
superimposed. (d) By time-shifting the reference wave, Fourier transform
spectroscopy allows to measure the spectrum at each pixel. The necessary
sampling for this measurement is significantly reduced when the multi-
wavelength nature of the hologram (e) is considered, as the different spectral
components are only partially overlapping. Thus, an HHG spectrum with
$N_{\lambda}$ harmonics requires at most $N_{\lambda}$ and typically less
measurements to fully recover a spectrally-resolved image (f).

##  2 Theory

Fourier-transform spectroscopic holography (Figure 1) can be understood as an
extension of FTH towards multiple wavelengths. For a single wavelength
$\lambda_{i}$, the far-field diffraction pattern of an FTH sample can be
expressed as
$I_{\lambda_{i}}(k)=|\tilde{p}_{i}|^{2}+|\tilde{r}_{i}|^{2}+\tilde{p}_{i}\tilde{r}_{i}^{*}+\tilde{p}_{i}^{*}\tilde{r}_{i},$
where $\tilde{p}$ and $\tilde{r}$ are the electric fields at the detection
plane due to the probed sample and reference, respectively, while $k$ is the
spatial frequency of the scattered wave. The fields $\tilde{p}$ and
$\tilde{r}$ can be related to the fields $p$ and $r$ at the sample plane by
optical propagation. For Fraunhofer diffraction, the propagation is given by a
scaled Fourier transform, where the scale (i.e., the extent of the diffraction
pattern on the camera) is proportional to the wavelength of illumination.
Consequently, a Fourier transform ($\mathcal{F}$)

$\mathcal{F}\\{I_{\lambda_{i}}\\}(x)=\ p_{i}\ast p_{i}^{*}+r_{i}\ast
r_{i}^{*}+p_{i}\ast r_{i}^{*}+p_{i}^{*}\ast r_{i}$ (1)

provides direct access to the holograms ($p_{i}\ast r_{i}^{*}$ and
$p_{i}^{*}\ast r_{i}$) in the real-space coordinates $x$. It is well known
that the resolution of these holograms is limited by the size of the reference
$r$, and furthermore that the hologram appears separated from the
autocorrelation terms ($p_{i}\ast p_{i}^{*}$ and $r_{i}\ast r_{i}^{*}$) at a
distance $\Delta$ from the center that is equal to to the separation of $p$
and $r$. In terms of pixels of the reconstructed image, the separation becomes
$\Delta_{\rm px}=2\Delta D/\lambda z$, where $z$ is the distance and $D$ the
size of the detector. Notably, the numerical hologram separation $\Delta_{\rm
px}$ is inversely proportional to the wavelength as a consequence of the
wavelength scaling of the Fraunhofer propagation operator.

While the precise position of the hologram in the frame is not important in
monochromatic FTH, it becomes relevant for spectrally-resolved applications.
If an FTH hologram is recorded using polychromatic radiation, a spatial
Fourier transform of the data yields an image where the individual holograms
are shifted according to the illuminating wavelength. This effect has been
harnessed to perform single-shot spectrally resolved FTH using the individual
high harmonics of a HHG light source [25]. However, this requires a specific
experimental geometry: in order to fully separate the individual holograms,
the sample width needs to be reduced significantly (or equivalently the
$\Delta$ must be increased) by a factor that depends on the wavelength
difference of neighbouring high harmonics. This drastically reduces the
achievable field of view and limits the spectral resolution, altogether
preventing the application to continuous (attosecond) EUV spectra. Therefore,
another method is necessary to resolve the individual spectral components.

Fourier transform spectroscopy (FTS) is a powerful method that provides such a
capability: analogous to FTH, interference of two pulses delayed in the time
domain gives access to the spectrum by an inverse Fourier transform. FTS-based
methods are already commonly used to perform hyperspectral imaging in the
visible and infrared ranges (e.g., [26, 27]). In the last decade, FTS at
extreme ultraviolet wavelengths has been enabled by the generation of phase-
locked EUV pulse pairs by phase-locked laser pulses [24, 28]. Integrated in a
coherent diffractive imaging experiment, spatially-resolved FTS yields
monochromatic diffraction patterns that can be used for numerical image
reconstruction [29]. In order to satisfy the Nyquist-Shannon sampling theorem
and to accurately measure the complete spectrum, however, typically a few
hundred interferometric diffraction patterns at different delays must be
measured.

Here, we show that the individual limitations of FTH and FTS can be lifted
when they are combined in a HHG-based interferometric measurement (cf. Fig.
1). In particular, we propose to combine the FTSH measurement with prior
knowledge of the experimental parameters. The generated wavelengths of a HHG
source are generally known in advance to a good accuracy, e.g., from a prior
spectroscopy measurement or from a measurement of the fundamental laser
parameters. This knowledge can already significantly reduce the sampling
requirements, but a further reduction is possible: Combined with knowledge of
the FTH mask structure, it is possible to calculate which spectral components
can contribute to each pixel in the FTH hologram. For comb-like HHG spectra
spanning an octave or more, it is generally found that only a small number of
harmonics contribute to each pixel in the multi-wavelength hologram.

It is instructive to consider the following example: an FTH sample containing
a reference structure with negligible width and a sample with width $W$ (along
the sample-reference direction) with a sample-reference separation $\Delta=2W$
is illuminated by a HHG spectrum containing the odd harmonics 13 to 29 of an
1030 nm driving laser. In this case, the hologram of the 29th harmonic has
significant overlap with the 27th, successively less overlap with the 25th,
23rd and 21st harmonic holograms, and no overlap with holograms of the 19th
and lower high harmonics. In fact, no more than 5 high harmonics overlap at
any point in the multi-wavelength hologram. This suggests that the sampling
requirements in an FTS measurement can be reduced.

In FTSH, the reference beam is delayed by a time $\tau$ with respect to the
probe beam. This leads to a measured intensity given by
$I(k,\tau)=\sum_{i}^{N_{\lambda}}|\tilde{p}_{i}+\tilde{r}_{i}e^{-i2\pi
c\tau/\lambda_{i}}|$. By applying a Fourier transform, we can express the
measured multi-wavelength hologram $M(x^{\prime},T)$ as

$\begin{split}M(x^{\prime},\tau)=\sum_{i}^{N_{\lambda}}&\ p_{i}\ast
p_{i}^{*}+r_{i}\ast r_{i}^{*}+(p_{i}\ast r_{i}^{*})e^{-i2\pi
c\tau/\lambda_{i}}+c.c.,\end{split}$ (2)

where $c.c.$ indicates the complex conjugate and $c$ is the speed of light.
The coordinate $x^{\prime}$ is Fourier conjugate to the pixel coordinates of
the camera, and can be related to the real-space coordinate by considering the
wavelength-dependent resolution and shift $\Delta_{\rm px}$. It deserves
emphasis that the _complex-valued_ data $\sum_{i}^{N_{\lambda}}(p_{i}\ast
r_{i}^{*})e^{-i2\pi c\tau/\lambda_{i}}$ can be directly retrieved from the
diffraction pattern measured at a single delay $\tau$ by spatially isolating
the multi-wavelength hologram (see Fig. 1c, e).

The spectroscopic reconstruction problem can therefore be posed as follows:
Given $M(x^{\prime},\tau_{j})$ for $j\in[1,N_{\tau}]$, find the complex-valued
$p_{i}^{*}\ast r_{i}$ for each $i\in[1,N_{\lambda}]$ at each $x^{\prime}$. For
suitably chosen delays $\tau_{j}$, this problem is well posed, and can be
solved by matrix inversion or least-squares methods. Ignoring prior knowledge
of the shape of the holography sample there are at most $N_{\lambda}+1$
unknown complex values for each $x^{\prime}$, indicating that at most
$N_{\tau}=N_{\lambda}+1$ measurements are required. This is reduced by 1 for
sufficiently large $\Delta$, since the autocorrelation terms can be filtered
out spatially. A more dramatic reduction is achieved if the partial overlap of
the monochromatic holograms is considered, as explained by the previous
example.

Figure 2: Fast Fourier-transform spectroscopic holography using minimal
sampling. a) Scanning electron microscopy image of the FTH sample, which was
made by focused ion beam milling of a gold-coated silicon nitride membrane. b)
Exemplary amplitude $|M(x^{\prime},T)|$ of the multi-wavelength hologram,
showing that the different spectral components are partially overlapping. c)
Extracted spectral intensity from the high resolution FTS scan, a dataset with
$N_{\tau}=11$ shots and $N_{\tau}=6$ shots. The fundamental wavelength is 1030
nm. For comparison, the data was normalized to the total observed intensity.
The small differences between the different datasets can be attributed to
intensity drift and shot-to-shot variation in the HHG light source. d)
Spectroscopic holography results from the $N_{\tau}=11$ dataset for the 6
brightest high harmonics. e) The corresponding spectroscopic holography
results from the $N_{\tau}=6$ dataset.

##  3 Spectroscopic holography

In order to record FTSH data at EUV wavelengths, we employ an ultrastable
birefringent common-path interferometer [30, 24] to split the output of a 1
kHz, 35 fs, 1030 nm laser (Light Conversion Pharos, compressed by cascaded
nonlinear compression in 800 mbar argon [31]) into phase-locked pulses with a
controllable delay. By tilting one wedge of the common-path interferometer,
these pulses are focused to two separate spots in a Krypton gas cell, where
they generate two phase-locked EUV beams. After a 200 nm Al filter to block
the fundamental laser light, the high-harmonic spectrum of both pulses spans
from 80 to 29 nm. At 94 cm after the HHG, we use a broadband curved multilayer
mirror ($f=25$ cm) to image the pulse pair onto the sample region. In the
focus plane, the two EUV beams are separated by roughly 70 µm, and the
diameter of the individual beams is roughly 25 µm. To facilitate a more
homogeneous illumination of the sample, we place the FTH sample $\approx$5 mm
downstream of the focus. The EUV camera (Andor IkonL, 13.5 micron pixel size)
is placed 10 cm behind the sample.

In this geometry, we now record full FTSH data. Analysis of a high-resolution
FTS scan with a time-step of 23 as over a range of 11.6 fs confirms that the
HHG spectrum consists of 11 high harmonics with frequencies that match to the
odd harmonics 13 to 33. The oscillation period of these harmonics ranges from
262 as to 103 as. In the following, this dataset will serve as a reference and
allows to benchmark the analysis results when only using few-delay subsets of
the full data.

Next, we therefore subsample the full FTS scan in order to demonstrate the
minimum sampling requirements for full spectroscopic image reconstruction.
Specifically, we now consider interferometric diffraction patterns at
$N_{\tau}=N_{\lambda}=11$ equidistant delays in the range 0 - 1.7 fs (step
size 190 as), which corresponds to a half cycle of the driving laser pulse. By
the Nyquist theorem, the spectral resolution of the Fourier transform of this
data exactly matches with the spacing of the high harmonics, namely twice the
fundamental frequency. However, the sampling frequency is much lower than the
highest frequency in the data, and we emphasize that the spectral amplitudes
cannot be reconstructed using a Fourier transform. Instead, we solve (2) using
the Newton conjugate-gradient method implemented in RegPy [32, 33], a Python-
based toolbox for implementing and solving (potentially ill-posed) inverse
problems. While the spectral reconstruction can also be achieved using other
methods, this toolbox provides a number of advanced capabilities that we will
exploit later.

For the case that $N_{\tau}=N_{\lambda}$, the reconstruction problem is fully
constrained: for each pixel $x^{\prime}$ we have $N$ complex-valued
measurements and use these to determine $N$ complex-valued amplitudes. This is
also reflected in the measurement results: as shown in Fig. 2b-c, the
$N_{\tau}=11$ measurement allows to accurately extract the spectrally-resolved
holograms from the data. This method already provides a powerful advantages
over the full FTS scan: As a full FTS scan of the diffraction pattern requires
to sample the highest harmonic with 2 points per oscillation and must also
resolve the individual harmonics, the number of samples must be at least twice
the order of the highest harmonic (i.e., 66 measurements for our spectrum with
harmonics up to the 33rd order), the total number of measurements is
dramatically reduced. In practice, the sampling advantage is around an order
of magnitude, as full FTS measurements in the EUV typically use hundreds of
delays. In addition to a shorter measurement time, this implies that the short
scan is much less sensitive to experimental drift, such as in the EUV beam
pointing.

Next, we exploit the intrinsic wavelength-sensitivity of FTH, which leads to
the relative displacement of the individual spectral components that was
discussed earlier. Based on the known width (30 µm) and separation (100 µm) of
the sample and reference, we construct a simple mask to constrain which
spectral components can contribute to each pixel $x^{\prime}$. From this
analysis, we find that at no point more than 6 high harmonics contribute to
the signal (and 99% of the pixels include contributions of $\leq 5$ harmonics
only). We therefore set $N_{\tau}=6$ and use interference patterns spaced in
delay by $1.7/5=0.34$ fs to again cover half cycle of the driving laser pulse.
The mask can be implemented conveniently into the RegPy reconstruction method.
As shown in Fig. 2d, this approach enables spectroscopic imaging at a total
measurement time of only 94 s.

Figure 3: Comparison of image quality for minimally sampled FTSH compared to
FTS with $N_{\tau}=75$ covering a delay range of 3.48 fs. This is the minimum
delay range that allows to separate the individual high harmonics by a Fourier
transform. (a) Experimental signal-to-noise ratio (SNR), calculated by
comparing the average logo intensity to the background. (b) SNR scaled by
$1/\sqrt{T}$, where $T$ is the total measurement time. (c-e) Exemplary
spectrally-resolved hologram for the 19th harmonic at 54 nm. Although the
overall best SNR is achieved from the long integration, the SNR is only
$\approx$2x better than that of the $N_{\tau}=11$ analysis, while the total
exposure time was 7x longer. Assuming that the SNR of the minimally sampled
reconstructions is limited by photon shot noise, minimally sampled FTSH is
expected to yield a 2x higher SNR than the full FTS measurement for equal data
acquisition time.

These results indicate that the sampling requirements in spatially-resolved
Fourier transform spectroscopy, specifically when applied in combination with
Fourier transform holography, is not determined by the well-known Nyquist-
Shannon theorem applied to the HHG source spectrum, but it is rather
determined by the total number of unknowns, which is in turn dependent on the
FTH sample structure, the spectral resolution, and the bandwidth of the EUV
spectrum. For typical HHG spectra and FTH sample structures, this leads to a
significant reduction of the sampling requirements. In practice, the number of
delay steps needed can be controlled by the separation $\Delta$. Increasing
$\Delta$ at fixed sample width leads to a reduced overlap of the individual
spectral components and consequentially a reduced number of steps. However,
the increased $\Delta$ also implies that a larger total field of view must be
considered in FTH. For diffraction-limited imaging at high numerical aperture
(see below), this can also result in a reduction of the final image
resolution. Thus, in general a trade-off between image acquisition speed,
image quality and experimental parameters such as the illumination beamline
must be found.

To quantify the benefit of minimal sampling FTSH compared to high-resolution
FTS, we carry out an analysis of the effective signal-to-noise (SNR) ratio of
the reconstructed images in Fig. 3. Overall, we find that a longer total
integration time leads to an increase in the SNR, which agrees with
expectation, because the effects of shot-noise and shot-to-shot HHG variations
are reduced. However, the improvement from $N_{\tau}=11$ to the full FTS
result is disproportional to the increase in measurement time. Assuming that
the noise is limited by shot-noise, the image quality can be compared if the
SNR is scaled by $1/\sqrt{T}$, where $T$ is the total measurement time (see
Fig. 3c). From this analysis, we find that both minimally-sampled FTSH methods
(11 delays without mask and 6 delays with applied mask, respectively) perform
on average 2x better than the full FTS measurement. We tentatively attribute
this difference to a reduced sensitivity to drifts in the EUV spectrum and
beam pointing for the shorter FTSH measurement.

The presented FTSH method relies on prior knowledge of the high-harmonic
frequencies. Although these can be commonly estimated from the spectrum of the
near-infrared driving pulse, the high-harmonic spectrum can be shifted due to
various effects such as pulse chirp and intensity blue-shift [34, 35].
Therefore, we have performed an analysis of frequency miss-estimation on the
reconstructed electric fields, from which we find that a $\pm 1$% error in the
fundamental frequency can lead to a root-mean-square difference up to 20%.
Such errors can be completely avoided, however, by performing an extended
measurement with at least one extra delay. The number of independent pixels is
normally much larger than the number of spectral components ($10^{5}$ pixels
and $11$ frequency components in our work), and thus this extra image provides
enough data to determine the high-harmonic frequencies precisely through
iterative minimization. With regard to the measurement calibration, also sub-
wavelength accuracy is required in the determination of the delays. This can
for example be achieved using a stable interferometer [30] or through separate
interferometric delay calibration [36]. An advantage of the minimally-sampled
FTSH scheme for typical HHG spectra is that the necessary delay range covers
less than one micrometer, thereby reducing the requirements on the
interferometer.

Figure 4: Diffraction-limited spectroscopic imaging based upon the
$N_{\tau}=6$ dataset with a total measurement time of 1m 34s. (a) Far-field
interference pattern $|\tilde{p}_{\lambda}\tilde{r}_{\lambda}^{*}|$ of the
19th harmonic at 54 nm. (b) The far-field diffraction pattern (Airy pattern)
of the reference structure at 54 nm. (c) Reconstructed far-field diffraction
pattern $|\tilde{p}_{\lambda}|$ of the object, as extracted from (a). (d)
Reconstructed images for the 15th to 25th harmonics, respectively.

##  4 Diffraction-limited spectromicroscopy

A limitation of FTH is the trade-off between resolution (requiring small
reference structures) and signal strength (requiring large references). Hence,
many FTH experiments use subsequent iterative phase retrieval to improve the
spatial resolution and SNR [37, 17, 18]. In the following, we will demonstrate
how iterative image reconstruction methods can be used to achieve diffraction-
limited resolution in FTSH. In a first step, it is useful to separate the
delay-dependent and delay-independent contributions to the diffraction
pattern. The delay-independent part has two contributions:
$\sum_{i}^{N_{\lambda}}|\tilde{p}_{i}|^{2}$ and
$\sum_{i}^{N_{\lambda}}|\tilde{r}_{i}|^{2}$. These terms each represent a
conventional broadband diffraction pattern, for which it is known that image
reconstruction is possible with strong prior knowledge of the spectral
response of the object [38]. Reconstruction of the individual spectral
components, however, is not generally possible. Also, the incoherent sum of
both broadband diffraction patterns will complicate image reconstruction.
Consequentially, we will focus in the following on the spectrally-resolved
interference patterns that are retrieved from the FTSH analysis.

After applying a Fourier transform to the output of the FTSH analysis (such as
shown in Fig. 2d, e), the monochromatic interference-diffraction patterns can
be expressed as
$m_{\lambda}(k)=\tilde{p}_{\lambda}\tilde{r}_{\lambda}^{*}(k)$. As this data
already includes the phase, it is not necessary to implement phase retrieval
algorithms as in monochromatic coherent diffractive imaging. Instead, the
image reconstruction of $p_{\lambda}$ from $m_{\lambda}$ can be posed as a
deconvolution problem, where the point-spread function of the reference
($r_{\lambda}$) needs to be subtracted from the data. This can also be
observed in the extracted monochromatic far-field interference pattern (Fig.
4a), where a clear imprint of the Airy pattern due to the 12 micron reference
can be seen. We emphasize that the clear visibility of the Airy rings is due
to the direct illumination of both the reference and the sample using focused
EUV beams. The shape of the observed Airy pattern matches well to the expected
pattern (Fig. 4b).

To reconstruct the diffraction-limited image from the monochromatic
interferogram $m_{\lambda}(k)$ and the known reference $r_{\lambda}$, we use
the Newton conjugate-gradient method to iteratively minimize
$||\tilde{p}_{\lambda}\tilde{r}_{\lambda}^{*}(k)-m_{\lambda}(k)||$. In order
to prevent overfitting to the noise, we use Morozov’s discrepancy principle,
i.e. we terminate the algorithm as soon as the residuals are on the order of
the noise floor. This is commonly achieved after 11 Newton steps. Using this
procedure, we find accurate image reconstructions for up to 6 high-harmonics,
see Fig. 4. Compared to the initial holograms that were limited by the 12
micron diameter reference, this enables a dramatic increase in the resolution
to a final resolution of approximately 1.5 micron. This resolution corresponds
to the numerical aperture of the observed reference diffraction. In following
experiments, the experimental resolution can further be improved by the
implementation of smaller or structured references, which can be illuminated
efficiently because the EUV beams are focused to the reference and sample
individually. Also, the resolution and SNR can be improved by averaging of
several consecutive measurements (see supplementary information). It is also
worth noting that this image reconstruction procedure does not depend on
knowledge of the object support.

In summary, we have investigated the potential of Fourier-transform
spectroscopic holography (FTSH), an interferometric technique that allows
minute-scale spectromicroscopy at extreme ultraviolet photon energies. In the
FTH geometry, prior knowledge of the illumination wavelengths allows a
dramatic reduction in the sampling requirements. We have identified two
sampling strategies for minimally sampled FTSH: (i) without assuming prior
knowledge on the overlap of monochromatic holograms, the number of
measurements can be reduced to the number of wavelength components in the
illumination, and (ii) by incorporating such prior knowledge, a further
reduction can be achieved. For our HHG spectrum and FTH sample, this reduces
the number of measurements to 6, which is a typical value. Overall, this
approach enables an order of magnitude reduction in the required sampling
compared to full FTS sampling by the Nyquist-Shannon theorem, thereby speeding
up the measurement dramatically. In the spectrally resolved holograms, we find
that a higher SNR is achieved in equivalent time. This is made possible by the
shorter measurement duration that reduces the effect of long-term drift and
instability in the experimental setup. Finally, although the measurement is
based on holography and intrinsically measures the convolution between sample
and reference, FTSH does enable diffraction-limited resolution. We expect that
this method will contribute strongly to the implementation of table-top
extreme ultraviolet spectromicroscopy, for example in time-resolved
experiments where a fast measurement is necessary to facilitate systematic
studies.

Funding This work was funded by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) - 432680300/SFB 1456, projects B01 and C03. F.Z.
and S.W. acknowledge support from the European Research Council (ERC-CoG
project 3D-VIEW, 864016).

The experiments were performed at ARCNL, a public-private partnership between
the University of Amsterdam (UvA), Vrije Universiteit Amsterdam (VU),
Rijksuniversiteit Groningen (RUG), the Dutch Research Council (NWO), and the
semiconductor equipment manufacturer ASML.

Disclosures The authors declare no conflicts of interest.

## References

  * [1] F. Krausz and M. Ivanov, “Attosecond physics,” Rev. Mod. Phys. 81, 163 (2009).
  * [2] L. Loetgering, S. Witte, and J. Rothhardt, “Advances in laboratory-scale ptychography using high harmonic sources [Invited],” Opt Express 30, 4133–4164 (2022).
  * [3] R. Battistelli, D. Metternich, M. Schneider, L.-M. Kern, K. Litzius, J. Fuchs, C. Klose, K. Gerlinger, K. Bagschik, C. M. Günther, D. Engel, C. Ropers, S. Eisebitt, B. Pfau, F. Büttner, and S. Zayko, “Coherent x-ray magnetic imaging with 5 nm resolution,” Optica 11, 234 (2024).
  * [4] M. D. Seaberg, D. E. Adams, E. L. Townsend, D. A. Raymondson, W. F. Schlotter, Y. Liu, C. S. Menoni, L. Rong, C.-C. Chen, J. Miao, H. C. Kapteyn, and M. M. Murnane, “Ultrahigh 22 nm resolution coherent diffractive imaging using a desktop 13 nm high harmonic source,” Optics Express 19, 22470–22479 (2011).
  * [5] M. Schultze, K. Ramasesha, C. Pemmaraju, S. Sato, D. Whitmore, A. Gandman, J. S. Prell, L. J. Borja, D. Prendergast, K. Yabana, D. M. Neumark, and S. R. Leone, “Attosecond band-gap dynamics in silicon,” Science 346, 1348–1352 (2014).
  * [6] F. Siegrist, J. A. Gessner, M. Ossiander, C. Denker, Y.-P. Chang, M. C. Schröder, A. Guggenmos, Y. Cui, J. Walowski, U. Martens, J. K. Dewhurst, U. Kleineberg, M. Münzenberg, S. Sharma, and M. Schultze, “Light-wave dynamic control of magnetism,” Nature 571, 240–244 (2019).
  * [7] R. Geneaux, H. J. B. Marroux, A. Guggenmos, D. M. Neumark, and S. R. Leone, “Transient absorption spectroscopy using high harmonic generation: a review of ultrafast X-ray dynamics in molecules and solids,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 377, 20170463 (2019).
  * [8] E. de Vos, S. Neb, A. Niedermayr, F. Burri, M. Hollm, L. Gallmann, and U. Keller, “Ultrafast Transition from State-Blocking Dynamics to Electron Localization in Transition Metal $\beta$-Tungsten,” Phys Rev Lett 131, 226901 (2023).
  * [9] T. Heinrich, H.-T. Chang, S. Zayko, K. Rossnagel, M. Sivis, and C. Ropers, “Electronic and Structural Fingerprints of Charge-Density-Wave Excitations in Extreme Ultraviolet Transient Absorption Spectroscopy,” Phys. Rev. X 13, 021033 (2023).
  * [10] C. Möller, H. Probst, G. S. M. Jansen, M. Schumacher, M. Brede, J. K. Dewhurst, M. Reutzel, D. Steil, S. Sharma, and S. Mathias, “Verification of ultrafast spin transfer effects in iron-nickel alloys,” Commun Phys 7, 1–6 (2024).
  * [11] H. Probst, C. Möller, M. Schumacher, T. Brede, J. K. Dewhurst, M. Reutzel, D. Steil, S. Sharma, G. S. M. Jansen, and S. Mathias, “Unraveling femtosecond spin and charge dynamics with extreme ultraviolet transverse MOKE spectroscopy,” Physical Review Research 6, 013107 (2024).
  * [12] A. S. Johnson, J. V. Conesa, L. Vidas, D. Perez-Salinas, C. M. Günther, B. Pfau, K. A. Hallman, R. F. Haglund, S. Eisebitt, and S. Wall, “Quantitative hyperspectral coherent diffractive imaging spectroscopy of a solid-state phase transition in vanadium dioxide,” Sci. Adv. 7, eabf1386 (2021).
  * [13] A. S. Johnson, D. Perez-Salinas, K. M. Siddiqui, S. Kim, S. Choi, K. Volckaert, P. E. Majchrzak, S. Ulstrup, N. Agarwal, K. Hallman, R. F. Haglund, C. M. Günther, B. Pfau, S. Eisebitt, D. Backes, F. Maccherozzi, A. Fitzpatrick, S. S. Dhesi, P. Gargiani, M. Valvidares, N. Artrith, F. de Groot, H. Choi, D. Jang, A. Katoch, S. Kwon, S. H. Park, H. Kim, and S. E. Wall, “Ultrafast X-ray imaging of the light-induced phase transition in VO2,” Nat Phys 19, 215–220 (2023).
  * [14] I. McNulty, J. Kirz, C. Jacobsen, E. H. Anderson, M. R. Howells, and D. P. Kern, “High-Resolution Imaging by Fourier Transform X-ray Holography,” Science 256, 1009–1012 (1992).
  * [15] H. He, U. Weierstall, J. C. H. Spence, M. Howells, H. A. Padmore, S. Marchesini, and H. N. Chapman, “Use of extended and prepared reference objects in experimental Fourier transform x-ray holography,” Appl Phys Lett 85, 2454–2456 (2004).
  * [16] S. Marchesini, S. Boutet, A. E. Sakdinawat, M. J. Bogan, S. Bajt, A. Barty, H. N. Chapman, M. Frank, S. P. Hau-Riege, A. Szöke, C. Cui, D. A. Shapiro, M. R. Howells, J. C. H. Spence, J. W. Shaevitz, J. Y. Lee, J. Hajdu, and M. M. Seibert, “Massively parallel X-ray holography,” Nature Photon 2, 560–563 (2008).
  * [17] V. T. Tenner, K. S. E. Eikema, and S. Witte, “Fourier transform holography with extended references using a coherent ultra-broadband light source,” Opt. Express, OE 22, 25397–25409 (2014).
  * [18] J. Geilhufe, B. Pfau, C. M. Günther, M. Schneider, and S. Eisebitt, “Achieving diffraction-limited resolution in soft-X-ray Fourier-transform holography,” Ultramicroscopy 214, 113005 (2020).
  * [19] W. Eschen, S. Wang, C. Liu, R. Klas, M. Steinert, S. Yulin, H. Meißner, M. Bussmann, T. Pertsch, J. Limpert, and J. Rothhardt, “Towards attosecond imaging at the nanoscale using broadband holography-assisted coherent imaging in the extreme ultraviolet,” Communications Physics 4, 1–7 (2021).
  * [20] S. Zayko, O. Kfir, M. Heigl, M. Lohmann, M. Sivis, M. Albrecht, and C. Ropers, “Ultrafast high-harmonic nanoscopy of magnetization dynamics,” Nature Communications 12, 6337 (2021).
  * [21] P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, and F. Pfeiffer, “High-Resolution Scanning X-ray Diffraction Microscopy,” Science 321, 379–382 (2008).
  * [22] B. Zhang, D. F. Gardner, M. H. Seaberg, E. R. Shanblatt, C. L. Porter, R. Karl, C. A. Mancuso, H. C. Kapteyn, M. M. Murnane, and D. E. Adams, “Ptychographic hyperspectral spectromicroscopy with an extreme ultraviolet high harmonic comb,” Opt Express 24, 18745–18754 (2016).
  * [23] L. Loetgering, X. Liu, A. C. C. D. Beurs, M. Du, G. Kuijper, K. S. E. Eikema, and S. Witte, “Tailoring spatial entropy in extreme ultraviolet focused beams for multispectral ptychography,” Optica, OPTICA 8, 130–138 (2021).
  * [24] G. S. M. Jansen, D. Rudolf, L. Freisem, K. S. E. Eikema, and S. Witte, “Spatially resolved Fourier transform spectroscopy in the extreme ultraviolet,” Optica 3, 1122–1125 (2016).
  * [25] G. O. Williams, A. I. Gonzalez, S. Künzel, L. Li, M. Lozano, E. Oliva, B. Iwan, S. Daboussi, W. Boutu, H. Merdji, M. Fajardo, and P. Zeitoun, “Fourier transform holography with high harmonic spectra for attosecond imaging applications,” Opt Lett 40, 3205–3208 (2015).
  * [26] S. G. Kalenkov, G. S. Kalenkov, and A. E. Shtanko, “Hyperspectral holography: an alternative application of the Fourier transform spectrometer,” JOSA B 34, B49–B55 (2017).
  * [27] A. Genco, C. Cruciano, M. Corti, K. E. McGhee, B. Ardini, L. Sortino, L. Hüttenhofer, T. Virgili, D. G. Lidzey, S. A. Maier, A. Bassi, G. Valentini, G. Cerullo, and C. Manzoni, “k-Space Hyperspectral Imaging by a Birefringent Common-Path Interferometer,” ACS Photonics 9, 3563–3572 (2022).
  * [28] Y. Meng, C. Zhang, C. Marceau, A. Y. Naumov, P. B. Corkum, and D. M. Villeneuve, “Octave-spanning hyperspectral coherent diffractive imaging in the extreme ultraviolet range,” Opt Express 23, 28960–28969 (2015).
  * [29] G. S. M. Jansen, A. C. C. d. Beurs, X. Liu, K. S. E. Eikema, and S. Witte, “Diffractive shear interferometry for extreme ultraviolet high-resolution lensless imaging,” Opt Express 26, 12479–12489 (2018).
  * [30] D. Brida, C. Manzoni, and G. Cerullo, “Phase-locked pulses for two-dimensional spectroscopy by a birefringent delay line,” Opt. Lett. 37, 3027–3029 (2012).
  * [31] M.-S. Tsai, A.-Y. Liang, C.-L. Tsai, P.-W. Lai, M.-W. Lin, and M.-C. Chen, “Nonlinear compression toward high-energy single-cycle pulses by cascaded focus and compression,” Science Advances 8, eabo1945 (2022).
  * [32] M. Hanke, “Regularizing properties of a truncated newton-cg algorithm for nonlinear inverse problems,” Numerical Functional Analysis and Optimization 18, 971–993 (1997).
  * [33] S. Maretzke, M. Bartels, M. Krenkel, T. Salditt, and T. Hohage, “Regularized newton methods for x-ray phase contrast and general imaging problems,” Opt. Express 24, 6490–6506 (2016).
  * [34] Z. Chang, A. Rundquist, H. Wang, I. Christov, H. C. Kapteyn, and M. M. Murnane, “Temporal phase control of soft-x-ray harmonic emission,” Physical Review A 58, R30–R33 (1998).
  * [35] H. J. Shin, D. G. Lee, Y. H. Cha, J.-H. Kim, K. H. Hong, and C. H. Nam, “Nonadiabatic blueshift of high-order harmonics from Ar and Ne atoms in an intense femtosecond laser field,” Physical Review A 63, 053407 (2001).
  * [36] Y. Meng, C. Zhang, C. Marceau, A. Y. Naumov, P. B. Corkum, and D. M. Villeneuve, “Interferometric time delay correction for Fourier transform spectroscopy in the extreme ultraviolet,” J Mod Optic 63, 1661–1667 (2016).
  * [37] F. Capotondi, E. Pedersoli, M. Kiskinova, A. Martin, M. Barthelmess, and H. N. Chapman, “A scheme for lensless x-ray microscopy combining coherent diffraction imaging and differential corner holography,” Opt. Express 20, 25152–25160 (2012).
  * [38] J. Huijts, S. Fernandez, D. Gauthier, M. Kholodtsova, A. Maghraoui, K. Medjoubi, A. Somogyi, W. Boutu, and H. Merdji, “Broadband coherent diffractive imaging,” Nat Photonics 14, 618–622 (2020).

2024_EUV_Holo


# Spectral Evolution of Ultraluminous X-ray Pulsar NGC 300 ULX-1

Mason Ng MIT Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Ronald A.
Remillard MIT Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  James F.
Steiner Smithsonian Astrophysical Observatory, 60 Garden St., Cambridge, MA
02138, USA  Deepto Chakrabarty MIT Kavli Institute for Astrophysics and Space
Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Dheeraj R. Pasham MIT Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA

(Accepted October 10, 2022)

###### Abstract

We report on results from a one-year soft X-ray observing campaign of the
ultraluminous X-ray pulsar NGC 300 ULX-1 by the Neutron star Interior
Composition Explorer (NICER) during 2018–2019. Our analysis also made use of
data from Swift/XRT and XMM-Newton in order to model and remove contamination
from the nearby eclipsing X-ray binary NGC 300 X-1. We constructed and fitted
a series of 5-day averaged NICER spectra of NGC 300 ULX-1 in the 0.4–4.0 keV
range to evaluate the long-term spectral evolution of the source, and found
that an absorbed power-law model provided the best fit overall. Over the
course of our observations, the source flux (0.4–4.0 keV; absorbed) dimmed
from $2\times 10^{-12}$ to below $10^{-13}{\rm\,erg\,s^{-1}\,cm^{-2}}$ and the
spectrum softened, with the photon index going from $\Gamma\approx 1.6$ to
$\Gamma\approx 2.6$. We interpret the spectral softening as reprocessed
emission from the accretion disk edge coming into view while the pulsar was
obscured by the possibly precessing disk. Some spectral fits were
significantly improved by the inclusion of a disk blackbody component, and we
surmise that this could be due to the pulsar emerging in between obscuration
episodes by partial covering absorbers. We posit that we observed a low-flux
state of the system (due to line-of-sight absorption) punctuated by the
occasional appearance of the pulsar, indicating short-term source variability
nested in longer-term accretion disk precession timescales.

stars: neutron – stars: oscillations (pulsations) – X-rays: binaries – X-rays:
individual (NGC 300 ULX-1)

††journal: ApJ††facilities: NICER, XMM, Swift††software: Astropy (Astropy
Collaboration et al., 2013, 2018), NumPy and SciPy (Virtanen et al., 2020),
Matplotlib (Hunter, 2007), IPython (Perez & Granger, 2007), tqdm (da Costa-
Luis et al., 2022), HEASoft 6.29333http://heasarc.gsfc.nasa.gov/ftools (NASA
High Energy Astrophysics Science Archive Research Center (2014), HEASARC)

##  1 Introduction

Ultraluminous X-ray sources (ULXs) are extremely bright, off-nuclear objects
that have luminosities (assuming isotropic emission) of at least
$L_{X}>10^{39}{\rm\,erg\,s^{-1}}$ (see Kaaret et al., 2017, for a review),
exceeding the Eddington limit of a 1.4 $M_{\odot}$ neutron star ($L_{X}\approx
2\times 10^{38}{\rm\,erg\,s^{-1}}$), and are believed to be powered by super-
Eddington accretion onto compact objects (Gladstone et al., 2009; Sutton et
al., 2013; Bachetti et al., 2014; Fabrika et al., 2021; Qiu & Feng, 2021).
Hundreds of ULXs and over a thousand ULX candidates have been discovered thus
far (Pasham et al., 2014, 2015; Earnshaw et al., 2019; Kovlakas et al., 2020;
Walton et al., 2022), and they were originally thought to be entirely due to
sub-Eddington accreting intermediate-mass black hole binaries based on the
observed luminosities (Colbert & Mushotzky, 1999; Körding et al., 2002;
Strohmayer & Mushotzky, 2009).

However, the discovery of coherent pulsations in M82 ULX-2 demonstrated
unequivocal evidence that the compact object in at least some ULXs is a
neutron star (Bachetti et al., 2014). In fact, it has been proposed that
neutron star accretors in ULXs comprise a significant fraction of the ULX
population (Mushtukov et al., 2015a; King & Lasota, 2016; Koliopanos et al.,
2017; Walton et al., 2018a). Since M82 ULX-2, we have discovered ULX pulsars
that show persistent pulsations (“persistent” ULX pulsars; Fürst et al., 2016;
Israel et al., 2017; Carpano et al., 2018; Sathyaprakash et al., 2019;
Rodríguez Castillo et al., 2020) or pulsations that manifest in outbursts
(“transient” ULX pulsars; Tsygankov et al., 2017; Wilson-Hodge et al., 2018;
Vasilopoulos et al., 2020). The discovery of ULX pulsars have forced a
revision in our understanding of the emission mechanisms in these extreme
regimes. Proposed mechanisms include the presence of magnetar-like magnetic
fields of $B\sim 10^{14}$ G (Mushtukov et al., 2015b), or the beaming of a
weak-field ($B\sim 10^{11}$ G) pulsar accreting at super-Eddington rates
(Kluzniak & Lasota, 2015; King & Lasota, 2016). In this work, we considered in
particular NGC 300 ULX-1, a high-mass X-ray binary hosting a pulsar (Binder et
al., 2011b; Heida et al., 2019).

NGC 300 ULX-1 (hereafter called ULX-1) hosts a ULX pulsar originally
associated with supernova impostor SN 2010da (Binder et al., 2011b, 2020).
Pulsations with an average value of $31.6$ s were discovered during
simultaneous XMM-Newton and NuSTAR observations in December 2016 (Carpano et
al., 2018). Monitoring with Chandra, NICER, XMM-Newton, and Swift showed that
the spin period evolved from 126 s to 18 s over 4 years (Vasilopoulos et al.,
2018). ULX-1 has also been shown to exhibit timing “anti-glitches” (i.e.,
sudden spin-down) during its spin-up due to accretion (Ray et al., 2019).
Vasilopoulos et al. (2019) has also reported on the long-term spin evolution
of ULX-1, and showed evidence that the pulsar had been experiencing a constant
mass accretion rate (inferred from roughly constant luminosity) even though
the X-ray flux decreased by a factor of $\sim 50$. This behavior is similar to
that seen in the long-term X-ray flux variations of ULX pulsar NGC 7793 P13
(Fürst et al., 2021). For both ULX-1 and NGC 7793 P13, the authors conjectured
that this behavior was due to obscuration of the pulsar possibly due to an
outflow (Kosec et al., 2018) or Lense-Thirring precession of the inner
accretion disk (Middleton et al., 2018, 2019; Khan et al., 2022). This
supported the work by Carpano et al. (2018), who found a large difference in
absorption column density between X-ray observations in 2010 (XMM-Newton) and
2016 (XMM-Newton and NuSTAR), while keeping other spectral parameters tied.
Optical observations by the Very Large Telescope/X-shooter have revealed the
companion object to be a red supergiant (Heida et al., 2019).

In this work, we used data from NICER, Swift, and XMM-Newton to characterize
the X-ray spectral evolution of ULX-1. However, a complication in studying
ULX-1 with non-imaging detectors is the presence of NGC 300 X-1 (hereafter
called X-1), a black hole-Wolf Rayet 32.8 hr eclipsing binary, which is 1.4′
away from ULX-1, and is well within the same NICER field-of-view (see Figure
1). Thus to account for the emission from the nearby X-1 in order to isolate
the ULX-1 emission, we used Swift/XRT measurements to create an orbital timing
model with which we applied to dense high signal-to-noise XMM-Newton/EPIC data
to characterize the X-1 spectrum. We generated mock spectral data, based on
the X-1 spectral parameters, folded through the NICER spectral response. To
isolate the spectrum of ULX-1, we subtracted both the modeled component from
X-1 and the background spectrum (mostly due to particle interactions with the
NICER detectors), determined with the “3C50” model (Remillard et al., 2022),
from each raw spectrum selected in this investigation.

In this paper, we presented the spectral evolution of NGC 300 ULX-1 across
twelve months of observations, completing the analysis of all available NICER
data, by constructing spectra from 5-day averaged data. We took advantage of
NICER’s sensitive coverage in soft X-rays to characterize the source as it
faded by a factor of $\sim 50$ across the observation span. In § 2, we outline
the observations taken for ULX-1 and expound on the data analysis. We present
the results in § 3, which will be discussed fully in § 4.

##  2 Observations and Data Analysis

###  2.1 Swift/XRT Observations

We made use of Swift/XRT 0.3–10 keV observations in the Photon Counting (PC)
readout mode taken over 2016 April 14 to 2019 May 11 with ObsIDs
498340[02-89], 88651001, and 88810002. As detailed below, we were looking for
a known periodicity on the order of 32 hr, so the 2.5 s time resolution was
sufficient, and we required imaging and spectroscopic capabilities (Burrows et
al., 2005). These data were processed by the standard processing pipeline with
the good time interval (GTI) expression criteria $\texttt{ELV}\geq 28$,
$\texttt{BR\\_EARTH}\geq 120$, $\texttt{SUN\\_ANGLE}\geq 45$,
$\texttt{ANG\\_DIST}\leq 0.15$, and $\texttt{MOON\\_ANGLE}\geq 14$.

Figure 1: Integrated image of the NGC 300 field from Swift/XRT observations.
The image was generated on SAOImageDS9 and smoothed with a Gaussian kernel
with radius 3. Events of X-1 (yellow; left) and ULX-1 (blue; right) were taken
from extraction regions 30″ and 35″ in radii, respectively. We used a nearby
120″ region for the background of both sources (magenta). The NICER field-of-
view (30 arcmin2; white) is also shown for reference, centered on NGC 300
ULX-1. The best known source position for NGC 300 X-1 is $\alpha_{\rm
J2000}=00^{\rm h}55^{\rm m}10\fs 00$ and $\delta_{\rm J2000}=-37\arcdeg
42\arcmin 12\farcs 2$ (Binder et al., 2011a), and for NGC 300 ULX-1 it is
$\alpha_{\rm J2000}=00^{\rm h}55^{\rm m}04\fs 85$ and $\delta_{\rm
J2000}=-37\arcdeg 41\arcmin 43\farcs 5$ (Binder et al., 2011b).

As shown in Figure 1, we used 30″ and 120″ extraction regions for X-1 and the
background, respectively, across all the observations; light curves were then
extracted with XSELECT, and rate corrections were applied to the light curves
with xrtlccorr to account for bad pixels, or point spread function
corrections. We also applied barycenter corrections with barycorr with the JPL
DE405 solar system ephemeris (Standish, 1998). For the Swift/XRT events, we
used swxpc0to12s6_20130101v014.rmf and swxpc0to12s6_20010101v013.arf for the
redistribution matrix (rmf) and ancillary response (arf) files, respectively.

###  2.2 XMM-Newton Observations

We used XMM-Newton/EPIC 0.3–10.0 keV observations of the NGC 300 field which
spanned 4 days, from 2016 December 17 to 2016 December 20, with ObsIDs
0791010101 and 0791010301. The XMM-Newton Science Analysis System (SAS)
v20.0.0 was used for the data reduction, in conjunction with the most up-to-
date calibration files as of 2022 June. The EPIC data were reprocessed with
epproc and emproc with standard filtering criteria to obtain new event lists.
The three EPIC cameras were operated in the imaging data mode (with the medium
filter), and we imposed standard filtering criteria where
$\textsc{pattern}\leq 12$ for MOS cameras, and $\textsc{pattern}\leq 4$ and
$\textsc{flag}=0$ for the PN camera. We also filtered the event files for
background flares, where we visually inspected the 10.0–12.0 keV light curves
to determine the threshold, and applied rate cuts of 0.35, 0.15, and 0.20 c/s
for PN, MOS1, and MOS2 cameras, respectively.

The source events for X-1 were extracted using circular regions with radius
25″ for all instruments, and the background was extracted from a nearby
source-free region with radius 50″. We also performed barycenter corrections
to the cleaned event files with the SAS task barycen, again with the JPL DE405
solar system ephemeris (Standish, 1998). The response files for PN and MOS
event files were generated with rmfgen and arfgen.

###  2.3 NICER Observations

In this subsection, we describe data processing steps that were undertaken
with special care to optimize spectral analyses for a source as faint as ULX-1
(0.2–2.0 count s-1 for 50 FPMs in the 0.3–4.0 keV range). We also note that
NICER is an external payload attached onto the International Space Station
(ISS) with 56 co-aligned X-ray concentrators and focal plane modules (FPMs)
with silicon drift detectors, of which 52 are operable. The data were
processed with version 8 of the NICER Data Analysis Software (NICERDAS), in
conjunction with version 6.29 of HEASOFT. We first describe steps to define
GTIs, which established the time windows for spectral extractions and
background predictions, and then we applied filters to background-subtracted
spectra to screen out GTIs that showed signs of inaccuracy in the background
prediction.

We investigated NICER observations of ULX-1 performed from 2018 May 1 through
2019 May 2 (ObsIDs 1034200102 through 2034200205). These data were first
brought to a uniform gain calibration by running the HEASARC tool nicerl2 on
each observation. The calibration file was “nixtiflightpi20170601v006.fits”,
and we note that the last 3 calibrations in this sequence (i.e., v005, v006,
and v007) produced only subtle differences, with essentially no changes to the
energy conversions in the range suitable for spectral fitting (i.e., 0.3–12.0
keV). The nicerl2 task re-runs the data pipeline to determine the time and
equivalent photon energy of each recorded event, while also regenerating the
“unfiltered” and “cleaned” event lists for the GTIs found in the reprocessing
steps. In the present investigation, we adopted all of the default settings
for geometric filters: pointing offset from the target ($<0\fdg 015)$, absence
in the South Atlantic Anomaly, and minimum angles to the Earth limb and to the
solar-illuminated Earth limb ($<15\arcdeg$, and $<30\arcdeg$, respectively).
However, the three detector-related filters - maxima for overshoot event
rates, undershoot rates, and the relationship between the overshoot rate and
the magnetic shielding index - were effectively disabled by specifying
impossibly high values for each (i.e., 15000) on the nicerl2 command line. The
reason for this was to increase the initial exposure time given consideration,
with the intent to filter for data quality at a later time, using background-
subtracted spectra.

After reprocessing, a list of GTI times was tabulated by running nimaketime,
independently, on each observation, using the same combination of filter
parameters used previously. However, when defining the final GTIs for spectral
extraction purposes, two additional steps were added to facilitate accuracy in
the background predictions. We used the “3C50” model (Remillard et al., 2022),
an empirical model based on spectral libraries formulated from NICER
observations of the designated background fields (RXTE_BKGD fields). The
adjustments to GTI times were intended to limit the dynamic range of 3C50
parameter values, within a GTI, so that non-linearities in the background
behavior might be contained when the background model is applied. First, it
was found that transitions in ISS day vs. night can reduce the effectiveness
of the “nz” parameter (raw count rate in 0.0–0.25 keV) for predicting the soft
X-ray excess tied to the optical light leak, when such transitions occur
within a GTI. In particular, solar reflections off ISS structures in advance
of ISS sunrise can cause an early onset of noise increase, while detector
noise bleeding after ISS sunset can broaden the noise curve to affect times
after a day/night transition. To avoid these effects, we ran nimaketime twice,
adding constraints, SUNSHINE=1 and SUNSHINE=0, to each trial, respectively. An
initial table of GTI times (MET START and MET STOP) was obtained by converting
each GTI file from FITS to ASCII format, while adding a third column to
distinguish the SUNSHINE value, and then combining all of the results into one
table. GTI sequences can sometimes show brief gaps caused by isolated packet
loss or by noise in nimaketime selection parameters, and such gaps can be
integrated over without significant consequences. We chose to ignore any gap
of 2 s or less, when there was no SUNSHINE transition, and we masked out
exposure times $\pm 30$ s from time-adjacent GTIs that are associated with a
SUNSHINE transition.

The other step to improve background predictability is to limit the time range
given to each 3C50 model prediction. This is motivated by the fact that NICER
background count rate and spectral shape routinely vary in complex ways over
the course of the ISS orbit, and also from orbit to orbit. Thus, a single
background prediction should not be made for a large interval in ISS orbital
phase. On the other hand, the 3C50 model parameters for the stage 1 library
(i.e., $ibg$ and $hrej$ \- in-focus 15.0–18.0 keV events and rejected particle
events near the detector edge, respectively) suffer from low count rates,
creating an opposing motivation to integrate as long as possible, to limit the
effects of Poisson noise. Noting the four passages between the Earth equator
and the highest polar latitudes ($52\arcdeg$) in the ISS orbit (92 min), we
chose a target GTI interval of 300 s. Intervals of duration ($dt$) longer than
450 s were subdivided into $N$ GTIs, with $N=int(dt/300+0.5)$. The final GTI
table was further limited to durations of at least 50 s. The conversion of the
initial GTI table to the final one was done with a C program written for this
purpose. The final GTI table was indexed, and the index number was included in
the filenames of all downstream data products associated with a given GTI.

For ULX-1, we applied these steps to the selected archive, and this yielded
1057 GTIs and a total exposure of 300.8 ks. This was the starting point for
spectral extractions and the application of the 3C50 background model, each
conducted per GTI. Extractions of spectra and 3C50 model parameters values
were made on the basis of 50 selectable FPMs, while ignoring FPMs 14 and 34,
as explained in Remillard et al. (2022). For occasions when one or more
selectable FPMs were not operating, for a given GTI, then the photon spectra
and measurement parameters were linearly rescaled to 50 FPMs, prior to further
consideration.

Investigations of NICER spectra should be made with screening efforts to
ignore GTIs that have an inaccurate background prediction. Current screening
methods use background-subtracted spectra (noted as parameters with subscript
“net”), per GTI. Filtering is performed in diagnostic energy bands that have
very low effective area from the NICER optics, leading to the expectation that
the net count rate should be near zero. Filtering levels are then advised with
criteria depending on source brightness (Remillard et al., 2022).

In the current investigation, spectral extractions and background modeling
were performed for each of the 1057 GTIs. Filtering was then applied at level
3 (Remillard et al., 2022), selecting $|hbg_{net}|<0.05$ and $|S0_{net}|<2.0$,
where $hbg_{net}$ is the background-subtracted count rate at 13.0–15.0 keV,
and $S0_{net}$ is the analogous rate at 0.2–0.3 keV. This yielded 908 GTIs for
ULX-1 with a total exposure time of 263.0 ks, which is 87% of the total. These
were the parent set of GTIs used for the 5-day binning analyses, described in
Section 3.2.

In this work, we concentrated on the spectral analysis, which was carried out
with XSPEC 12.12.1 (Arnaud, 1996). We restricted the energy range to 0.4–4.0
keV due to low source counts in the higher energy band ($>4.0$ keV) and in
order to adopt a uniform energy range across all spectra. The spectra were
each grouped using the optimal binning scheme and rebinned to have at least 1
count per bin (Kaastra & Bleeker, 2016). The fit statistic employed was the
Cash statistic, suitable for the low count rates we are dealing with (Cash,
1979; Kaastra, 2017). We also applied a 1% systematic error to the spectra as
recommended by the NICER
team111https://heasarc.gsfc.nasa.gov/docs/nicer/analysis_threads/cal-
recommend/.

##  3 Results

We show the flux history of X-1, and ULX-1 from the Swift/XRT 0.2–10.0 keV
events, as well as illustrate the timeline of XMM-Newton and NICER
observations, in Figure 2. X-1 exhibited a constant flux in the long-term,
whereas ULX-1 clearly showed a declining flux over the course of the NICER
observations.

Figure 2: 0.2–10.0 keV Swift/XRT background-subtracted and exposure-corrected
light curve of NGC 300 X-1 (red) and NGC 300 ULX-1 (blue), binned with 1 hr
time bins. The grey shaded regions indicate the time span of the XMM-Newton
data (left) and NICER data (right) used in this work.

###  3.1 NGC 300 X-1

We attempted to look for variability in the hardness ratio (1.0–2.0
keV/0.4–1.0 keV) that would suggest spectral changes. We constructed 5-day
bins with the Swift/XRT data set and did not find any evidence of systematic
variability of X-1 within the uncertainties (not plotted). We similarly did
not find any systematic trend with the hardness-intensity diagram (with
intensity defined over 0.3–10.0 keV).

In order to derive an orbital timing ephemeris to construct orbital phase-
resolved spectra for X-1, we calculated a Lomb-Scargle periodogram using the
Swift/XRT events since we were dealing with an irregularly sampled light curve
(Lomb, 1976; Scargle, 1982; Horne & Baliunas, 1986). First, we noted that
Carpano et al. (2019) combined 86 Swift/XRT PC mode observations over 2006
September 5 to 2018 June 19, 7 XMM-Newton observations over 2000 December to
2016 December, and 5 Chandra observations over 2006 June and 2014 November, to
derive an orbital period for X-1 of $32.7932\pm 0.0029$ h ($1\sigma$). In
constructing our ephemeris, we then focused our search around $32.79$ h. To
determine the best-fit orbital frequency for our observation, we constructed a
grid of trial frequency values and calculated

$\chi^{2}=\sum\limits_{i=1}^{N}\frac{(x_{i}-\bar{x})^{2}}{\sigma_{i}^{2}},$
(1)

where $N$ corresponds to the number of bins in the folded profile, $x_{i}$ is
the count rate in the i-th bin, $\sigma$ is the associated uncertainty, and
$\bar{x}$ is the mean count rate across the profile. The uncertainty was
calculated from the flux randomization and random subset selection methods
(Peterson et al., 1998; Welsh, 1999; Peterson et al., 2004), which accounted
for the sampling and measurement uncertainties of the light curves, where we
ended up with $\sigma=0.002$ $\mu$Hz. Thus we found $P_{\rm orb}=32.788\pm
0.007$ h, which is formally consistent with the value derived by Carpano et
al. (2019). The final folded orbital profile is shown in Figure 3. In order to
derive $T_{0}$, the time of eclipse minimum, we fitted the ramp-and-step model
for X-1. The model describes the folded orbital light curve with a constant
count rate for the off-eclipse and on-eclipse components, and a linear count-
rate evolution during ingress and egress (Wachter et al., 2000; Iaria et al.,
2018). We thus obtained $T_{0}={\rm MJD\,}57493.476\pm 0.017$ (TDB) by adding
the phase offset due to the eclipse minimum (multiplied by the best-fit
orbital period) to the arrival time of the first photon.

Figure 3: Folded exposure-corrected 0.3–10.0 keV orbital profile with $P_{\rm
orb}=32.788\pm 0.007$ h and $T_{0}=$ MJD 57493.476 $\pm 0.017$ (TDB). The
orbital profile has been shifted so that phase 0 corresponds to the eclipse
minimum. The orbital profile for the black hole-Wolf Rayet binary NGC 300 X-1
was constructed using events from the Swift/XRT observations that span a time
interval that included XMM-Newton and NICER data. Two orbital cycles were
plotted for clarity.

With this derived orbital timing ephemeris for X-1, we defined the “on-
eclipse” and “off-eclipse” segments of X-1, where orbital phases 0.15–0.70
refer to the “off-eclipse” and orbital phases 0.70–1.15 (modulo 1) refer to
the “on-eclipse”, where the black hole is obscured by the dense winds of the
companion Wolf-Rayet star. There were too few events to do detailed modeling
of the ingress and egress states. We assumed that the spectral behavior did
not change significantly throughout the Swift/XRT observation period, which
can be inferred from roughly constant spectral hardness in the Swift/XRT
spectral data (see first paragraph).

Events for each of these segments were extracted using XSELECT and finally
binned into a spectrum, for each of the XMM-PN, XMM-MOS1, and XMM-MOS2
instruments, imposing a grouping with a minimum of 1 count per bin and with
the ‘optimal binning’ method (Kaastra & Bleeker, 2016). The “on-eclipse” and
“off-eclipse” spectra for each of the 3 instruments were simultaneously fit
with the C-statistic, i.e., a total of 6 spectra. We found that the model
comprising of an absorbed disk blackbody with a Comptonization component and
two multiplicative constants best described the spectra, where we employed
SIMPL to model the Comptonization component of the plasma surrounding the
black hole binary in a simple, self-consistent way (Steiner et al., 2009). The
two multiplicative constants accounted for two different physical effects -
the first parameterizes the ratio of the X-1 flux during “on-eclipse” (set to
1) and “off-eclipse”; the second parameterized the calibration differences
between the three XMM-Newton instruments, with the constant for XMM-PN fixed
at 1. Single-component models (with a power law or blackbody) and simple two-
component models (absorbed power law and blackbody) inadequately described the
X-1 spectra.

We found a reasonable fit with a reduced C-stat of 1.26 (550 d.o.f.), with a
photon index of $\Gamma=2.32_{-0.04}^{+0.04}$, a scattered fraction of
$f=0.35$, inner disk temperature of $T_{\rm in}=0.112_{-0.013}^{+0.017}$ keV,
on-eclipse and off-eclipse hydrogen column densities of
$n_{H}=0.082_{-0.017}^{+0.027}\times 10^{22}{\rm\,cm^{-2}}$ and
$n_{H}=0.10_{-0.02}^{+0.03}\times 10^{22}{\rm\,cm^{-2}}$, respectively, and a
$E=0.949_{-0.013}^{+0.012}$ keV emission line characterized by a Gaussian with
width $\sigma=0.089_{-0.014}^{+0.016}$ keV. The results of the fit are
summarized in Table 1.

In order to characterize the phase-resolved background contribution from X-1
to use with the NICER data, we used the best-fitting spectral model (see Table
1) and generated NICER response-folded spectra (using fakeit in XSPEC) that
was applied to every spectrum that was extracted for each GTI (defined in §
2.3). We made use of the centroid time of the GTIs to determine the orbital
phase and hence whether the GTI was in the on-eclipse or off-eclipse phase,
which was valid given the short length of the GTIs (average of 290 s) relative
to the much longer orbital period for X-1 ($\sim 32.8$ h). When modeling the
X-1 spectra as an additional component of ”background” for ULX analyses, we
assumed no effective area losses for the X-1 spectrum, which was reasonable
given that X-1 and ULX-1 are 1.4′ apart, and NICER’s vignetting profile
suggested little attenuation in the
response222https://heasarc.gsfc.nasa.gov/docs/nicer/data_analysis/
workshops/NICER-CalStatus-Markwardt-2021.pdf (Okajima et al., 2016).

Table 1: Best-fit spectral parameters from simultaneously fitting 0.3–10.0 keV XMM-Newton spectra for the “off-eclipse” and “on-eclipse” components of NGC 300 X-1. The spectra are described with a phenomenological model with the model tbnew(simpl(diskbb+gaussian))*constant in XSPEC parlance.  Parameter | Value | Units
---|---|---
Tied |  |
Photon Index | $2.32_{-0.04}^{+0.04}$ |  $\cdots$
Scattered Fraction | $0.35_{-0.10}^{+0.10}$ |  $\cdots$
$T_{\rm in}$ | $0.112_{-0.013}^{+0.017}$ | keV
norm (diskbb) | $104_{-61}^{+151}$ |  $\cdots$
Constants |  |
Flux | $0.50_{-0.02}^{+0.03}$ |  $\cdots$
MOS1 | $1.08_{-0.03}^{+0.03}$ |  $\cdots$
MOS2 | $1.14_{-0.03}^{+0.03}$ |  $\cdots$
On-eclipse |  |
$n_{H}$ | $0.082_{-0.017}^{+0.027}$ | $10^{22}{\rm\,cm^{-2}}$
Off-eclipse |  |
$n_{H}$ | $0.10_{-0.02}^{+0.03}$ | $10^{22}{\rm\,cm^{-2}}$
LineE (keV) | $0.949_{-0.013}^{+0.012}$ | keV
Sigma (keV) | $0.089_{-0.014}^{+0.016}$ | keV
norm (Gauss) | $1.40_{-0.18}^{+0.21}\times 10^{-5}$ | photons s-1 cm2
C-stat/d.o.f. | 693.92/550 |  $\cdots$

###  3.2 NGC 300 ULX-1

The nominal 3C50 background predictions for ULX-1 and the NGC 300 X-1 emission
model were then combined to create a ‘total background’ measurement for the
NGC 300 field to isolate the emission from ULX-1. The spectra from the
individual GTIs were then grouped into 5-day bins in order to ascertain the
long-term, average behavior of the source and to increase the signal-to-noise.
The 0.4–12.0 keV light curve in Figure 4 shows a gradual decrease in flux by a
factor of $\sim 50$ over the one-year span of the data. In Figure 5, we show
the hardness-intensity diagram, where we represent spectral hardness with the
soft color, defined as the ratio of 1.0–2.0 keV flux to that of 0.4–1.0 keV
flux. The hardness-intensity diagram shows a clear dimming of the source as
the source softened (i.e., decreasing soft color).

Figure 4: 5-day averaged 0.4–12.0 keV light curve of ULX-1 over the span of 12
months of NICER monitoring (2018 May 1 to 2019 May 2), where the particle and
X-1 background contributions have been subtracted. The flux from ULX-1 broadly
varied by a factor of 25 throughout the span of the NICER observations. Figure
5: Hardness-intensity diagram of ULX-1 with NICER data where the intensity
(0.4–12.0 keV) is plotted as a function of the soft color ((1.0–2.0
keV)/(0.4–1.0 keV)). The intensity is derived using the 5-day averaged
spectra. The count rates in each energy band have been corrected for the
contributions of the particle and X-1 background (see text for details). There
is a clear correlation between the broadband intensity and soft color, where
the dimming of the source is accompanied by a steeper spectrum in soft X-rays
(i.e., 0.4–2.0 keV).

To understand the long-term evolution of the source, we analyzed the 5-day
binned spectra. The source spectra from the individual GTIs in each 5-day bin
were simply summed together in counts space with the mathpha tool, whereas the
background spectra in the same bin were calculated via a exposure time-
weighted average of the individual spectra (i.e., X-1 plus particle
background) in rate space also with mathpha. We ended up with 39 spectra, and
we performed spectral fits with an absorbed power law and found that it was a
good fit overall. We fixed the value of the hydrogen column density at
$n_{H}=0.11\times 10^{22}{\rm\,cm^{-2}}$ (Carpano et al., 2018). The results
of the spectral evolution can be seen in Figure 7 with tbabs(powerlaw), where
panel (a) shows the evolution of the power law photon index ($\Gamma$), panel
(b) is the power law normalization (in
$10^{-4}{\rm\,photons\,keV^{-1}\,cm^{-2}}$), panel (c) shows the reduced
C-statistic (C-stat/number of degrees of freedom), and panel (d) plots the
0.4–4.0 keV absorbed flux in ${\rm erg\,s^{-1}\,cm^{-2}}$. We also tried
fitting the spectra with absorbed blackbody models (with bbodyrad and diskbb)
and found that they were poor fits with reduced C-statistic values above 2.

The spectral evolution in Fig. 7 matches the hardness-intensity diagram from
Fig. 5, where in Fig. 7a, the photon index ($\Gamma$) increased from about
$\Gamma\sim 1.6$, and reached a maximum value of around $\Gamma\sim 4$. The
flux meanwhile dropped by a factor of over 50, representing dramatic dimming
over the course of the observations. The change in the spectral shape between
the bright and faint states of ULX-1 can be seen in Fig. 6, where the bright
state is shown by the red points (MJD 58249) and the faint state is
represented in blue (MJD 58489).

Figure 6: 0.4–4.0 keV spectra corresponding to the brighter state of ULX-1
(MJD 58249) and the fainter state (MJD 58489). The absorption column density
was fixed at $n_{H}=0.11\times 10^{22}{\rm\,cm^{-2}}$, and the best-fit power
law indices were $\Gamma=1.69_{-0.09}^{+0.09}$ and $\Gamma=4.0_{-0.5}^{+0.7}$
for the brighter and fainter states, respectively.

While an absorbed power law was generally a good fit to the spectra throughout
the observation span (see Fig. 7c), some spectra were significantly improved
by the introduction of an additional thermal disk blackbody (diskbb in XSPEC
parlance) component. We evaluated the improvement of one model relative to the
other by considering the difference in the Akaike information criterion (AIC;
Akaike, 1974), i.e., $\Delta{\rm AIC}_{12}={\rm AIC}_{1}-{\rm AIC}_{2}$, where

${\rm AIC}=-2\text{ln}\mathcal{L}+2d,$ (2)

where $\mathcal{L}$ is the maximum likelihood of the given model
($\mathcal{L}\propto e^{-C/2}$; $C$ for C-stat), and $d$ is the number of
model parameters. We can also reject the competing model at a confidence level
by evaluating

$P_{1}/P_{2}\approx e^{-\Delta{\rm AIC}_{12}/2};$ (3)

for example, $\Delta{\rm AIC}_{12}=10$ corresponds to $P_{1}/P_{2}\approx
0.0067$ and is the commonly adopted threshold for decisively distinguishing
two models (Liddle, 2007; Tan & Biswas, 2012; Arcodia et al., 2018). The
corresponding MJDs, C-statistic values from the fits, as well as the
$\Delta{\rm AIC}_{12}$ values, are reported below in Table 2. We only reported
data points where $\Delta{\rm AIC}_{12}>10$.

Table 2: List of 5-day binned spectra whose fits were significantly improved with an additional thermal disk blackbody model. We only reported the C-statistic values from fits where the inclusion of a disk blackbody component resulted in a $\Delta{\rm AIC}_{12}>10$.  MJD | Temperature | $F_{\rm diskBB}/F_{\rm total}$ | C-stat/d.o.f. | C-stat/d.o.f. | $\Delta{\rm AIC}_{12}$
---|---|---|---|---|---
| keV | (0.4–4.0 keV) | (PL only) | (PL + diskBB) |
58239 | $0.30_{-0.04}^{+0.04}$ | 0.15 | 75.1/56 | 53.0/54 | 18.2
58244 | $0.34_{-0.03}^{+0.03}$ | 0.20 | 87.4/60 | 45.8/58 | 35.6
58254 | $0.27_{-0.04}^{+0.04}$ | 0.16 | 69.1/55 | 46.0/53 | 17.9
58454 | $0.19_{-0.05}^{+0.06}$ | 0.27 | 53.9/49 | 29.8/47 | 19.6
58484 | $0.19_{-0.05}^{+0.05}$ | 0.20 | 64.0/50 | 27.4/48 | 32.1
58504 | $0.21_{-0.04}^{+0.05}$ | 0.30 | 62.0/53 | 34.4/51 | 23.5
Figure 7: Spectral evolution of NGC 300 ULX-1 with an absorbed power law model
as a function of time (MJD), using 1 year of NICER observations with the
contribution from X-1 accounted for (details in text). a): Power law index
$\Gamma$; b) power law normalization, in units of $10^{-4}$ photons/keV/cm2 at
1 keV; c) ratio of the C-statistic to the number of degrees of freedom
(d.o.f.); d) 0.4–4.0 keV absorbed flux in ${\rm erg/s/cm^{2}}$.

##  4 Discussion

We have characterized the 1-year spectral evolution of the dimming ULX pulsar
NGC 300 ULX-1 with NICER. Due to the contaminating emission by the nearby
(1.4′) eclipsing X-ray binary, NGC 300 X-1, we used Swift/XRT observations and
high-quality XMM-Newton data to generate orbital phase-dependent NGC 300 X-1
(background) spectra to subtract from the raw NICER spectrum of the NGC 300
field. We showed that over time, the softening of ULX-1 accompanied the
dimming down to almost undetectable levels in Figs. 4 and 7, where the flux
dropped by a factor of over 50, and the photon index varied from
$\Gamma\approx 1.6$ to $\Gamma\approx 2.6$ with a significant scatter (see
Fig. 7).

Ray et al. (2019) investigated the spectral evolution of the pulsed component
over an observation epoch that overlapped with the first 40 days of our
observation span (MJD 58240 to 58280) and found similar results to this work.
However, in their spectral fits, they varied $n_{H}$ along with $\Gamma$,
which introduced a degeneracy in the spectral fitting that made it non-trivial
to evaluate the contribution of the power law. Their work also supported our
approach of increasing the bin size (to five days) as the source was evolving
sufficiently slowly over five days.

The picture that has been emerging for the evolution of ULX-1 is of a pulsar
that is being obscured as a result of outflows and/or Lense-Thirring
precession causing the plane of the accretion disk to precess along the line-
of-sight (Carpano et al., 2018; Middleton et al., 2018; Kosec et al., 2018;
Vasilopoulos et al., 2019; Binder et al., 2020). The Lense-Thirring precession
is thought to arise from super-Eddington accretion ($>20\dot{M}_{\rm Edd}$)
building up a large scale-height wind cone outside of the magnetosphere
(Middleton et al., 2018; Vasilopoulos et al., 2019; King & Lasota, 2020).
Other precession mechanisms are also possible, although they are not likely to
dominate the Lense-Thirring precession (Middleton et al., 2018; Vasilopoulos
et al., 2018; Carpano et al., 2018; Walton et al., 2018b; Pan et al., 2022).
As a result of the obscuration, however, we would also expect to see an
increased spectral hardening due to the soft photons being scattered by dust
in the accretion disk. On the contrary, we observed a softening of ULX-1 with
decreasing intensity over time. We surmise that this could be due to
reprocessed emission off the accretion disk that increasingly intrudes on our
line-of-sight while the pulsar was obscured. The NICER monitoring campaign
allowed us to observe this long term dimming episode as manifested in the
observed flux (see Fig. 7d).

We showed representative spectra of the bright (red; MJD 58249) and faint
(blue; MJD 58489) states of ULX-1 in Fig. 6. During the 2010 outburst from the
source, Swift/XRT monitoring showed ULX-1 having a very hard spectrum, with
$\Gamma\sim 0$ and a 0.3–10.0 keV unabsorbed flux of $8.6\times
10^{-13}{\rm\,erg\,s^{-1}\,cm^{-2}}$, five times fainter compared to the
brightest state NICER observed around MJD 58239 (Binder et al., 2011b). Four
months after the initial outburst, Chandra observations showed that the flux
dropped by a factor of $\sim 20$ to $4.5\times
10^{-14}{\rm\,erg\,s^{-1}\,cm^{-2}}$ while there was no obvious change in the
spectral shape (within uncertainties for the power law index) (Binder et al.,
2011b). Further Chandra observations (fitting over 0.35–8.0 keV) in 2014
spaced six months apart showed a brightening of the source (by a factor of
$\sim 10$) while there was little evidence of a spectral change, where
$\Gamma=0.0\pm 0.4$ and $\Gamma=1.2\pm 0.8$ in the brighter and fainter
states, respectively (Binder et al., 2016). While we cannot directly compare
the photon indices given the difference in the energy range used in the
spectral fitting, we also note a similar “brightening with accompanying
spectral hardening” behavior in the NICER data (see Fig. 5. As mentioned
above, Carpano et al. (2018) attributed the spectral differences to a large
difference in absorption column density between 2010 (XMM-Newton) and 2016
(XMM-Newton and NuSTAR) observations of ULX-1.

Our findings provided qualitative confirmation of the common interpretation of
the source evolution, where previous work by Carpano et al. (2018) and
Vasilopoulos et al. (2019) have shown that the decreasing flux of ULX-1 can be
explained by increasing obscuration along the line-of-sight. However, in our
work, we opted to fix the interstellar line-of-sight absorption at $0.11\times
10^{22}{\rm\,cm^{-2}}$ (as in Carpano et al. (2018)), and used the photon
index ($\Gamma$) as our proxy to characterize the ULX-1 spectral evolution.
While this is not a physical model of the source, given the low number of
counts in the NICER observations, we would not be able to distinguish between
simple spectral models (e.g., absorbed power law) and more complex spectral
models that include variable absorption and a scattering component (model
parameter degeneracy). We also remark that in our work, we are reporting
results from a one-year monitoring campaign of ULX-1 with amply dense
coverage. While the energy range we employed for the spectral fitting is
narrow (0.4–4.0 keV), we are focusing on the long-term spectral evolution of
the source rather than diagnosing the source state at any given time.

While we did not have sufficient statistics to compare the partial covering
absorber models to the absorbed power law model in a robust manner like in
Carpano et al. (2018), we explored the possibility of short-term variability
of ULX-1 that could be explained by a partial covering absorber. There were
some 5-day binned spectra whose fits were significantly improved by including
an additional thermal disk blackbody component. The significance was assessed
using differences in the AIC values between the absorbed power law and the
absorbed power law and disk blackbody models (see Table 2); we surmise that
this provided some suggestive evidence of the ‘appearance’ of the pulsar in
between obscuration episodes by partial covering absorbers, thus opening up
the detection of X-ray pulsations and thermal emission.

It is intriguing to note from Table 2 that significant fit improvements from
the addition of a disk blackbody traced the detection of X-ray pulsations
somewhat – the additional thermal component is well-detected in some time
intervals prior to MJD 58349, and then pulsations disappeared until MJD 58454,
around which pulsations were detected once again until around MJD 58484
(Vasilopoulos et al., 2019). Thus we suggest a link between the appearance of
a disk blackbody component and the presence of X-ray pulsations. Ray et al.
(2019) also looked at the long-term spin frequency evolution of the ULX pulsar
with the data at MJD 58155, as well as between MJD 58239 and MJD 58411; they
found that the pulsar was spinning up, and stopped detecting significant
pulsations after around MJD 58350. Intriguingly, we found a significant
thermal component in the ULX-1 spectrum around MJD 58504, but no pulsations
seemed to have been detected around that epoch down to a significance of
$4.5\sigma$ (Vasilopoulos et al., 2019).

Optical spectroscopy has identified the companion to be a red supergiant
(RSG), which suggested that the orbital period of the ULX-1 system was at
least 0.8–2.1 yr (Heida et al., 2019). It is noted that NGC 7793 P13 is the
only other known ULX pulsar with a supergiant companion (spectral
classification B9Ia) (Motch et al., 2011, 2014). Townsend & Charles (2020)
derived an empirical relationship between the orbital period ($P_{\rm orb}$)
and superorbital period ($P_{\rm sup}$) using a sample of high mass X-ray
binaries and pulsar ULX systems, where $P_{\rm sup}=(22.9\pm 0.1)P_{\rm orb}$.
If we assume the orbital period for ULX-1 to be 0.8–2.1 yr, then the
superorbital period (i.e., precession period) could be at least 18–41 years,
and could be even higher for larger orbital periods (Heida et al., 2019;
Townsend & Charles, 2020). This means that it could be decades before
sustained pulsation detections from the ULX pulsar will come into full view
again; and we are currently observing a persistent low-lying flux level
punctuated by occasional detections of the source (and pulsations) due to the
partial covering absorber (as mentioned through Table 2).

The observed behavior from ULX-1 suggests that unlike several known transient
ULX pulsars (Tsygankov et al., 2017; Wilson-Hodge et al., 2018; Vasilopoulos
et al., 2020) that show pulsations in outbursts and variable mass accretion
rates, ULX-1 is more like the persistent ULXs (Bachetti et al., 2014; Fürst et
al., 2016; Israel et al., 2017; Sathyaprakash et al., 2019; Rodríguez Castillo
et al., 2020) given the long-term behavior we see, where the mass accretion
rate is constant (Vasilopoulos et al., 2019). This behavior is also similar to
that seen in NGC 7793 P13 (Fürst et al., 2021), where long-term monitoring
(spanning about four years) of the source evolution showed that while the
source X-ray flux was dimming, the period derivative was not correlated with
the X-ray flux, suggesting that a variable accretion rate was not responsible
for the dimming X-ray flux.

Finally, we remark that in the spectrum for the “off-eclipse” in NGC 300 X-1,
we found a Gaussian emission line with centroid energy
$E=0.949_{-0.013}^{+0.012}$ keV and width $\sigma=0.089_{-0.014}^{+0.016}$
keV, which could be associated with Fe-L emission or with the Ne X line
(Degenaar et al., 2013; Bult et al., 2021). It could also be due to Fe XVIII
from photoionized material originating in the outer edge of the accretion
disk, as seen in Cyg X-2 (Vrtilek et al., 1986; Chiappetti et al., 1990), or
the complex blending of Fe and O emission lines (Vrtilek et al., 1988).

We thank the anonymous referee for a careful review of the manuscript. We also
thank Georgios Vasilopoulos for initial discussions that led to refining this
work. M.N., R.R., D.P., and this work were supported by NASA under grant
80NSSC19K1287. This research has made use of data and/or software provided by
the High Energy Astrophysics Science Archive Research Center (HEASARC), which
is a service of the Astrophysics Science Division at NASA/GSFC and the High
Energy Astrophysics Division of the Smithsonian Astrophysical Observatory.

## References

  * Akaike (1974) Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716
  * Arcodia et al. (2018) Arcodia, R., Campana, S., Salvaterra, R., & Ghisellini, G. 2018, A&A, 616, A170, doi: 10.1051/0004-6361/201732322
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bachetti et al. (2014) Bachetti, M., Harrison, F. A., Walton, D. J., et al. 2014, Nature, 514, 202, doi: 10.1038/nature13791
  * Binder et al. (2011a) Binder, B., Williams, B. F., Eracleous, M., et al. 2011a, ApJ, 742, 128, doi: 10.1088/0004-637X/742/2/128
  * Binder et al. (2011b) Binder, B., Williams, B. F., Kong, A. K. H., et al. 2011b, ApJ, 739, L51, doi: 10.1088/2041-8205/739/2/L51
  * Binder et al. (2016) —. 2016, MNRAS, 457, 1636, doi: 10.1093/mnras/stw119
  * Binder et al. (2020) Binder, B. A., Carpano, S., Heida, M., & Lau, R. 2020, Galaxies, 8, 17, doi: 10.3390/galaxies8010017
  * Bult et al. (2021) Bult, P., Strohmayer, T. E., Malacaria, C., Ng, M., & Wadiasingh, Z. 2021, arXiv e-prints, arXiv:2103.12556.  https://arxiv.org/abs/2103.12556
  * Burrows et al. (2005) Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, Space Sci. Rev., 120, 165, doi: 10.1007/s11214-005-5097-2
  * Carpano et al. (2019) Carpano, S., Haberl, F., Crowther, P., & Pollock, A. 2019, IAU Symposium, 346, 187, doi: 10.1017/S1743921318007615
  * Carpano et al. (2018) Carpano, S., Haberl, F., Maitra, C., & Vasilopoulos, G. 2018, MNRAS, 476, L45, doi: 10.1093/mnrasl/sly030
  * Cash (1979) Cash, W. 1979, ApJ, 228, 939, doi: 10.1086/156922
  * Chiappetti et al. (1990) Chiappetti, L., Treves, A., Branduardi-Raymont, G., et al. 1990, ApJ, 361, 596, doi: 10.1086/169223
  * Colbert & Mushotzky (1999) Colbert, E. J. M., & Mushotzky, R. F. 1999, ApJ, 519, 89, doi: 10.1086/307356
  * da Costa-Luis et al. (2022) da Costa-Luis, C., Larroque, S. K., Altendorf, K., et al. 2022, tqdm: A fast, Extensible Progress Bar for Python and CLI, v4.64.0, Zenodo, doi: 10.5281/zenodo.6412640
  * Degenaar et al. (2013) Degenaar, N., Miller, J. M., Wijnands, R., Altamirano, D., & Fabian, A. C. 2013, ApJ, 767, L37, doi: 10.1088/2041-8205/767/2/L37
  * Earnshaw et al. (2019) Earnshaw, H. P., Roberts, T. P., Middleton, M. J., Walton, D. J., & Mateos, S. 2019, MNRAS, 483, 5554, doi: 10.1093/mnras/sty3403
  * Fabrika et al. (2021) Fabrika, S. N., Atapin, K. E., Vinokurov, A. S., & Sholukhova, O. N. 2021, Astrophysical Bulletin, 76, 6, doi: 10.1134/S1990341321010077
  * Fürst et al. (2016) Fürst, F., Walton, D. J., Harrison, F. A., et al. 2016, ApJ, 831, L14, doi: 10.3847/2041-8205/831/2/L14
  * Fürst et al. (2021) Fürst, F., Walton, D. J., Heida, M., et al. 2021, A&A, 651, A75, doi: 10.1051/0004-6361/202140625
  * Gladstone et al. (2009) Gladstone, J. C., Roberts, T. P., & Done, C. 2009, MNRAS, 397, 1836, doi: 10.1111/j.1365-2966.2009.15123.x
  * Heida et al. (2019) Heida, M., Lau, R. M., Davies, B., et al. 2019, ApJ, 883, L34, doi: 10.3847/2041-8213/ab4139
  * Horne & Baliunas (1986) Horne, J. H., & Baliunas, S. L. 1986, ApJ, 302, 757, doi: 10.1086/164037
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Iaria et al. (2018) Iaria, R., Gambino, A. F., Di Salvo, T., et al. 2018, MNRAS, 473, 3490, doi: 10.1093/mnras/stx2529
  * Israel et al. (2017) Israel, G. L., Papitto, A., Esposito, P., et al. 2017, MNRAS, 466, L48, doi: 10.1093/mnrasl/slw218
  * Kaaret et al. (2017) Kaaret, P., Feng, H., & Roberts, T. P. 2017, ARA&A, 55, 303, doi: 10.1146/annurev-astro-091916-055259
  * Kaastra (2017) Kaastra, J. S. 2017, A&A, 605, A51, doi: 10.1051/0004-6361/201629319
  * Kaastra & Bleeker (2016) Kaastra, J. S., & Bleeker, J. A. M. 2016, A&A, 587, A151, doi: 10.1051/0004-6361/201527395
  * Khan et al. (2022) Khan, N., Middleton, M. J., Wiktorowicz, G., et al. 2022, MNRAS, 509, 2493, doi: 10.1093/mnras/stab3049
  * King & Lasota (2016) King, A., & Lasota, J.-P. 2016, MNRAS, 458, L10, doi: 10.1093/mnrasl/slw011
  * King & Lasota (2020) —. 2020, MNRAS, 494, 3611, doi: 10.1093/mnras/staa930
  * Kluzniak & Lasota (2015) Kluzniak, W., & Lasota, J. P. 2015, MNRAS, 448, L43, doi: 10.1093/mnrasl/slu200
  * Koliopanos et al. (2017) Koliopanos, F., Vasilopoulos, G., Godet, O., et al. 2017, A&A, 608, A47, doi: 10.1051/0004-6361/201730922
  * Körding et al. (2002) Körding, E., Falcke, H., & Markoff, S. 2002, A&A, 382, L13, doi: 10.1051/0004-6361:20011776
  * Kosec et al. (2018) Kosec, P., Pinto, C., Walton, D. J., et al. 2018, MNRAS, 479, 3978, doi: 10.1093/mnras/sty1626
  * Kovlakas et al. (2020) Kovlakas, K., Zezas, A., Andrews, J. J., et al. 2020, MNRAS, 498, 4790, doi: 10.1093/mnras/staa2481
  * Liddle (2007) Liddle, A. R. 2007, MNRAS, 377, L74, doi: 10.1111/j.1745-3933.2007.00306.x
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447, doi: 10.1007/BF00648343
  * Middleton et al. (2019) Middleton, M. J., Fragile, P. C., Ingram, A., & Roberts, T. P. 2019, MNRAS, 489, 282, doi: 10.1093/mnras/stz2005
  * Middleton et al. (2018) Middleton, M. J., Fragile, P. C., Bachetti, M., et al. 2018, MNRAS, 475, 154, doi: 10.1093/mnras/stx2986
  * Motch et al. (2011) Motch, C., Pakull, M. W., Grisé, F., & Soria, R. 2011, Astronomische Nachrichten, 332, 367, doi: 10.1002/asna.201011501
  * Motch et al. (2014) Motch, C., Pakull, M. W., Soria, R., Grisé, F., & Pietrzyński, G. 2014, Nature, 514, 198, doi: 10.1038/nature13730
  * Mushtukov et al. (2015a) Mushtukov, A. A., Suleimanov, V. F., Tsygankov, S. S., & Poutanen, J. 2015a, MNRAS, 454, 2539, doi: 10.1093/mnras/stv2087
  * Mushtukov et al. (2015b) —. 2015b, MNRAS, 447, 1847, doi: 10.1093/mnras/stu2484
  * NASA High Energy Astrophysics Science Archive Research Center (2014) (HEASARC) NASA High Energy Astrophysics Science Archive Research Center (HEASARC). 2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, record ascl:1408.004.  http://ascl.net/1408.004
  * Okajima et al. (2016) Okajima, T., Soong, Y., Balsamo, E. R., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, T. Takahashi, & M. Bautz, 99054X, doi: 10.1117/12.2234436
  * Pan et al. (2022) Pan, Y. Y., Li, Z. S., Zhang, C. M., & Zhong, J. X. 2022, MNRAS, 513, 6219, doi: 10.1093/mnras/stac1365
  * Pasham et al. (2015) Pasham, D. R., Cenko, S. B., Zoghbi, A., et al. 2015, ApJ, 811, L11, doi: 10.1088/2041-8205/811/1/L11
  * Pasham et al. (2014) Pasham, D. R., Strohmayer, T. E., & Mushotzky, R. F. 2014, Nature, 513, 74, doi: 10.1038/nature13710
  * Perez & Granger (2007) Perez, F., & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21, doi: 10.1109/MCSE.2007.53
  * Peterson et al. (1998) Peterson, B. M., Wanders, I., Horne, K., et al. 1998, PASP, 110, 660, doi: 10.1086/316177
  * Peterson et al. (2004) Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ, 613, 682, doi: 10.1086/423269
  * Qiu & Feng (2021) Qiu, Y., & Feng, H. 2021, ApJ, 922, 91, doi: 10.3847/1538-4357/ac3569
  * Ray et al. (2019) Ray, P. S., Guillot, S., Ho, W. C. G., et al. 2019, ApJ, 879, 130, doi: 10.3847/1538-4357/ab24d8
  * Remillard et al. (2022) Remillard, R. A., Loewenstein, M., Steiner, J. F., et al. 2022, AJ, 163, 130, doi: 10.3847/1538-3881/ac4ae6
  * Rodríguez Castillo et al. (2020) Rodríguez Castillo, G. A., Israel, G. L., Belfiore, A., et al. 2020, ApJ, 895, 60, doi: 10.3847/1538-4357/ab8a44
  * Sathyaprakash et al. (2019) Sathyaprakash, R., Roberts, T. P., Walton, D. J., et al. 2019, MNRAS, 488, L35, doi: 10.1093/mnrasl/slz086
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835, doi: 10.1086/160554
  * Standish (1998) Standish, E. M. 1998, JPL Planetary and Lunar Ephemerides, DE405/LE405, JPL Interoffice Memo 312.F-98-048 (Pasadena: NASA Jet Propulsion Laboratory)
  * Steiner et al. (2009) Steiner, J. F., Narayan, R., McClintock, J. E., & Ebisawa, K. 2009, PASP, 121, 1279, doi: 10.1086/648535
  * Strohmayer & Mushotzky (2009) Strohmayer, T. E., & Mushotzky, R. F. 2009, ApJ, 703, 1386, doi: 10.1088/0004-637X/703/2/1386
  * Sutton et al. (2013) Sutton, A. D., Roberts, T. P., & Middleton, M. J. 2013, MNRAS, 435, 1758, doi: 10.1093/mnras/stt1419
  * Tan & Biswas (2012) Tan, M. Y. J., & Biswas, R. 2012, MNRAS, 419, 3292, doi: 10.1111/j.1365-2966.2011.19969.x
  * Townsend & Charles (2020) Townsend, L. J., & Charles, P. A. 2020, MNRAS, 495, 139, doi: 10.1093/mnrasl/slaa078
  * Tsygankov et al. (2017) Tsygankov, S. S., Doroshenko, V., Lutovinov, A. A., Mushtukov, A. A., & Poutanen, J. 2017, A&A, 605, A39, doi: 10.1051/0004-6361/201730553
  * Vasilopoulos et al. (2018) Vasilopoulos, G., Haberl, F., Carpano, S., & Maitra, C. 2018, A&A, 620, L12, doi: 10.1051/0004-6361/201833442
  * Vasilopoulos et al. (2019) Vasilopoulos, G., Petropoulou, M., Koliopanos, F., et al. 2019, MNRAS, 488, 5225, doi: 10.1093/mnras/stz2045
  * Vasilopoulos et al. (2020) Vasilopoulos, G., Ray, P. S., Gendreau, K. C., et al. 2020, MNRAS, 494, 5350, doi: 10.1093/mnras/staa991
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Vrtilek et al. (1986) Vrtilek, S. D., Kahn, S. M., Grindlay, J. E., Helfand, D. J., & Seward, F. D. 1986, ApJ, 307, 698, doi: 10.1086/164455
  * Vrtilek et al. (1988) Vrtilek, S. D., Swank, J. H., & Kallman, T. R. 1988, ApJ, 326, 186, doi: 10.1086/166079
  * Wachter et al. (2000) Wachter, S., Smale, A. P., & Bailyn, C. 2000, ApJ, 534, 367, doi: 10.1086/308754
  * Walton et al. (2022) Walton, D. J., Mackenzie, A. D. A., Gully, H., et al. 2022, MNRAS, 509, 1587, doi: 10.1093/mnras/stab3001
  * Walton et al. (2018a) Walton, D. J., Fürst, F., Heida, M., et al. 2018a, ApJ, 856, 128, doi: 10.3847/1538-4357/aab610
  * Walton et al. (2018b) Walton, D. J., Bachetti, M., Fürst, F., et al. 2018b, ApJ, 857, L3, doi: 10.3847/2041-8213/aabadc
  * Welsh (1999) Welsh, W. F. 1999, PASP, 111, 1347, doi: 10.1086/316457
  * Wilson-Hodge et al. (2018) Wilson-Hodge, C. A., Malacaria, C., Jenke, P. A., et al. 2018, ApJ, 863, 9, doi: 10.3847/1538-4357/aace60