# Hybrid SBI or How I Learned to Stop Worrying and Learn the Likelihood

Chirag Modi<EMAIL_ADDRESS>Center for Computational
Astrophysics, Flatiron Institute, New York, NY 10010, USA  Center for
Computational Mathematics, Flatiron Institute, New York, NY 10010, USA  Oliver
H. E. Philcox<EMAIL_ADDRESS>Simons Society of Fellows, Simons
Foundation, New York, NY 10010, USA  Center for Theoretical Physics, Columbia
University, New York, NY 10027, USA

###### Abstract

We propose a new framework for the analysis of current and future cosmological
surveys, which combines perturbative methods (PT) on large scales with
conditional simulation-based implicit inference (SBI) on small scales. This
enables modeling of a wide range of statistics across all scales using only
small-volume simulations, drastically reducing computational costs, and avoids
the assumption of an explicit small-scale likelihood. As a proof-of-principle
for this hybrid simulation-based inference (HySBI) approach, we apply it to
dark matter density fields and constrain cosmological parameters using both
the power spectrum and wavelet coefficients, finding promising results that
significantly outperform classical PT methods. We additionally lay out a
roadmap for the next steps necessary to implement HySBI on actual survey data,
including consideration of bias, systematics, and customized simulations. Our
approach provides a realistic way to scale SBI to future survey volumes,
avoiding prohibitive computational costs.

Current and upcoming cosmological surveys such as the Dark Energy
Spectroscopic Instrument [DESI; 1, 2, 3], the ESA Euclid satellite mission
[4], Vera C. Rubin Observatory (LSST) [5], and the NASA Nancy Grace Roman
Space Telescope [Roman; 6, 7], will probe the Universe across a vast array of
scales providing high-resolution measurements with the potential to constrain
cosmological parameters to unprecedented precision. The past iteration of
cosmological data has been analyzed primarily through traditional analytical
methods such as perturbation theory (PT); whilst this has been highly
successful [8, 9, 10], the approach is necessarily limited to low-order
clustering statistics (such as the two- and three-point functions), and to
linear and quasi-linear scales. If we wish to fully exploit the information
contained within upcoming large scale structure surveys, new methods are
needed.

Recently, simulation-based inference (SBI) has emerged as a promising
alternative to overcome the limitations of traditional methods [11, 12, 13,
14, 15]. SBI uses numerical simulations to model observables and their summary
statistics, before these are combined with neural density estimation methods
to efficiently infer cosmological parameters. As a result, it can be applied
to any summary statistic that can be simulated and evaluated, including those
deep in the non-linear regime, inaccessible to PT, and those whose likelihoods
are far from Gaussian. A major limitation of SBI is that its training requires
a large number of high-fidelity simulations, spanning a wide range of
cosmologies. This makes it a computationally expensive exercise; a drawback
that is further exacerbated by the increasing simulation volumes and
resolutions needed to model upcoming data. The largest simulation suite
currently available for training SBI for galaxy clustering (quijote
simulations [16]) has moderately low mass resolution and is only
$1\,h^{-3}\mathrm{Gpc}^{3}$ in volume, which is smaller than than the SDSS
III-BOSS survey [17]. Creating analogous simulations appropriate for the next
generation of cosmological surveys is thus likely to be _computationally
prohibitive_.

The primary challenge in scaling up simulations stems from the simultaneous
requirements of large volumes, high resolution, and accurate small-scale
dynamics. On the largest scales, however, physics is well understood
analytically, thus simulation-based approaches are not strictly necessary. In
this regime, the density field is close to Gaussian and with quasi-linear
perturbations that can be well modeled using perturbation theory. Furthermore,
traditional statistics like the power spectrum and bispectrum are close to
optimal on such scales [18], thus any more nuanced analysis schemes will not
yield significant additional information. It is only the small, non-linear
scales where higher-order statistics dominate and thus SBI is needed. The
conclusion is clear: if we can combine perturbative methods for analyzing
large scales with simulation-based methods for analyzing small scales, we can
mitigate the computational cost, and thus make feasible the full analysis of
upcoming generations of surveys.

In this work, we propose such a framework for hybrid simulation-based
inference (HySBI), and demonstrate its utility in predicting cosmological
parameters from summary statistics of the matter density field. This is
depicted schematically in Fig. 1.

Figure 1: Schematic overview of this work. Large-scale statistics are modeled
using classical perturbative techniques, whilst small-scale summaries make use
of simulation-based inference, trained on small sub-volume simulations,
obviating the need to run costly high-resolution simulations of the entire
survey volume.

##  I Simulation-based Inference

SBI begins with a training dataset composed of ($\bm{\theta},\bm{x}$) pairs
where $\bm{\theta}$ are the underlying model parameters to be inferred (here,
cosmology parameters and any relevant nuisance parameters), and $\bm{x}$ is
the corresponding data-vector generated from a simulator (such as an $N$-body
realization). Analysis usually proceeds by training a conditional neural
density estimator $q_{\phi}$ (with parameters $\phi$) to learn: (a) the
likelihood model, $p(\bm{x}|\bm{\theta})$, by maximizing the conditional log-
probability of the data $\bm{x}$ conditioned on the parameters $\bm{\theta}$
(known as a neural likelihood estimator; NLE); or (b) the posterior itself,
$p(\bm{\theta}|\bm{x})$, by maximizing the conditional log-probability of the
model parameters $\bm{\theta}$ conditioned on the data $\bm{x}$, over the
training data (known as a neural posterior estimator; NPE).

When performing inference on a test observation $\bm{x}^{\prime}$ when using
NPE, we can directly query the trained estimator $q_{\phi^{*}}$ to generate
samples from the posterior, i.e $\bm{\theta}\sim
q_{\phi^{*}}(\bm{\theta}|\bm{x}^{\prime})$. In NLE, we instead need to combine
the learnt likelihood model $q_{\phi^{*}}(\bm{x}|\bm{\theta})$ with the prior
$p(\bm{\theta})$, then use Markov chain Monte Carlo (MCMC) algorithms to
generate samples from this posterior. NPE is often preferred over NLE for two
reasons: (1) often the dimensionality of the posterior distribution is smaller
than that of the likelihood which can lead to asymptotically better learning;
(2) NPE does not require the additional step of running MCMC. In the HySBI
framework however, we use NLE, since it enables us to combine the learnt
neural likelihood on small scales with the analytic likelihood on large
scales.

##  II Hybrid Simulation-based Inference

Consider the scenario in which the data-vector $\bm{x}$ can be split into two
components – large scales $\bm{x}_{L}$, and small scales $\bm{x}_{S}$ – such
that $\bm{x}=\\{\bm{x}_{L},\bm{x}_{S}\\}$. Such a split is natural for the
power spectrum (where $\bm{x}_{L}$ would comprise measurements with
wavenumbers below some transition scale, i.e. $\\{P(k):k\leq k_{*}\\}$);
however, we note that the following discussion is applicable to any
combination of large- and small-scale statistics (e.g., power spectra on large
scales and wavelets on small scales). The data likelihood can then be
decomposed as the product of a large-scale likelihood and the conditional
likelihood of small-scales given large-scales:

$p(\bm{x}|\bm{\theta})=p(\bm{x}_{L},\bm{x}_{S}|\bm{\theta})=p(\bm{x}_{L}|\bm{\theta})p(\bm{x}_{S}|\bm{x}_{L},\bm{\theta}).$
(1)

In HySBI, we propose to use classical statistics such as the power spectrum
and bispectrum on large scales for which the likelihood,
$p(\bm{x}_{L}|\bm{\theta})$, can be modeled analytically with perturbation
theory [e.g., 19, 8, 9]. As such, only the likelihood term on small scales,
$p(\bm{x}_{S}|\bm{x}_{L},\bm{\theta})$, needs to be learnt with
simulations.111This differs from previous hybrid approaches, such as Ref. [20,
21], who use small-volume simulations to model the signal, but assume a
Gaussian likelihood; this fails for many higher-order statistics. We stress
that, in contrast to global SBI approaches, this likelihood term depends not
only on the model parameters $\bm{\theta}$ but also on the large scale
statistic $\bm{x}_{L}$.

As discussed above, the main advantage of learning only the small-scale
likelihood with simulations is that, in principle, it can be done by
simulating only a small fractional volume at full fidelity instead of the
whole simulation box, creating a surrogate likelihood that mitigates the
increasing computational cost of scaling SBI to upcoming cosmological surveys.
For a sufficiently large sub-volume, we still have enough modes on small
scales that the associated sample variance is sub-dominant compared to other
errors such as data noise and stochasticity in training. However, two new
issues arise: (1) to correctly learn the conditional dependence
$p(\bm{x}_{S}|\bm{x}_{L},\bm{\theta})$, one needs to be able to estimate the
correct large-scale statistics $\bm{x}_{L}$ for the exact large-scale
realization in which the small-scale Universe forms with only coarse
simulations on the full volume and full fidelity simulations only on the small
volume; (2) a given sub-volume suffers from increased sample variance and
super-sample effects from the large box that will impact the evolution in a
given sub-volume, thus adding noise to the small scale statistics $\bm{x}_{S}$
used for training SBI. Thus we learn a _surrogate likelihood_ for the small
scale data which has the same mean but larger volume and hence is more
conservative (effectively trading bias from poor quality simulations with
variance from a smaller simulation box). In this first work, we present a
proof-of-principle of our methodology and focus on studying the impact of this
super-sample variance. In the discussion below, we will comment on other
subtleties associated with this approach.

##  III Data Analysis

### Setup –

In order to present our methodology and its advantages, this work considers
constraining the $\Omega_{m}$ and $\sigma_{8}$ parameters from the three-
dimensional dark matter distribution. Our motivation to use dark matter
density field rather than a more realistic halo or galaxy sample for this
first exercise is to avoid the need to infer non-cosmological parameters (such
as those of the halo-occupation distribution), allowing us to focus solely on
the impact of super-sample variance.

To train the small-scale likelihood we use the quijote latin-hypercube suite
[16] which contains 2000 $N$-body simulations varying five cosmology
parameters: $\\{\Omega_{\mathrm{m}},\,\sigma_{8},\,\Omega_{\rm
b},\,n_{s},h\\}$. These simulations evolve $1024^{3}$ cold dark matter
particles in a $1\,h^{-3}\mathrm{Gpc}^{3}$ box with the TreePM Gadget-III code
[22]. To mimic a small-volume simulation, we split each simulation into eight
sub-volumes with side-length $500\,h^{-1}\mathrm{Mpc}$. 1500 simulations are
used for training, 200 for validation, and the remaining 300 for testing. We
emphasize that we use the sub-volumes to estimate small scale statistics only
for _training HySBI_ and in tests, we always analyze the statistics estimated
from the whole-volume as simulated data.

### Large-scale statistics –

Here, we set $\bm{x}_{L}$ equal to the observed power spectrum of dark matter
density field, $\hat{P}(k)$, on scales $k\in[0.007,0.15]\,h\mathrm{Mpc}^{-1}$,
obtained from the full simulation volume of $1\,h^{-3}\mathrm{Gpc}^{3}$ with
$\Delta k=2\pi/1000\,h\,\mathrm{Mpc}^{-1}$, i.e. the fundamental mode of the
box. This scale cut is chosen to balance two competing factors: (a)
perturbation theory predictions are accurate only on large-scales (favoring
lower $k_{\mathrm{max}}$); (b) to suppress sample-variance in the small-scale
sub-box statistics, we require many modes and thus larger $k_{\mathrm{max}}$.
Power spectra are estimated with nbodykit, using cloud-in-cell interpolation
on an $N=256$ side mesh [23].

We adopt a Gaussian model for the large-scale likelihood, as specified by the
Effective Field Theory of Large Scale Structure:

$\displaystyle-2\log p(\bm{x}_{L}|\bm{\theta})$ $\displaystyle=$
$\displaystyle\sum_{k}\left[\frac{P_{\rm loop}(k)-2c_{s}^{2}P_{\rm
ct}(k)-\hat{P}(k)}{\sigma_{P}(k)}\right]^{2},$ (2)

where the variance term is specified by the linear theory power spectrum,
evaluated at the true cosmology of the test simulations. Our model for the
theoretical power spectrum is given by 1-loop standard perturbation theory
($P_{\rm loop}(k)$) in addition to a renormalization contribution
$-2c_{s}^{2}P_{\rm ct}(k)$, depending on the counter-term $c_{s}^{2}$, which
must be treated as a free parameter and marginalized over [e.g., 24,
25].222The PT models are computed via loop integrals over the linear power
spectrum (with e.g., class-pt [26]), itself computed using Boltzmann solvers
such as class or camb. For the sake of speed, here we instead train a simple
neural network emulator to predict $P_{\mathrm{loop}}$ and $P_{\mathrm{ct}}$
as a function of cosmological parameters (given a grid of training data from
class-pt). We have validated that this does not impact the final inference but
speeds up our MCMC chains by a factor $\approx 100$.

### Small scale statistics –

We consider two different choices for the small-scale statistic $\bm{x}_{S}$:
(1) the power spectrum for $k\leq 0.5\,h\,\mathrm{Mpc}^{-1}$; (2) wavelet
coefficients.333Since the wavelets are not localized in Fourier space, it is
hard to define a strict cut-off. The coefficients used here are constructed
using a mother wavelet whose weight falls to zero at the Nyquist frequency
$k=0.8\,h\,\mathrm{Mpc}^{-1}$. In both cases, the conditional likelihood model
$p(\bm{x}_{S}|\bm{x}_{L},\bm{\theta})$, is learnt using SBI, given the above
$\bm{x}_{L}$ definition. To train the likelihood model, we estimate the small-
scale statistics using only the sub-volumes of the simulation, one at a time;
this emulates the realistic setting when it is practical only to compute
reduced-volume simulations. We thus have eight estimates of $\bm{x}_{S}$ from
each simulation which generically differ due to sample and super-sample
variance effects. To test our methodology, however, we generate $\bm{x}_{S}$
from the full simulation volume of $1\,h^{-3}\mathrm{Gpc}^{3}$ such that all
the observational data is used.

Here, we measure the power spectrum on scales
$k\in[0.15,0.5]\,h\,\mathrm{Mpc}^{-1}$ with $\Delta k$ fixed to the
fundamental mode of the small sub-volume, i.e.
$2\pi/500\,h\,\mathrm{Mpc}^{-1}$; this allows for sub-volume simulations and
full-volume ‘data’ to be robustly compared. Wavelet coefficients are defined
for a statistically homogeneous real field $F$ as
$S_{0}=\langle|F|^{p}\rangle$ and $S_{1}=\langle|I\ast\psi_{j}|^{p}\rangle$
where the dilation index $j$ varies as $q/Q$ with $q\in\\{0,\dots,J\times
Q-1\\}$. The parameter $J$ is the total number of octaves i.e. doublings in
scale of the mother wavelet, and $Q$ is the quality factor, which controls the
number of scales per octave. We use the package galactic-
wavelets444https://github.com/bregaldo/galactic_wavelets to estimate these
wavelet coefficients. Our mother wavelet is inspired from [27] and we take
$J=3$, $Q=4$, and $p\in\\{0.5,1,1.5,2\\}$, resulting in a total of four
$S_{0}$ and 48 $S_{1}$ wavelet coefficients [28].

### Non-periodic boundary conditions –

Being cut-outs from a larger simulation, the sub-volumes do not have periodic
boundary conditions, which must be accounted for when estimating the small-
scale power spectrum and wavelet coefficients. For power spectrum, this is
achieved via a (square) window-function matrix $M$ such that the true small
scale power is $P(k)=[M^{-1}\tilde{P}](k)$, where $\tilde{P}$ is the measured
small-scale power spectrum on the volume, assuming periodicity. This matrix is
estimated by generating Gaussian white-noise fields over the full-volume with
power injected only in a single $k-$bin, and then measuring the leakage of
this power into other bins when power spectrum is evaluated on a sub-volume.

Performing similar corrections for wavelet coefficients is easier since they
are localized in configuration space. Here, we achieve this by masking the
pixels near the boundary such that for any wavelet centered on a pixel within
the mask, the fraction of weight that leaks outside the sub-volume is less
than a pre-defined small threshold. This implies that larger the wavelets, the
more pixels are masked. As such, the threshold determines the bias-variance
trade-off in our estimator, and is here set to 0.02.

### SBI Implementation –

We use the sbi555https://github.com/mackelab/sbi package to train masked auto-
regressive flows as conditional neural density estimators and learn the
likelihood, $q_{\phi}(\bm{x}_{S}|\bm{x}_{L},\bm{\theta})\sim
p(\bm{x}_{S}|\bm{x}_{L},\bm{\theta})$. To minimize stochasticity in training,
we train 400 networks for each statistic by varying hyperparameters such as
the network width, number of layers, learning rate and batch size. We use the
Weights-and-Biases666https://wandb.ai/site package for the hyperparameter
search. After training, we collect the ten neural density estimators with best
validation loss and use them as an ensemble.

### Posterior Inference –

For inference over a new observation $\bm{x}^{\prime}$, we combine the learnt
small-scale likelihood with the Gaussian likelihood for large scale
$p(\bm{x}_{L}|\bm{\theta})$ (via Eq. 1), and a uniform prior on cosmology
parameters (dictated by the bounds of the quijote latin-hypercube). The
counter-term $c_{s}$ has a Gaussian prior $\mathcal{N}(0,10)$ [e.g., 8]. To
generate samples from the posterior, we use the emcee package [29] for MCMC
sampling and run 20 chains with 11,000 samples each.

##  IV Results

Figure 2: Posterior distribution of $\Omega_{m}$ and $\sigma_{8}$ inferred
from a power spectrum analysis using classical perturbative methods on large-
scales (PT), full simulation-based inference (NLE), and our new hybrid scheme
(HySBI). We also show HySBI results including small-scale wavelet statistics
(purple). Results are shown for a representative test simulation (quiojte
LH-872). HySBI strongly outperforms classical techniques, and is able to
robustly infer cosmological parameters without large-volume simulations.

In Fig. 2, we show the results of HySBI using the power spectrum and wavelet
coefficients for inferring $\Omega_{m}$ and $\sigma_{8}$ on one of the test
simulations777Results are similar on all the 160 simulations in the test
dataset.. To minimize the impact of super-sample variance, the SBI component
of HySBI is trained using all the eight sub-volumes of the simulations i.e. we
take the mean of the small-scale statistics estimated from all the individual
sub-volumes for training.888The mean of all subvolume statistics well
approximates the statistic computed on the total volume; differences appear
only from boundary effects, which are small if the characteristic scale of
$\bm{x}$ is small compared to the subbox size. As discussed in the previous
section, the test data is generated from the full $1\,h^{-3}\mathrm{Gpc}^{3}$
simulation volume.

For comparison, we show also the results of the PT analysis which is
restricted to the large scales ($k\in[0.007,0.15]\,h^{-1}\mathrm{Mpc}$), and a
global SBI power spectrum analysis (hereafter NLE) that covers the same scales
as HySBI ($k\in[0.007,0.5]\,h^{-1}\mathrm{Mpc}$). As expected, HySBI
significantly outperforms the classical PT analysis, clearly demonstrating the
utility of our method. That said, the NLE power spectrum constraints on
$\sigma_{8}$ are somewhat tighter than those of HySBI over the same range;
this is because, in our scale-decoupled formalism, HySBI also infers and
marginalizes over the large-scale PT counter-term $c_{s}^{2}$, unlike global
SBI. This can in principle be improved upon by learning a prior of $c_{s}^{2}$
conditioned on small scale statistics and cosmology, but we do not explore
that in this work. Finally, HySBI with wavelet coefficients on small scales
outperforms all the power spectrum analyses, emphasizing the advantage of
using higher order statistics.

Figure 3:  As Fig. 2, but comparing the power spectrum HySBI constraints when
the SBI is trained using one, two, four, and eight sub-volumes of every
simulation. This reduces super-sample variance in the small-scale training
data, and is found to somewhat enhance $\sigma_{8}$ constraints. Figure 4: As
Fig. 3, but using wavelet coefficients as the small-scale statistic. The
improvement in $\sigma_{8}$ with number of sub-volumes is more muted in this
case.

To study the impact of super-sample covariance, we repeat the HySBI analysis
training the small-scale likelihood using different number of sub-volumes
(each of size $500^{3}\,h^{-3}\mathrm{Mpc}^{3}$). The corresponding
constraints on $\Omega_{m}$ and $\sigma_{8}$ are shown in Fig. 3 for a single
representative simulation using the power spectrum as the summary statistic
for both small and large scales. Averaging over all test simulations, the
uncertainties when using four, two and one sub-volumes are inflated by $\sim
6,\,9,\,11\%$ ($\sim 45,\,80,\,130\%$) for $\Omega_{m}$ ($\sigma_{8}$)
compared to those using all eight sub-volumes. The inference of $\sigma_{8}$
is thus more sensitive to super-sample effects than $\Omega_{m}$; this matches
expectations, since most of the constraining power on the latter comes from
the large scales which are modeled with PT, whilst $\sigma_{8}$ constraints
are most informed by small scales.

In Fig. 4, we show the same analysis using wavelet coefficients on small
scales. In this case, we find that constraints are inflated by $\sim
20,\,24,\,50\%$ ($\sim 40,\,80,\,100\%$) for $\Omega_{m}$ ($\sigma_{8}$); here
the $\Omega_{m}$ analysis more sensitive to super-sample effects than before,
since the wavelets carry significant information on the matter density on all
scales. On the other hand, $\sigma_{8}$ constraints are slightly less
sensitive for wavelet coefficients than the power spectrum since wavelets
include statistics corresponding to $\delta^{p}$ for $p=0.5,1,1.5$ (where
$\delta$ is the overdensity field), whilst the power spectrum only measures
squared statistics, i.e. $\delta^{2}$. The former lead to lower sample
variance across the different sub-volumes.

##  V Discussion and Outlook

Simulation-based inference provides a promising avenue through which to
overcome the limitations of traditional analytic methods for cosmological
analysis, and make use of the upcoming tranche of small-scale data. However,
the computational cost of generating the high-fidelity and high-volume
simulations needed to train SBI at the precision required for the next
generation of cosmological surveys will likely be prohibitive.

To mitigate this challenge, we propose Hybrid simulation-based inference
(HySBI); this combines perturbative methods on large scales (using quasi-
optimal statistics such as the power spectrum and bispectrum), with
simulation-based techniques on small scales (including arbitrary higher-order
statistics which cannot be modeled analytically). This is achieved by
decomposing the data likelihood as the product of a large-scale likelihood
modeled analytically with perturbation theory, and a conditional likelihood of
small-scales given large-scales – the latter is learnt using simulations. This
approach can be lead to significant computational gains, since we need to
simulate only a small fraction of the whole data volume to model the whole.

A natural drawback is that modeling small-scales using only one (or a few)
sub-volumes of a simulation adds super-sample effects from the large scale
modes in the global simulation box. As a result, we learn a conservative
surrogate likelihood for small scales which has a larger variance (but removes
the bias that would result from simulating the entire volume at low
resolution). In this first work, we have presented a proof-of-principle of
HySBI and studied the impact of this super-sample variance by inferring
$\sigma_{8}$ and $\Omega_{m}$ using the three-dimensional distribution of the
dark-matter density field, analyzed with the large-scale power spectrum, and
small-scale power spectrum or wavelet coefficients. Upon splitting the entire
$1h^{-3}\mathrm{Gpc}^{3}$ simulation box into eight sub-volumes, we found that
using even one of the eight sub-volumes to train the small-scale likelihood
leads to significant gains over state-of-the-art power spectrum analyses.

Performing a global SBI analysis leads to slightly better results than our
hybrid scheme; this occurs since we have enforced that the large-scale
likelihood is not influenced by small-scales, and thus marginalized over a
necessary perturbative counterterm parameter appearing on large scales. Such
an assumption must be altered for realistic observables such as galaxies,
whence there are nuisance parameters in both PT and SBI (encoding bias and
halo-occupation distributions). These parameters however are not independent
and should be inferred jointly. This will be explored in future work.

To fully realize the computational advanced promised by HySBI, two technical
advances must be realized. As mentioned briefly above, a key ingredient of our
SBI framework is accurate, high-fidelity small sub-volume training simulations
which “know” about the large-volume Universe in which they are embedded, i.e.
those from which we can measure both the small-scale statistics, $\bm{x}_{S}$,
and the corresponding realization of the large-scale data, $\bm{x}_{L}$ but
without running the full-fidelity simulations on the full volume for the
latter. We plan on developing these simulations by building upon promising
techniques that are already available. For example, the s-COLA method creates
coarse PT simulations on large-scales (which provide the global $\bm{x}_{L}$
measurements), and evolves small sub-volumes with high-resolution particle-
mesh simulations [30, 31]. The FlowPM simulations, which use multi-grid
schemes to estimate small-scale forces independently in different sub-volumes,
can also be adapted to this end [32]. A third avenue is by building upon the
zoom-in simulations for galaxy formation [33]. Crucially, sub-volumes can be
independently simulated in these frameworks, thus we can straightforwardly run
a single small-volume simulation with a known (coarsened) large-volume
realization. Since these can be run on a single GPU, the resulting pipeline
will not require multi-node communications and can enjoy GPU accelerations;
this will lead to significant computational gains even if we simulate all the
sub-volumes of a given simulation to minimize super-sample effects [34, 32].

A second technical hurdle to be overcome is the the inclusion of systematic
effects such as survey masks, fibre collisions et cetera. These often couple
small- and large-scales and will be considered in future works. Finally, we
currently assume a clear division between the large and small scale statistics
(such that $\bm{x}_{S}$ can be measured from sub-boxes), which excludes some
types of statistics such as squeezed bispectra; similar techniques can likely
be developed for such modified analyses (often of use in inflationary studies)
[e.g., 35, 36, 37, 38].Though much work remains to be done, it seems likely
that hybrid approaches such as HySBI can create a promising symbiosis of old
and new techniques, significantly enhancing the science output of future
surveys without requiring prohibitive computational costs.

###### Acknowledgements.

OHEP is a Junior Fellow of the Simons Society of Fellows and thanks the office
of Chuck Schumer for immigration assistance. CM would like to thank Dr.
Strangelove for the title, and Ben Wandelt, Michael Eickenberg, and Mikhail
Ivanov for useful discussions. We would especially like to thank Bruno
Régaldo-Saint Blancard for setting up and sharing the code for estimating
wavelet coefficients with erosion. We additionally thank Sebastian Wagner-
Carena for insightful comments on the draft. This work was partly conceived
during the “Beyond 2-point Challenge” workshop at the University of Arizona in
February 2023.

## References

  * Collaboration  _et al._ [2016a] D. Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L. E. Allen, C. A. Prieto, J. Annis, S. Bailey, C. Balland, _et al._ , The DESI Experiment Part I: Science,Targeting, and Survey Design, arXiv:1611.00036 [astro-ph] (2016a), arXiv:1611.00036 [astro-ph] .
  * Collaboration  _et al._ [2016b] D. Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L. E. Allen, C. A. Prieto, J. Annis, S. Bailey, C. Balland, O. Ballester, _et al._ , The DESI Experiment Part II: Instrument Design, arXiv:1611.00037 [astro-ph] (2016b), arXiv:1611.00037 [astro-ph] .
  * Abareshi  _et al._ [2022] B. Abareshi, J. Aguilar, S. Ahlen, S. Alam, D. M. Alexander, R. Alfarsy, L. Allen, C. A. Prieto, _et al._ , Overview of the Instrumentation for the Dark Energy Spectroscopic Instrument (2022), arXiv:2205.10939 [astro-ph] .
  * Laureijs  _et al._ [2011] R. Laureijs, J. Amiaux, S. Arduini, J.-L. Auguères, J. Brinchmann, R. Cole, M. Cropper, C. Dabin, L. Duvet, A. Ealet, _et al._ , Euclid Definition Study Report, arXiv e-prints , arXiv:1110.3193 (2011).
  * Ivezić  _et al._ [2019] Ž. Ivezić, S. M. Kahn, J. A. Tyson, B. Abel, E. Acosta, R. Allsman, D. Alonso, Y. AlSayyad, S. F. Anderson, J. Andrew, and et al., LSST: From Science Drivers to Reference Design and Anticipated Data Products, ApJ 873, 111 (2019), arXiv:0805.2366 .
  * Spergel  _et al._ [2015] D. Spergel, N. Gehrels, C. Baltay, D. Bennett, J. Breckinridge, M. Donahue, A. Dressler, B. S. Gaudi, T. Greene, O. Guyon, _et al._ , Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 Report (2015).
  * Wang  _et al._ [2022] Y. Wang, Z. Zhai, A. Alavi, E. Massara, A. Pisani, A. Benson, C. M. Hirata, L. Samushia, D. H. Weinberg, _et al._ , The High Latitude Spectroscopic Survey on the Nancy Grace Roman Space Telescope, The Astrophysical Journal 928, 1 (2022).
  * Philcox and Ivanov [2022] O. H. E. Philcox and M. M. Ivanov, BOSS DR12 full-shape cosmology: $\Lambda$CDM constraints from the large-scale galaxy power spectrum and bispectrum monopole, Phys. Rev. D 105, 043517 (2022), arXiv:2112.04515 [astro-ph.CO] .
  * D’Amico  _et al._ [2022a] G. D’Amico, Y. Donath, M. Lewandowski, L. Senatore, and P. Zhang, The BOSS bispectrum analysis at one loop from the Effective Field Theory of Large-Scale Structure, (2022a), arXiv:2206.08327 [astro-ph.CO] .
  * Chen  _et al._ [2022] S.-F. Chen, Z. Vlah, and M. White, A new analysis of galaxy 2-point functions in the BOSS survey, including full-shape information and post-reconstruction BAO, JCAP 02 (02), 008, arXiv:2110.05530 [astro-ph.CO] .
  * Alsing  _et al._ [2018] J. Alsing, B. Wandelt, and S. Feeney, Massive optimal data compression and density estimation for scalable, likelihood-free inference in cosmology, arXiv:1801.01497 [astro-ph] (2018), arXiv:1801.01497 [astro-ph] .
  * Alsing  _et al._ [2019] J. Alsing, T. Charnock, S. Feeney, and B. Wandelt, Fast likelihood-free cosmology with neural density estimators and active learning, Monthly Notices of the Royal Astronomical Society 488, 4440 (2019).
  * Jeffrey  _et al._ [2021] N. Jeffrey, J. Alsing, and F. Lanusse, Likelihood-free inference with neural compression of DES SV weak lensing map statistics, Monthly Notices of the Royal Astronomical Society 501, 954 (2021).
  * Hahn  _et al._ [2022] C. Hahn, M. Eickenberg, S. Ho, J. Hou, P. Lemos, E. Massara, A. Modi Chirag, Moradinezhad Dizgah, B. Régaldo-Saint Blancard, and M. M. Abidi, SimBIG: A Forward Modeling Approach To Analyzing Galaxy Clustering, (2022).
  * Cranmer  _et al._ [2020] K. Cranmer, J. Brehmer, and G. Louppe, The frontier of simulation-based inference, Proceedings of the National Academy of Sciences 117, 30055 (2020), https://www.pnas.org/content/117/48/30055.full.pdf .
  * Villaescusa-Navarro  _et al._ [2020] F. Villaescusa-Navarro, C. Hahn, E. Massara, A. Banerjee, A. M. Delgado, D. K. Ramanah, T. Charnock, E. Giusarma, Y. Li, E. Allys, _et al._ , The Quijote Simulations, The Astrophysical Journal Supplement Series 250, 2 (2020).
  * Dawson  _et al._ [2013] K. S. Dawson, D. J. Schlegel, C. P. Ahn, S. F. Anderson, É. Aubourg, S. Bailey, R. H. Barkhouser, J. E. Bautista, A. Beifiori, A. A. Berlind, _et al._ , The Baryon Oscillation Spectroscopic Survey of SDSS-III, The Astronomical Journal 145, 10 (2013).
  * Cabass  _et al._ [2023] G. Cabass, M. Simonović, and M. Zaldarriaga, Cosmological Information in Perturbative Forward Modeling, (2023), arXiv:2307.04706 [astro-ph.CO] .
  * Ivanov  _et al._ [2020] M. M. Ivanov, M. Simonović, and M. Zaldarriaga, Cosmological parameters from the BOSS galaxy power spectrum, Journal of Cosmology and Astroparticle Physics 2020, 042.
  * Kokron  _et al._ [2021] N. Kokron, J. DeRose, S.-F. Chen, M. White, and R. H. Wechsler, The cosmology dependence of galaxy clustering and lensing from a hybrid N-body–perturbation theory model, Mon. Not. Roy. Astron. Soc. 505, 1422 (2021), arXiv:2101.11014 [astro-ph.CO] .
  * Modi  _et al._ [2020] C. Modi, S.-F. Chen, and M. White, Simulations and symmetries, MNRAS 492, 5754 (2020), arXiv:1910.07097 [astro-ph.CO] .
  * Springel [2005] V. Springel, The Cosmological simulation code GADGET-2, Mon. Not. Roy. Astron. Soc. 364, 1105 (2005), arXiv:astro-ph/0505010 .
  * Hand  _et al._ [2018] N. Hand, Y. Feng, F. Beutler, Y. Li, C. Modi, U. Seljak, and Z. Slepian, Nbodykit: An Open-source, Massively Parallel Toolkit for Large-scale Structure, The Astronomical Journal 156, 160 (2018).
  * Baumann  _et al._ [2012] D. Baumann, A. Nicolis, L. Senatore, and M. Zaldarriaga, Cosmological Non-Linearities as an Effective Fluid, JCAP 07, 051, arXiv:1004.2488 [astro-ph.CO] .
  * Carrasco  _et al._ [2012] J. J. M. Carrasco, M. P. Hertzberg, and L. Senatore, The Effective Field Theory of Cosmological Large Scale Structures, JHEP 09, 082, arXiv:1206.2926 [astro-ph.CO] .
  * Chudaykin  _et al._ [2020] A. Chudaykin, M. M. Ivanov, O. H. E. Philcox, and M. Simonović, Nonlinear perturbation theory extension of the Boltzmann code CLASS, Phys. Rev. D 102, 063533 (2020), arXiv:2004.10607 [astro-ph.CO] .
  * Lanusse  _et al._ [2012] F. Lanusse, A. Rassat, and J. L. Starck, Spherical 3D isotropic wavelets, A&A 540, A92 (2012), arXiv:1112.0561 [astro-ph.CO] .
  * Eickenberg  _et al._ [2022] M. Eickenberg, E. Allys, A. Moradinezhad Dizgah, P. Lemos, E. Massara, M. Abidi, C. Hahn, S. Hassan, B. Regaldo-Saint Blancard, S. Ho, S. Mallat, J. Andén, and F. Villaescusa-Navarro, Wavelet Moments for Cosmological Parameter Estimation, arXiv e-prints , arXiv:2204.07646 (2022), arXiv:2204.07646 [astro-ph.CO] .
  * Foreman-Mackey  _et al._ [2013] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The MCMC Hammer, Publ. Astron. Soc. Pac. 125, 306 (2013), arXiv:1202.3665 [astro-ph.IM] .
  * Tassev  _et al._ [2015] S. Tassev, D. J. Eisenstein, B. D. Wandelt, and M. Zaldarriaga, sCOLA: The N-body COLA Method Extended to the Spatial Domain, (2015), arXiv:1502.07751 [astro-ph.CO] .
  * Leclercq  _et al._ [2020] F. Leclercq, B. Faure, G. Lavaux, B. D. Wandelt, A. H. Jaffe, A. F. Heavens, W. J. Percival, and C. Noûs, Perfectly parallel cosmological simulations using spatial comoving Lagrangian acceleration, Astron. Astrophys. 639, A91 (2020), arXiv:2003.04925 [astro-ph.CO] .
  * Modi  _et al._ [2021] C. Modi, F. Lanusse, and U. Seljak, FlowPM: Distributed TensorFlow implementation of the FastPM cosmological N-body solver, Astronomy and Computing, Volume 37, article id. 100505. 37, 100505 (2021).
  * Nadler  _et al._ [2023] E. O. Nadler, P. Mansfield, Y. Wang, X. Du, S. Adhikari, A. Banerjee, A. Benson, E. Darragh-Ford, Y.-Y. Mao, S. Wagner-Carena, R. H. Wechsler, and H.-Y. Wu, Symphony: Cosmological Zoom-in Simulation Suites over Four Decades of Host Halo Mass, ApJ 945, 159 (2023), arXiv:2209.02675 [astro-ph.CO] .
  * Li  _et al._ [2022] Y. Li, L. Lu, C. Modi, D. Jamieson, Y. Zhang, Y. Feng, W. Zhou, N. P. Kwan, F. Lanusse, and L. Greengard, pmwd: A Differentiable Cosmological Particle-Mesh $N$-body Library, (2022), arXiv:2211.09958 [astro-ph.IM] .
  * Goldstein  _et al._ [2022] S. Goldstein, A. Esposito, O. H. E. Philcox, L. Hui, J. C. Hill, R. Scoccimarro, and M. H. Abitbol, Squeezing fNL out of the matter bispectrum with consistency relations, Phys. Rev. D 106, 123525 (2022), arXiv:2209.06228 [astro-ph.CO] .
  * Cabass  _et al._ [2022a] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović, and M. Zaldarriaga, Constraints on Single-Field Inflation from the BOSS Galaxy Survey, Phys. Rev. Lett. 129, 021301 (2022a), arXiv:2201.07238 [astro-ph.CO] .
  * Cabass  _et al._ [2022b] G. Cabass, M. M. Ivanov, O. H. E. Philcox, M. Simonović, and M. Zaldarriaga, Constraints on multifield inflation from the BOSS galaxy survey, Phys. Rev. D 106, 043506 (2022b), arXiv:2204.01781 [astro-ph.CO] .
  * D’Amico  _et al._ [2022b] G. D’Amico, M. Lewandowski, L. Senatore, and P. Zhang, Limits on primordial non-Gaussianities from BOSS galaxy-clustering data, (2022b), arXiv:2201.11518 [astro-ph.CO] .


# Discovery of a long thermonuclear X-ray burst from the ultra-compact binary
4U 1850–087

Yongqi Lu  Key Laboratory of Stars and Interstellar Medium, Xiangtan
University, Xiangtan 411105, Hunan, P.R. China Zhaosheng Li, Yuanyue Pan
<EMAIL_ADDRESS><EMAIL_ADDRESS>Zhaosheng Li Key Laboratory of Stars
and Interstellar Medium, Xiangtan University, Xiangtan 411105, Hunan, P.R.
China  Wenhui Yu  Key Laboratory of Stars and Interstellar Medium, Xiangtan
University, Xiangtan 411105, Hunan, P.R. China  Yuanyue Pan  Key Laboratory of
Stars and Interstellar Medium, Xiangtan University, Xiangtan 411105, Hunan,
P.R. China  Maurizio Falanga International Space Science Institute (ISSI),
Hallerstrasse 6, 3012 Bern, Switzerland  Physikalisches Institut, University
of Bern, Sidlerstrasse 5, 3012 Bern, Switzerland

###### Abstract

We report the detection of a long X-ray burst triggered on MJD 60171.65 from
the ultra-compact binary 4U 1850–087 by the Monitor of All-sky X-ray Image and
Neutron Star Interior Composition Explorer (NICER). We analyse the NICER data
observed in between MJD 60095.19–60177.43, including one observation covered
part of the long X-ray burst tail, i.e., 0.15–3.8 hr after the trigger. The
persistent spectra are quite similar and well described by a combination of
multi-color disk blackbody, with the inner temperature of 0.5 keV, and a
thermally comptonized continuum with the asymptotic power-law photon index of
$\Gamma\sim 2.2$, and electron temperature of $kT_{\rm e}\sim 20-30$ keV. The
persistent fluxes were around $3.8\times 10^{-10}~{}{\rm
erg~{}cm^{-2}~{}s^{-1}}$, corresponding to a local accretion rate of
$1\%~{}\dot{m}_{\rm Edd}$. Part of time-resolved burst spectra show a clear
deviation from the blackbody model, which can be improved by considering the
enhanced persistent emission due to the Poynting-Robertson drag, or the
reflected disk emission illuminated by the burst. From the burst flux during
the cooling tail, we estimate the burst duration, $\tau\approx 0.78$ hr, the
burst fluence, $E_{\mathrm{b}}\approx 4.1\times 10^{41}$ ergs, and the
ignition column depth, $y_{\rm ign}\approx 3.5\times 10^{10}~{}{\rm
g~{}cm^{-2}}$. We propose that the long X-ray burst is powered by unstable
burning of pure helium in deep layer. Moreover, we identify significant 1 keV
emission lines in the burst spectra, which may originate from the surrounding
disk.

##  1 Introduction

In a neutron star low-mass X-ray binary (NS LMXB), NS accretes matter from its
companion star ($M<1M_{\odot}$) through the Roche–lobe overflow, behaving as a
persistent or transient source (Lasota, 2001). The accreted matter, mainly
composed of helium or mixed of hydrogen and helium, can be consumed through
stable, marginally stable or unstable burning, which dependents on the
accretion rate and the abundance (Keek et al., 2012; Li et al., 2021). The
unstable thermonuclear burning is also known as type I X-ray burst (hereafter
X-ray burst, see Strohmayer & Bildsten, 2006; Galloway & Keek, 2021). X-ray
bursts usually release total energy about 1039 erg and ignite at a typical
column depth of $y_{\rm ign}\sim 10^{8}\mathrm{~{}g~{}cm^{-2}}$. It lasts
$\sim$ 10–100 s, depending on the thickness of the burning layer (see Lewin et
al., 1993; Galloway et al., 2008, 2020). Intermediate-duration bursts can
release $\sim 10^{41}\mathrm{~{}erg}$ in a longer duration of $\sim$ 100–1000
s and are powered by the unstable burning of helium in deep layer at an
ignition column depth of $y_{\rm ign}\sim 10^{10}\mathrm{~{}g~{}cm^{-2}}$ (see
in’t Zand et al., 2005; Cumming et al., 2006; Falanga et al., 2008; Keek et
al., 2010). Most intermediate-duration bursts occurred in the low accretion
rate, that is lower than 1% of the Eddington limit (Falanga et al., 2008;
Alizai et al., 2023). Superbursts have the longest duration of $>10^{3}$ s and
releasing energies of $\sim 10^{42}\mathrm{~{}erg}$, which are powered by the
unstable burning of carbon for its mass fraction reaching 15–30% (Cumming &
Bildsten, 2001; Strohmayer & Brown, 2002). However, the X-ray bursts in 4U
0614+091 and SAX J1712.6–3739 with duration of hours were classified as
intermediate-duration bursts since the accretion rates were around 1% of the
local Eddington limit and cannot produce carbon sufficiently (see e.g.,
Kuulkers et al., 2010; in’t Zand et al., 2019; Alizai et al., 2023).

Several types of spectral features have been identified during X-ray bursts.
The 1 keV emission lines have been observed in IGR J17062–6143 (Degenaar et
al., 2013; Keek et al., 2017; Bult et al., 2021a) and SAX J1808.4–3658 (Bult
et al., 2021b), which were interpreted as a reflection feature from accretion
disk, possibly associated with Fe-L band transitions. The 1 keV emission line
was also observed in the bursts from 4U 1820–30 during the expansion phase
(Strohmayer et al., 2019). The reflection features were observed in normal or
long X-ray bursts (Strohmayer & Brown, 2002; Ballantyne & Strohmayer, 2004;
Keek et al., 2014; Zhao et al., 2022; Lu et al., 2023).

4U 1850–087 is a persistent X-ray source, which was first detected as an X-ray
burster by Swank et al. (1976) and immediately confirmed that the compact star
is an NS. The possible 20.6 min orbital period of 4U 1850–087 suggested it
belongs to an ultracompact X-ray binary (UCXB; Homer et al., 1996). 4U1850–087
is located in the globular cluster NGC 6712 at a distance of 8 kpc (Tremou et
al., 2018). The XMM-Newton spectra showed clear O K, Fe L, and Ne K edges at
0.54, 0.71, and 0.87 keV, respectively (Sidoli et al., 2005). Several long
X-ray bursts in 4U 1850–087 were detected by the Monitor of All-sky X-ray
Image (MAXI) and Swift which have been suggested as the power of unstable
burning of pure helium in deep layer due to their relatively high peak fluxes
and low accretion rates (e.g., in’t Zand et al., 2014; Serino et al., 2016).

In this work, we analyze a long X-ray burst detected by MAXI and NICER in 2023
August 15. In Section 2, we introduce the observations and describe the
properties of the persistent emissions and long X-ray burst. In Section 3, we
perform the spectra fitting for the persistent and time-resolved burst
emissions. We discuss the results in Section 4.

##  2 Observation and Data Reduction

MAXI and NICER are two instruments installed onboard the International Space
Station on 2009 July 15 and 2017 June 3, respectively. MAXI performs a full-
sky survey every 96 minutes, which can monitor the X-ray sources activities,
e.g., X-ray outbursts and X-ray bursts. From the MAXI novae
webpage111http://maxi.riken.jp/alert/novae/, we noticed a long X-ray burst
from 4U 1850–087 triggered on 2023 August 15 (MJD 60171.654514). We then
cross-checked the archived data from NICER, and found six contemporaneous
observations, including Obs. IDs 6617010101–6617010103 and
6403520101–6403520103 with the net unfiltered exposure time of 44 ks. The Obs.
ID 6403520101 started on MJD 60171.66085, 0.15 hrs after the trigger of MAXI,
partially covering the cooling tail of the long X-ray burst.

We used HEASOFT V6.31.1 and the NICER Data Analysis Software (NICERDAS) to
process the NICER data. The standard filtering criteria have been applied,
i.e., the pointing offset $<54^{\prime\prime}$, the Earth elevation angle
$>15^{\circ}$, the elevation angle with respect to the bright Earth limb $\leq
30^{\circ}$, and the instrument locating outside the South Atlantic Anomaly
(SAA). We used nicerl3-lc to extract the light curves in 0.5–10, 0.5–2, and
2–10 keV. The 64 s light curves in the energy range 0.5–10 keV are shown in
Figure 1. To have better coverage of the data between MJD 60171–60172, we
added the 2-20 keV light curves with 1 day binned from
MAXI/GSC222http://maxi.riken.jp/mxondem/ in Figure 1.

We show the 64 s light curves of 0.5–2, 2–10, and 0.5–10 keV from Obs. ID
6403520101 in Figure 2. The hardness ratio between 2–10 and 0.5–2 keV is
displayed in the bottom panel in Figure 2. The light curve contains three
nearly continuous segments, i.e., 0.15–0.25, 1.7–2.05, and 3.5–3.8 hr after
the MAXI trigger, which are marked as regions I, II, and III, respectively.
The count rate of Obs. ID 6403520101 decreased from 1250 to 120 cts/s in 3.6
hr. The count rate in region III is higher than the persistent count rate of
about 70 cts/s, therefore the three regions are all considered as part of the
burst. Meanwhile, the hardness ratio decreased from 1.3 to 1.1 in region I and
maintained at 0.8 and 0.4 in regions II and III, respectively, which indicates
the spectral soften during the burst cooling, see Sect. 3.2. We model the 1 s
burst light curve in the energy range 0.5–10 keV by an exponential function

$C(t)=C(t_{0})e^{-(t-t_{0})/\tau_{\rm LC}}+C_{0}.$ (1)

The model includes the normalization, $C(t_{0})$, the exponential decay time,
$\tau_{\rm LC}$, and a constant, $C_{0}$, which represents the asymptotic rate
of the burst. We fix the $t_{0}$ to the MAXI trigger time. The fitted model
gives an exponential decay time, $\tau=0.55\pm 0.01$ hr, with the
reduced-$\chi^{2}/{\rm dof}=1.35/2306$. We show the best fitted model in the
top panel of Figure 2.

We applied the Fast Fourier transform (FFT) method with Leahy normalization
(Leahy et al., 1983) to search for the burst oscillation in Obs. ID
6403520101. The power spectra were calculated in the 0.5–3, 3–10, and 0.5–10
keV energy ranges from the cleaned event file, respectively. We adopted a
moving window method with the window size $T$ of 8 s, and steps of 1 s. We
calculated the Fourier frequency between 50 and 2000 Hz with an interval of
0.25 Hz in each window (see, e.g. Bilous & Watts, 2019; Li et al., 2021). We
did not find any prominent signals during the burst in NICER observation, and
put an upper limit of the fractional root-mean-square amplitude about $5\%$.

Figure 1:  Light curves of 4U 1850–087 from MAXI and NICER observations. We
show the NICER light curve (64 s, 0.5–10 keV, black and red points) and the
MAXI light curve (one day averaged, 2–20 keV, blue points), where the red
points represent the NICER long X-ray burst’s data. We use the vertical grey
dashed line to represent the burst trigger time on MJD 60171.654514 from MAXI.
Figure 2: The burst light curves of 4U 1850–087 from NICER Obs. ID 6403520101.
From top to bottom, we show the 64 s NICER burst light curves in the energy
range 0.5–10 keV, 0.5–2 keV, and 2–10 keV, and the hardness ratio between
2.0–10 and 0.5–2 keV, respectively. In the top panel, we show the fitted
exponential model with the blue dashed line. The time intervals of the regions
I, II, and III are represented.

##  3 Spectral Analysis

We extracted the NICER spectra, ancillary response files (ARFs), and response
matrix files (RMFs) using nicerl3-spect, and the 3C50 background spectra were
also produced simultaneously (Remillard et al., 2021). We used XSPEC v12.13.0
to analyze the NICER spectra in 0.7–10 keV energy band (Arnaud, 1996). All
bolometric fluxes were estimated in the energy range 0.01-250 keV by using the
cflux model. The uncertainties were reported at 1 sigma confidence level.

###  3.1 Persistent Emissions

We fitted the persistent spectra from Obs. ID 6403520102-6403520103 and Obs.
ID 6617010101-6617010103 individually. All persistent spectra were grouped
with at least 50 counts per channel. The persistent spectra can be well
described by a combination of a disk blackbody plus a high-energy Comptonized
component (see e.g. Sidoli et al., 2006). We adopted the model
TBabs*(diskbb+nthcomp) to fit the spectra, where TBabs is the interstellar
absorption with abundances from Wilms et al. (2000); diskbb is a multi-disk
blackbody spectrum; nthcomp is a thermally Comptonized continuum (Zdziarski et
al., 1996; Życki et al., 1999). The model parameters include the asymptotic
power law index, $\Gamma$, the electron temperature, $kT_{\mathrm{e}}$, the
seed photon temperature, $kT_{\mathrm{bb,~{}seed}}$, the input type of seed
photons and the normalization for nthcomp; the inner disk temperature,
$kT_{\mathrm{in}}$, and normalization for diskbb; and the equivalent hydrogen
column, $N_{\rm H}$, for TBabs. We set the seed photons type from nthcomp to 1
and tied the $kT_{\mathrm{bb,~{}seed}}$ to $kT_{\mathrm{in}}$, which means the
input seed photons from the accretion disk. Considering the long exposure time
for few spectra, the spectra were fitted well from the reduced $\chi^{2}$,
$\chi^{2}_{\nu}\sim 0.9-1.3$. The best-fitted parameters are listed in Table
1.

The averaged $N_{\rm H}=0.606\pm{0.012}\times 10^{22}~{}{\rm cm}^{-2}$, is
consistent with the range of 0.4–0.6 $\times 10^{22}~{}{\rm cm}^{-2}$ from
Sidoli et al. (2006). The disk temperature, radius, and the power-law index of
the nthcomp component are $\sim$ 0.5 keV, 7 km, and 2.1, respectively. The
electron temperature is hard to constraints tightly and the uncertainties are
large. These persistent spectra showed quite similar profiles, implying that
the persistent emissions did not change significantly. Therefore, we adopted
the spectrum from Obs. ID 6403520102, which was close to the long X-ray burst,
to account for the persistent emission during the burst. We show the best-
fitted spectrum and residuals in Figure 3.

Figure 3: The un-absorbed persistent spectrum and residuals from NICER Obs. ID
6403520102 in 0.7–10 keV, which was regarded as the persistent emission during
the long X-ray burst. The red, yellow, and blue dashed dots present the model,
nthcomp component, and diskbb component, respectively.  Table 1:  Best fitted
parameters of all persistent spectra.

NICER | Exposure |  $N_{\rm H}$ | $\Gamma$ | $kT_{\rm e}^{\mathrm{a}}$ | $kT_{\rm in}$ | $R_{\rm in}^{\mathrm{b}}$ | $F_{\rm per}^{\mathrm{c}}$ |  $\chi_{\nu}^{2}(\mathrm{dof}$)
---|---|---|---|---|---|---|---|---
(Obs Id) | (s) | $(10^{21}\enspace\mathrm{cm^{-2}})$ |  | $\mathrm{(keV)}$ | $\mathrm{(keV)}$ | (km) | ($10^{-10}~{}{\rm erg~{}cm^{-2}~{}s^{-1}}$) |
6403520102 | 647 | 5.92$\pm{0.17}$ | $2.13^{+0.10}_{-0.18}$ | $18.46_{-17.30}^{+p}$ | $0.53\pm{0.04}$ | $6.76^{+1.32}_{-0.63}$ | $3.86\pm{0.03}$ | 0.89(294)
6403520103 | 1533 | 5.72$\pm{0.12}$ | $2.07^{+0.08}_{-0.09}$ | $28.50_{-27.50}^{+p}$ | $0.56\pm{0.03}$ | $6.08^{+0.38}_{-0.42}$ | $3.79\pm{0.02}$ | 1.05(403)
6617010101 | 8833 | 5.98$\pm{0.11}$ | $2.13^{+0.04}_{-0.04}$ | $31.31_{-30.05}^{+p}$ | $0.49\pm{0.02}$ | $7.57^{+0.27}_{-0.29}$ | $3.82\pm{0.02}$ | 1.29(643)
6617010102 | 4492 | 6.32$\pm{0.09}$ | $2.38^{+0.03}_{-0.08}$ | $18.31_{-17.25}^{+p}$ | $0.41\pm{0.03}$ | $7.14^{+1.03}_{-1.79}$ | $3.75\pm{0.02}$ | 1.24(454)
6617010103 | 3931(p) | 6.40$\pm{0.12}$ | $2.27^{+0.06}_{-0.08}$ | $14.70_{-10.67}^{+p}$ | $0.40\pm{0.04}$ | $8.71^{+1.21}_{-1.51}$ | $3.85\pm{0.02}$ | 1.13(438)

a The symbol $p$ means that the poorly constrained electron temperature is
pegged at the hard limit.

b We calculated the inner disk radius by assuming the angle of the disk,
$i=30^{\circ}$.

c The unabsorbed bolometric persistent flux in $0.1-250$ keV.

###  3.2 Time-resolved burst Spectroscopy

We extracted the time-resolved burst spectra with exposure time of 10 s in
region I and 100 s in regions II and III from the cleaned event file in Obs.
ID 6403520101. We grouped the burst spectra to guarantee that each channel has
at least 20 counts by using grappha. We first adopted the standard model,
TBabs $\times$ bbodyrad, to fit the burst spectra, where the parameters are
the temperature, $kT_{\rm bb}$, and normalization, $K=R_{\rm bb}/D_{\rm
10~{}kpc}$ where $D_{\rm 10~{}kpc}$ is the distance to the source in units of
10 kpc, for bbodyrad, and $N_{\rm H}$ for Tbabs. We fixed the $N_{\rm
H}=0.6\times 10^{22}~{}{\rm cm}^{-2}$ and took the persistent spectra from
Obs. ID 6403520102 as the burst background. The blackbody model can well
describe the burst spectra in regions II and III for $\chi_{\nu}^{2}$ close to
unity. However, for region I, we obtained the $\chi_{\nu}^{2}\sim 1.8-2.3$,
implying poor fits. We show the fitted parameter and $\chi_{\nu}^{2}$ as grey
stars in Figure 4. The improvements of the fit can be obtained by introducing
an extra component from the enhancement of persistent during burst, i.e., the
$f_{a}$-model (see Sect. 3.2.1), or a reflected emission from the accretion
disk, i.e., the reflection model (see Sect. 3.2.2).

####  3.2.1 Enhanced Persistent Emission

We used the $f_{\rm a}$-model, TBabs$\times$(bbodyrad + constant $\times$
(nthcomp + diskbb)), to fit the burst spectra, where constant accounts for the
enhancement of the persistent emission due to the Poynting-Robertson drag
(Walker, 1992; in’t Zand et al., 2013; Worpel et al., 2013, 2015). Here, we
assumed that only the amplitude of the persistent emission was changed if the
accretion rate increased during the burst. We adopted the 3c50 spectrum as the
instrumental background (Remillard et al., 2022). We fixed $N_{\rm H}$ to the
averaged value $0.6\times 10^{22}~{}{\rm cm}^{-2}$ obtained in Sect. 3.1, the
diskbb and nthcomp parameters to the best-fitted values from the persistent
spectra. We show the fitted parameter and $\chi_{\nu}^{2}$ as red points in
Figure 4. The $\chi_{\nu}^{2}$ of $f_{\rm a}$-model is around 1.0, which
implies the enhancement of the persistent emission can significantly improve
the fitted results. The blackbody temperature decreased from $\sim$ 1.5 keV
during region I to $\sim$ 0.55 keV during region III as the burst cooling. As
the temperature decreases, the blackbody flux decays from $1\times 10^{-8}$ to
$1.5\times 10^{-10}~{}{\rm erg~{}cm^{-2}~{}s^{-1}}$. The blackbody radius
keeps around 10 km during regions I, II, and III, implying the cooling of the
whole NS surface. The enhanced persistent emission was about five times the
pre-burst emission in region I and close to the pre-burst level in regions II
and III. We can not tell if this burst had experienced a photospheric radius
expansion from time-resolved spectroscopy due to the peak of the burst being
missed.

Figure 4: The time-resolved spectroscopy of the long X-ray burst observed by
NICER. All data points are belonged to part of the long X-ray burst. From top
to bottom, we show the bolometric flux of burst, $F_{\mathrm{burst}}$; the
blackbody temperature, $kT_{\mathrm{bb}}$, the blackbody radius, $R_{\rm bb}$,
which were calculated at a distance of 8 kpc; the enhanced persistent flux,
$F_{\mathrm{f_{a}}}$, or reflection flux, $F_{\mathrm{refl}}$; and the
goodness of fit per degree of freedom, $\chi_{\nu}^{2}$. The grey star, red
dot, and blue square display the results from the blackbody model, the $f_{a}$
model, and the reflection model, respectively. The blue dashed line represent
the persistent emission level. The time zero corresponds to about 0.15 hr
after the long burst starts.

####  3.2.2 The reflection model

The residual of the absorbed blackbody model can also be well described by
adding a reflection emission which is illuminated by the long X-ray burst. We
use the model relxillNS (García et al., 2022) to represent the reflection
component (see also Zhao et al., 2022; Lu et al., 2023; Yu et al., 2023). We
fitted the burst spectra by TBabs
$\times(\texttt{bbodyrad}+\texttt{relxillNS})$. We fixed all parameters to a
constant except the normalization. We tied the input black body temperature,
$kT_{\mathrm{bb,~{}in}}$, to the temperature of bbodyrad, and set the
reflection fraction parameter to -1, which means that this model only
represents the reflection emission, and the accretion disk is illuminated by
bbodyrad photon. We fixed the indices of emissivity, $q_{1}=q_{2}=3$, the
inner radius of the accretion disk, $R_{\mathrm{in}}=6R_{g}$, and the outer
radius of the accretion disk, $R_{\mathrm{out}}=400R_{g}$, where the $R_{g}$
is the gravitational radius. Since the spin of 4U 1850–087 is unknown, we
assume the spin parameter, $a$, of the model relxillNS to 0.1. We also tried
two typical inclination angles, $30^{\circ}$ and $60^{\circ}$, to fit the
spectra, and found smaller inclination provides slightly better
$\chi_{\nu}^{2}$. Therefore, we use $30^{\circ}$ as the inclination angle. We
let the ionization parameter, $\log\xi$, the iron abundance in the accretion
disk in units of solar abundance, $A_{\mathrm{Fe}}$, and the logarithmic
density of the accretion disk, $\log n$, free to vary. We found the mean
values of $\log\xi$, $A_{\mathrm{Fe}}$, and $\log n$ are around 2.8, 6, and 16
cm-3 during the burst, respectively. Therefore, we fixed these parameters to
the corresponding values. We only show the parameters of relxillNS in region
I, due to the reflection component being insignificant during the cooling tail
in regions I and II. The $\chi_{\nu}^{2}$ are close to 1 in region I which
means the reflection model can also explain the derivations. The best-fitting
results are shown in Figure 4 as blue squares. Compared with the $f_{\rm
a}$-model, the reflection model provides 5% lower blackbody temperatures, 5%
higher blackbody radii, and thus 15% lower bolometric fluxes.

###  3.3 The 1 keV emission line

During the time-resolved burst spectra analysis, we found a faint emission
line around 1 keV in region I with an exposure time of 10 s. In order to
improve the significance of the emission line, we extracted six spectra in
region I with an exposure time of around 50 s for each of them. We used the
$f_{\rm a}$-model to fit these spectra, same as section 3.2.1. We also
extracted longer-duration spectra in regions II and III to find this possible
1 keV feature, but no significant features were found. We also did not find
any prominent feature around 1 keV in the persistent spectra. We show the
residuals of six spectra in region I, and also the residuals of spectra
extracted from regions II and III as comparisons in Figure 5.

We added a Gaussian to account for the emission line. The Gaussian parameters
of these six spectra are shown in Table 2. We used the simftest script in
Xspec to evaluate the statistical significance of the Gaussian component. For
each spectrum, the estimation contains $10^{5}$ sets of simulated data sets.
We obtained the probability of the emission line arising by chance less than
$1.5\times 10^{-14}$, which means all lines’ significance are higher than
$7\sigma$. We list the centroid energy in Table 2. By using the reflection
model, we also observed this 1 keV feature, but slightly smaller significance.

Figure 5: The burst residuals obtained by fitting the $f_{\rm a}$-model. We show the error-weighted residuals from the region I-III. The 1 keV emission features are clearly shown in the spectra from region I but not in the spectra from regions II and III. Table 2: Gaussian line parameters.  Spectra in region I | $E_{\rm center}$ | Line Norm
---|---|---
| (keV) | $(\mathrm{photons~{}cm^{-2}~{}s^{-1}})$
1 | 0.970$\pm{0.015}$ | 0.034$\pm{0.005}$
2 | 0.984$\pm{0.012}$ | 0.037$\pm{0.005}$
3 | 0.991$\pm{0.012}$ | 0.032$\pm{0.005}$
4 | 0.996$\pm{0.013}$ | 0.030$\pm{0.005}$
5 | 0.983$\pm{0.014}$ | 0.035$\pm{0.005}$
6 | 0.984$\pm{0.012}$ | 0.037$\pm{0.005}$

##  4 Discussion and conclusion

###  4.1 Burst and persistent spectra

In this work, we report the long X-ray burst from 4U 1850–087 observed by
NICER and MAXI observations in 2023. The persistent spectra can be well
described by a combination of a multi-color disk blackbody with the inner
temperature around 0.5 keV, and a thermally comptonized continuum with the
asymptotic power-law photon index of $\Gamma\sim 2.2$, and electron
temperature $kT_{\rm e}\sim 20-30$ keV. The persistent spectra are similar
across all NICER observations, therefore we adopted a spectrum after the long
X-ray burst burst to represent the persistent emission during burst. We
extracted the light curve and identified a long X-ray burst from the Obs. ID
6403520101. To describe the burst evolution, we divide the observation into
regions I, II, and III according to the light curve. The spectra from regions
II and III were well described by a single blackbody model, but the region I
spectra showed clear derivations from the blackbody shape. The results can be
improved by adding the enhanced persistent emission due to the
Poynting–Robertson drag, the $f_{a}$-model, or the reflection emission
illuminated by the burst, the reflection model. These two models provided
similar burst parameters and $\chi_{\nu}^{2}$, and are difficult to
distinguish.

We use

$\begin{split}\dot{m}&=\frac{L_{\mathrm{{per}}}(1+z)}{4\pi R_{\rm
NS}^{2}(GM_{\rm NS}/R_{\rm NS})}\\\ &\approx 1.7\times
10^{3}\biggl{(}\frac{F_{\mathrm{{per}}}}{4\times
10^{-10}\mathrm{~{}ergs~{}cm^{-2}~{}s^{-1}}}\biggr{)}\biggl{(}\frac{d}{8\mathrm{~{}kpc}}\biggr{)}^{2}\\\
&\quad\times\ \biggl{(}\frac{M_{\rm
NS}}{1.4M_{\odot}}\biggr{)}^{-1}\biggl{(}\frac{1+z}{1.31}\biggr{)}\biggl{(}\frac{R\mathrm{{}_{NS}}}{10\mathrm{~{}km}}\biggr{)}^{-1}\mathrm{g~{}cm^{-2}}\mathrm{~{}s^{-1}},\end{split}$
(2)

to calculate the local accretion rate, where the $F_{\rm per}$ is the
persistent flux, the mass and radius of NS are $M_{\mathrm{NS}}=1.4M_{\odot}$
and $R_{\mathrm{NS}}=10$ km, the gravitational redshift on NS surface,
$z\approx 0.31$, the source distance, $d=8$ kpc. We obtained
$\dot{m}=(1.63\pm{0.01})\times 10^{3}~{}\mathrm{g~{}cm^{-2}~{}s^{-1}}$ for a
distance $8$ kpc (Paltrinieri et al., 2001). The local Eddington accretion
rate is $\dot{m}_{\rm Edd}=(8.8\times
10^{4})[1.7/(X+1)]~{}\mathrm{g~{}cm^{-2}~{}s^{-1}}$, where $X$ represents the
fraction of hydrogen in the accreted material. Considered 4U 1850-087 as an
UCXB, $\dot{m}\sim 1\%~{}\dot{m}_{\rm Edd}$.

###  4.2 Burst parameters

The decay of the long X-ray burst flux can be well described with the analytic
expression by Cumming & Macbeth (2004) and Cumming et al. (2006), leading to
constraints on the energy release per unit mass, $E_{17}$ in units of
$10^{17}\mathrm{~{}ergs~{}g^{-1}}$, and the ignition column depth $y_{12}$ in
units of $10^{12}\mathrm{~{}g~{}cm^{\mathbf{-2}}}$. We used the trigger time
from MAXI as the burst peak time which is about 0.15 hr earlier than the first
NICER observation after the burst. We obtained a poor
reduced-$\chi^{2}$/dof=2.94/47 by using all three regions. The burst fluxes in
Region III contribute the most significant residuals, which have been also
shown in other long bursts (see e.g., Cumming & Macbeth, 2004; Cumming et al.,
2006). Therefore we only fit the fluxes from Regions I and II. The fit is
improved to an acceptable reduced-$\chi^{2}$/dof=1.02/37. We show the fitted
$f_{\rm a}$-model bolometric flux in Figure 6. We estimated
$E_{17}=$0.718$\pm{0.001}$, $y_{12}=$0.035$\pm{0.001}$ and the burst peak flux
$F_{\mathrm{peak}}=$1.88$\pm{0.07}\times
10^{-8}\mathrm{~{}erg~{}cm^{-2}~{}s^{-1}}$. The burst peak flux is less than
the Eddington flux, that is $F_{\mathrm{Edd}}\approx 3.8\times
10^{-8}~{}\mathrm{erg~{}cm^{-2}~{}s^{-1}}$ in hydrogen poor environment, which
indicates the photospheric radius expansion may not appear in the burst.

We then use

$f_{\mathrm{{b}}}=\frac{4\pi y_{\rm ign}R_{\rm
NS}^{2}Q\mathrm{{}_{nuc}}}{{4\pi d^{2}(1+z)}},$ (3)

to calculate the burst fluence from the estimated ignition column depth
$y_{\rm ign}$, where $Q\mathrm{{}_{nuc}}\approx
1.31\mathrm{~{}MeV~{}nucleon^{-1}}$ for a given hydrogen fraction of $X=0$
(Goodwin et al., 2019). We obtained the burst fluence
$f_{\mathrm{{b}}}=$5.3$\pm{0.2}\times 10^{-5}\mathrm{~{}erg~{}cm^{-2}}$. We
then calculated the burst decay time $\tau=f_{\rm b}/F_{\rm peak}\approx
0.78\pm{0.04}$ hr, which is higher than the exponential decay time 0.55 hr
fitted from the light curve. All the burst parameters are reported in Table 3.
We added two intermediate-duration bursts’ parameters of 4U 1850-087 from MAXI
and Swift data in Table 3 as a comparison (Serino et al., 2016; in’t Zand et
al., 2014). The burst fluence and total released energy are rescaled to a
distance of 8 kpc. Compared to the other two bursts, the burst has a smaller
peak flux, a similar total burst emitted energy, and a longest burst decay
time. The predicted recurrence time is obtained via the relation $\Delta
t_{\rm rec}=y(1+z)/\dot{m}$. We found that the observed recurrence time,
$\Delta T_{\rm rec}$, of the long X-ray burst observed in 2015, which is
higher than the predicted value. The discrepancy may be caused by the not well
calibrated flux from MAXI, i.e., the burst peak flux is two times higher than
the Eddington limit. The predicted recurrence time of the long X-ray in our
work, 0.87 yr, is lower than the observed value, 7.76 yr, which means that
some long or normal X-ray bursts were likely missed during 2015-2023 due to
poorly coverage of the observations.

Figure 6: The burst flux fitted by the model which proposed in Cumming &
Macbeth (2004) during the burst decay. The burst flux was estimated from the
$f_{a}$ model. The black points represent the data in regions I and II. The
red points represent the region III data, which were excluded from the fit. We
use the trigger time in MAXI as the burst peak time. Table 3: The parameters
of long X-ray bursts from 4U 1850–087.

Burst date | $F_{\rm peak}$ | $E_{\rm b}$ | $f_{\rm b}$ | $\tau$ | $\Delta t_{\mathrm{rec}}\equiv y(1+z)/\dot{m}$ | $\Delta T_{\mathrm{rec}}^{\mathrm{a}}$
---|---|---|---|---|---|---
| $(10^{-8}~{}\mathrm{erg~{}cm^{-2}~{}s^{-1}})$ | $(10^{41}~{}\mathrm{erg})$ | $(10^{-5}~{}\mathrm{erg~{}cm^{-2}})$ | (hr) | (year) | (year)
2014-03-10b | 10.7$\pm{0.9}$ | $7.9\pm{0.7}~{}$ | $10.3\pm{0.9}$ | 0.27 | 1.70 | -
2015-05-09b | 7.8$\pm{0.5}$ | $14.8\pm{0.9}~{}$ | $19.3\pm{1.2}$ | 0.71 | 3.18 | 1.16
2015-11-11c | - | - | - | - | - | 0.51
2023-08-15 | 1.88$\pm{0.07}$ | 4.1$\pm{0.2}$ | 5.3$\pm{0.2}$ | 0.78$\pm{0.04}$ | 0.87$\pm{0.03}$ | 7.76

${}^{\mathrm{a}}~{}T_{\mathrm{rec}}$ means the time interval between the
current burst and the previous one. $T_{\mathrm{rec}}$ only presents the upper
limit of burst recurrence time due to the data gap of observations.

b The data of bursts in 2014 and 2015 were taken simply from in’t Zand et al.
(2014); Serino et al. (2016). Note that the distance of 8 kpc rather than 6.9
kpc was adopted to recalculate $E_{\rm b}$ and $f_{\rm b}$.

c The trigger time of bursts were taken from the MAXI novae webpage.

###  4.3 Burst fuel

The decay time of the burst from 4U 1850–087 on August 15, 2023, $\tau\sim
0.78$ h, which is within the longest intermediate-duration bursts and shortest
superbursts. The burst decay time in some longest intermediate-duration bursts
are 0.8 hr from SAX J1712.6–3739 (in’t Zand et al., 2019), and 0.71 hr from 4U
1850-087 (see Table 3, Serino et al., 2016). The burst duration of several
shortest superbursts were 0.7 and 1.0 hr from GX 17+2, 1.2 hr from Ser X–1,
and 2.1 hr from 4U 0614+091 (Keek et al., 2010). The accretion rate of GX 17+2
and Ser X–1 are $>20\%~{}\dot{m}_{\rm Edd}$, of 4U 0614+091 are
$<1\%~{}\dot{m}_{\rm Edd}$. The donor star of 4U 1850–08 is a white dwarf or
helium star, which means the accretion matter includes a large pile of helium.
At a low accretion rate $\sim 1\%~{}\dot{m}_{\rm Edd}$, the carbon is
difficult to accumulate (Cumming et al., 2006). Considering the burst decay
time is less than the superburst with a low accretion rate from 4U 0614+091
and the calculated recurrence time consistent with the observed recurrence
time between the other two intermediate-duration bursts from 4U 1850–087 (see
Table 3), we propose that this burst is the longest intermediate-duration
burst like other long X-ray bursts observed in 4U 1850–08 (in’t Zand et al.,
2014; Serino et al., 2016). We also calculated the ratio of the total
persistent flux between two bursts to the total energy emitted from burst
$\alpha=(F_{\mathrm{pers}}\times\Delta T_{\mathrm{rec}})/f_{\mathrm{b}}\sim
117$, which the shortest observed recurrence time of 0.51 hr in Table 3 is
used. When we assume all accreted fuel burned during the burst, we calculated
the theoretical
$\alpha=44~{}M_{\mathrm{NS}}R_{\mathrm{NS}}^{-1}(Q_{\mathrm{nuc}}/4.4\mathrm{MeV/nucleon})^{-1}\approx
121$ (Falanga et al., 2008). The calculated $\alpha$ from measurable
quantities is close to a pure helium burst, which is also support the burst as
an intermediate-duration burst.

###  4.4 The origin of the 1 keV emission line

We found an emission feature near 1 keV in the burst spectrum, which can be
improved by adding a Gaussian component. This line can only be detected in
region I, but not in regions II and III, and the persistent spectra, which
suggests that this line does not come from the instrumental effect. By using
the tool simfest, we estimated the significance of the emission line,
$>7\sigma$, which proves that it is not caused by random fluctuations.
Therefore, we propose that this line is caused by an astrophysical origin.

A weaker 1 keV emission line was also detected in an intermediate duration
burst of 4U 1850–087 by Swift (see in’t Zand, 2017). The similar emission line
have also been detected in IGR J17062–6143 (Bult et al., 2021a), SAX
J1808.4–3658 (Bult et al., 2019), and 4U 1820–30 (Strohmayer et al., 2019).
The line in IGR J17062–6143 has been explained by the Fe L shell transitions
or Ly$\alpha$ transition of Ne X at 1.022 keV in the accretion disk (Keek et
al., 2017; Degenaar et al., 2013) or the emit from photoionized plasma
illuminated by a blackbody continuum (Bult et al., 2021a). The 1 keV emission
line in 4U 1820–30 was detected in an PRE burst, which were possibly produced
in the PRE wind (Strohmayer et al., 2019). Considering we did not find the
photosphere expansion in time-resolved spectra, the 1 keV emission line is
likely produced in the accretion disk.

This work was supported by the Major Science and Technology Program of
Xinjiang Uygur Autonomous Region (No. 2022A03013-3). Z.S.L. and Y.Y.P. were
supported by National Natural Science Foundation of China (12273030,
12103042). This work made use of data from the High Energy Astrophysics
Science Archive Research Center (HEASARC), provided by NASA’s Goddard Space
Flight Center.

## References

  * Alizai et al. (2023) Alizai, K., Chenevez, J., Cumming, A., et al. 2023, MNRAS, 521, 3608, doi: 10.1093/mnras/stad374
  * 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
  * Ballantyne & Strohmayer (2004) Ballantyne, D. R., & Strohmayer, T. E. 2004, ApJ, 602, L105, doi: 10.1086/382703
  * Bilous & Watts (2019) Bilous, A. V., & Watts, A. L. 2019, ApJS, 245, 19, doi: 10.3847/1538-4365/ab2fe1
  * Bult et al. (2019) Bult, P., Jaisawal, G. K., Güver, T., et al. 2019, ApJ, 885, L1, doi: 10.3847/2041-8213/ab4ae1
  * Bult et al. (2021a) Bult, P., Altamirano, D., Arzoumanian, Z., et al. 2021a, ApJ, 920, 59, doi: 10.3847/1538-4357/ac18c4
  * Bult et al. (2021b) —. 2021b, ApJ, 907, 79, doi: 10.3847/1538-4357/abd54b
  * Cumming & Bildsten (2001) Cumming, A., & Bildsten, L. 2001, ApJ, 559, L127, doi: 10.1086/323937
  * Cumming & Macbeth (2004) Cumming, A., & Macbeth, J. 2004, ApJ, 603, L37, doi: 10.1086/382873
  * Cumming et al. (2006) Cumming, A., Macbeth, J., in ’t Zand, J. J. M., & Page, D. 2006, ApJ, 646, 429, doi: 10.1086/504698
  * Degenaar et al. (2013) Degenaar, N., Miller, J. M., Wijnands, R., Altamirano, D., & Fabian, A. C. 2013, The Astrophysical Journal, 767, L37, doi: 10.1088/2041-8205/767/2/l37
  * Falanga et al. (2008) Falanga, M., Chenevez, J., Cumming, A., et al. 2008, A&A, 484, 43, doi: 10.1051/0004-6361:20078982
  * Galloway & Keek (2021) Galloway, D. K., & Keek, L. 2021, Astrophysics and Space Science Library, 461, 209, doi: 10.1007/978-3-662-62110-3_5
  * Galloway et al. (2008) Galloway, D. K., Muno, M. P., Hartman, J. M., Psaltis, D., & Chakrabarty, D. 2008, The Astrophysical Journal Supplement Series, 179, 360
  * Galloway et al. (2020) Galloway, D. K., in’t Zand, J., Chenevez, J., et al. 2020, ApJS, 249, 32, doi: 10.3847/1538-4365/ab9f2e
  * García et al. (2022) García, J. A., Dauser, T., Ludlam, R., et al. 2022, ApJ, 926, 13, doi: 10.3847/1538-4357/ac3cb7
  * Goodwin et al. (2019) Goodwin, A. J., Heger, A., & Galloway, D. K. 2019, ApJ, 870, 64, doi: 10.3847/1538-4357/aaeed2
  * Homer et al. (1996) Homer, L., Charles, P. A., Naylor, T., et al. 1996, MNRAS, 282, L37, doi: 10.1093/mnras/282.3.L37
  * in’t Zand (2017) in’t Zand, J. 2017, in 7 years of MAXI: monitoring X-ray Transients, ed. M. Serino, M. Shidatsu, W. Iwakiri, & T. Mihara, 121.  https://arxiv.org/abs/1702.04899
  * in’t Zand et al. (2014) in’t Zand, J., Linares, M., & Markwardt, C. 2014, The Astronomer’s Telegram, 5972, 1
  * in’t Zand et al. (2005) in’t Zand, J. J. M., Cumming, A., van der Sluys, M. V., Verbunt, F., & Pols, O. R. 2005, A&A, 441, 675, doi: 10.1051/0004-6361:20053002
  * in’t Zand et al. (2019) in’t Zand, J. J. M., Kries, M. J. W., Palmer, D. M., & Degenaar, N. 2019, A&A, 621, A53, doi: 10.1051/0004-6361/201834270
  * in’t Zand et al. (2013) in’t Zand, J. J. M., Galloway, D. K., Marshall, H. L., et al. 2013, A&A, 553, A83, doi: 10.1051/0004-6361/201321056
  * Keek et al. (2014) Keek, L., Ballantyne, D. R., Kuulkers, E., & Strohmayer, T. E. 2014, ApJ, 797, L23, doi: 10.1088/2041-8205/797/2/L23
  * Keek et al. (2010) Keek, L., Galloway, D. K., in’t Zand, J. J. M., & Heger, A. 2010, ApJ, 718, 292, doi: 10.1088/0004-637X/718/1/292
  * Keek et al. (2012) Keek, L., Heger, A., & in’t Zand, J. J. M. 2012, ApJ, 752, 150, doi: 10.1088/0004-637X/752/2/150
  * Keek et al. (2017) Keek, L., Iwakiri, W., Serino, M., et al. 2017, ApJ, 836, 111, doi: 10.3847/1538-4357/836/1/111
  * Kuulkers et al. (2010) Kuulkers, E., in’t Zand, J. J. M., Atteia, J. L., et al. 2010, A&A, 514, A65, doi: 10.1051/0004-6361/200913210
  * Lasota (2001) Lasota, J.-P. 2001, New A Rev., 45, 449, doi: 10.1016/S1387-6473(01)00112-9
  * Leahy et al. (1983) Leahy, D. A., Darbro, W., Elsner, R. F., et al. 1983, ApJ, 266, 160, doi: 10.1086/160766
  * Lewin et al. (1993) Lewin, W. H. G., van Paradijs, J., & Taam, R. E. 1993, Space Sci. Rev., 62, 223, doi: 10.1007/BF00196124
  * Li et al. (2021) Li, Z., Pan, Y., & Falanga, M. 2021, ApJ, 920, 35, doi: 10.3847/1538-4357/ac1f15
  * Lu et al. (2023) Lu, Y., Li, Z., Pan, Y., et al. 2023, A&A, 670, A87, doi: 10.1051/0004-6361/202244984
  * Paltrinieri et al. (2001) Paltrinieri, B., Ferraro, F. R., Paresce, F., & De Marchi, G. 2001, AJ, 121, 3114, doi: 10.1086/321069
  * Remillard et al. (2021) Remillard, R. A., Loewenstein, M., Steiner, J. F., et al. 2021, arXiv e-prints, arXiv:2105.09901.  https://arxiv.org/abs/2105.09901
  * Remillard et al. (2022) —. 2022, AJ, 163, 130, doi: 10.3847/1538-3881/ac4ae6
  * Serino et al. (2016) Serino, M., Iwakiri, W., Tamagawa, T., et al. 2016, PASJ, 68, 95, doi: 10.1093/pasj/psw086
  * Sidoli et al. (2005) Sidoli, L., La Palombara, N., Oosterbroek, T., & Parmar, A. N. 2005, A&A, 443, 223, doi: 10.1051/0004-6361:20053438
  * Sidoli et al. (2006) Sidoli, L., Paizis, A., Bazzano, A., & Mereghetti, S. 2006, A&A, 460, 229, doi: 10.1051/0004-6361:20065759
  * Strohmayer & Bildsten (2006) Strohmayer, T., & Bildsten, L. 2006, in Compact stellar X-ray sources, Cambridge Astrophysics Series, No. 39, ed. W. Lewin & M. van der Klis (Cambridge: Cambridge University Press), 113–156
  * Strohmayer & Brown (2002) Strohmayer, T. E., & Brown, E. F. 2002, ApJ, 566, 1045, doi: 10.1086/338337
  * Strohmayer et al. (2019) Strohmayer, T. E., Altamirano, D., Arzoumanian, Z., et al. 2019, ApJ, 878, L27, doi: 10.3847/2041-8213/ab25eb
  * Swank et al. (1976) Swank, J. H., Becker, R. H., Pravdo, S. H., Saba, J. R., & Serlemitsos, P. J. 1976, IAU Circ., 3010, 1
  * Tremou et al. (2018) Tremou, E., Strader, J., Chomiuk, L., et al. 2018, ApJ, 862, 16, doi: 10.3847/1538-4357/aac9b9
  * Walker (1992) Walker, M. A. 1992, ApJ, 385, 642, doi: 10.1086/170969
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, The Astrophysical Journal, 542, 914, doi: 10.1086/317016
  * Worpel et al. (2013) Worpel, H., Galloway, D. K., & Price, D. J. 2013, ApJ, 772, 94, doi: 10.1088/0004-637X/772/2/94
  * Worpel et al. (2015) —. 2015, ApJ, 801, 60, doi: 10.1088/0004-637X/801/1/60
  * Yu et al. (2023) Yu, W., Li, Z., Lu, Y., et al. 2023, arXiv e-prints, arXiv:2312.16420, doi: 10.48550/arXiv.2312.16420
  * Zdziarski et al. (1996) Zdziarski, A. A., Johnson, W. N., & Magdziarz, P. 1996, MNRAS, 283, 193, doi: 10.1093/mnras/283.1.193
  * Zhao et al. (2022) Zhao, G., Li, Z., Pan, Y., et al. 2022, A&A, 660, A31, doi: 10.1051/0004-6361/202142801
  * Życki et al. (1999) Życki, P. T., Done, C., & Smith, D. A. 1999, MNRAS, 309, 561, doi: 10.1046/j.1365-8711.1999.02885.x