<figure><img src="image_1.png" /></figure>

# NuSTAR Observations of the Transient Galactic Black Hole Binary Candidate Swift J1858.6-0814: A New Sibling of V404 Cyg and V4641 Sgr?

Jeremy Hare<sup>1,2,8</sup>①, John A. Tomsick<sup>1</sup>①, Douglas J. K. Buisson<sup>3,4</sup>, Maïca Clavel<sup>5</sup>①, Poshak Gandhi<sup>4</sup>①, Javier A. García<sup>6,7</sup>①, Brian W. Grefenstette<sup>6</sup>①, Dominic J. Walton<sup>3</sup>①, and Yanjun Xu<sup>6</sup>①

## Abstract

Swift J1858.6-0814 was discovered by the Burst Alert Telescope on Swift on 2018 October 25. Here we report on the first follow-up Nuclear Spectroscopic Telescope Array (NuSTAR) observation of the source, which shows variability spanning two orders of magnitude in count rate on timescales of ~10–100 s. The power spectrum of the source does not show any quasi-periodic oscillations or periodicity, but has a large fractional rms amplitude of 147% ± 3%, exhibiting a number of large flares throughout the observation. The hardness ratio (defined as $R_{10–79\text{keV}}/R_{3–10\text{keV}}$) of the flares tends to be soft, while the source spans a range of hardness ratios during nonflaring periods. The X-ray spectrum of the source shows strong reflection features, which become more narrow and peaked during the nonflaring intervals. We fit an absorbed relativistic reflection model to the source spectra to place physical constraints on the system. Most notably, we find that the source exhibits a large and varying intrinsic absorbing column density ($N_H = (1.4–4.2) \times 10^{23} \text{ cm}^{-2}$). This large intrinsic absorption is further supported by the energy spectra extracted from two flares observed simultaneously by NuSTAR and the Neutron Star Interior Composition Explorer. We find that the inner accretion disk of the source has a low inclination, $i < 29^\circ$ (3$\sigma$ upper limit), while the iron abundance in the disk is close to solar, $A_{\text{Fe}} = 1.0 \pm 0.3$. We set a 90% confidence upper limit on the inner radius of the accretion disk of $r_{\text{in}} < 8r_{\text{ISCO}}$ and, by fixing $r_{\text{in}}$ to be at $r_{\text{ISCO}}$, a 90% confidence lower limit on the spin of the black hole of $a^* > 0.0$. Finally, we compare the properties of Swift J1858.6-0814 to those of V404 Cygni and V4641 Sgr, which both show rapid flaring and a strong and variable absorption.

*Unified Astronomy Thesaurus concepts: Stellar mass black holes (1611); X-ray transient sources (1852); Low-mass X-ray binary stars (939)*

## 1. Introduction

Since its launch in 2012, the Nuclear Spectroscopic Telescope Array (NuSTAR; Harrison et al. 2013) has played a pivotal role in studying, and sometimes identifying, the nature of Galactic hard X-ray transients discovered by, e.g., Swift's Burst Alert Telescope (BAT), the International Gamma-Ray Astrophysics Laboratory, or Monitor of All-sky X-ray Image. Among some of the most interesting sources include a new magnetar (Mori et al. 2013), several supergiant fast X-ray transients (e.g., Bhalerao et al. 2015; Ferrigno et al. 2019; Hare et al. 2019), and numerous low-mass X-ray binaries hosting either a neutron star (e.g., Homan et al. 2018; Jaisawal et al. 2018) or a black hole (BH; e.g., Xu et al. 2017; Beri et al. 2019). Following the detection and identification of a new BH transient candidate, NuSTAR has also helped to constrain the physical parameters of these systems through spectral fitting, using a combination of reflection models and NuSTAR's unprecedented sensitivity above 10 keV (see, e.g., Xu et al. 2018a, 2018b; Buisson et al. 2019).

As a BH transient undergoes an outburst it typically evolves through several spectral states, showing relatively slow variability on kilosecond- to day-long timescales (see, e.g., Remillard & McClintock 2006; Belloni & Motta 2016). These

outbursts usually start in the hard state, in which the X-ray spectrum is dominated by a hard power-law component. The source then transitions into the soft state, where the X-ray spectrum becomes dominated by the hot thermal emission from the accretion disk. Finally, the BH transient returns back to the hard state at the end of the outburst. While in the soft state, the accretion flow is expected to reach the innermost stable circular orbit (ISCO) of the BH, whose radius ($r_{\text{ISCO}}$) depends on the spin of the BH. Modeling the X-ray spectra during these spectral states with relativistic reflection models allows for constraints to be placed on the BH's spin.

While the majority of BH transients generally follow the standard progression through the spectral states outlined in the previous paragraph, there are a few outliers, such as V404 Cyg and V4641 Sgr. These systems exhibit large amplitude flares, with X-ray count rates rising by factors of $10^2–10^3$ on timescales of seconds to minutes, and reaching Eddington or even super-Eddington luminosities (Wijnands & van der Klis 2000; Revnivtsev et al. 2002a; Gandhi et al. 2017; Motta et al. 2017b; Walton et al. 2017). Accompanying these flares are significant changes in the shape of the X-ray spectra, including variations of the intrinsic absorbing column density, photon index, and reflection strength (see, e.g., Wijnands & van der Klis 2000; Motta et al. 2017b; Walton et al. 2017). The large-scale flaring behavior and rapid changes in the X-ray

<footnote n="1">¹Space Sciences Laboratory, 7 Gauss Way, University of California, Berkeley, CA 94720, USA; jeremy.hare@nasa.gov</footnote>

<footnote n="2">²NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="3">³Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="4">⁴Department of Physics and Astronomy, University of Southampton, Highfield, Southampton SO17 1BJ, UK</footnote>

<footnote n="5">⁵Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="6">⁶Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="7">⁷Dr. Karl Remeis-Observatory and Erlangen Centre for Astroparticle Physics, Sternwartstr. 7, D-96049 Bamberg, Germany
*Received 2019 November 4; revised 2019 December 27; accepted 2020 January 8; published 2020 February 12*</footnote>

<footnote n="8">⁸NASA Postdoctoral Program Fellow.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Swift-BAT light curve of Swift J1858.6–0814 (with 10 day averaged time bins) up to 2019 October 28 (to date the outburst is still ongoing). The thin vertical green and blue lines show times of the six NuSTAR observations of the source. The green line marks the observation reported here. The wide pink band shows the time in which the source was Sun constrained for NuSTAR.</figcaption></figure>

spectrum of these sources makes it difficult to characterize their spectral state; however, relativistic reflection modeling can still place constraints on their physical parameters (see, e.g., Walton et al. 2017).

Swift J1858.6–0814 (hereafter J1858), discovered as a Galactic ($l = 26°, 395, b = -5°, 351$) X-ray transient by Swift-BAT on 2018 October 25, is a new BH candidate exhibiting similar characteristics to V404 Cyg and V4641 Sgr (Krimm et al. 2018; Ludlam et al. 2018). The source was subsequently followed up by NuSTAR and the Neutron Star Interior Composition Explorer (NICER). The NICER data showed that the source exhibited large amplitude flares on timescales as short as ~10 s, the largest of which had a peak count rate of ~1000 cts s⁻¹ and lasted roughly 15 s (Ludlam et al. 2018). The NICER spectra were divided into high (>100 cts s⁻¹), moderate (20–100 cts s⁻¹), and low (<20 cts s⁻¹) intensity intervals, and fit with an absorbed thermal disk plus power-law model. The best-fit models found fairly low disk temperatures of ~0.2–0.3 keV, which are rather typical for BHs in the hard state (see, e.g., Reis et al. 2010; Reynolds & Miller 2013), while the power-law component was found to be very hard, Γ ~ 1 (Ludlam et al. 2018). Additionally, the NICER spectra showed both Fe L and K reflection features.

J1858’s longer wavelength counterpart was first detected as a variable UV source by the Ultraviolet and Optical Telescope on Swift (Swift-UVOT), and it was found that the source was coincident with a previously detected UKIRT Infrared Deep Sky Survey and Pan-STARRS source (Kennea & Krimm 2018). Additional optical follow-up found that the source had brightened by ~2.5 magnitudes over the source’s cataloged Pan-STARRS *r'* magnitude (Vasilopoulos et al. 2018). The source has also been detected in radio by the Arcminute Microkelvin Imager Large Array and appears to be variable, having a flux density of 300–600 μJy at 15.5 GHz (Bright et al. 2018).

The source’s outburst is still ongoing to date<fnref n="9" /> (i.e., 2019 October), having a mean Swift-BAT (Krimm et al. 2013) flux

of ≈14 mCrab (see Figure 1). Further, the source has now been observed six times by NuSTAR (see Figure 1). Here we report the results of the analysis of the first of these NuSTAR observations. In Section 2, we describe the details of the NuSTAR and NICER observations and data reduction, then in Section 3 we discuss the data analysis and results. Next, in Section 4, we discuss the physical parameters of this system and compare them to other similar systems, namely, V404 Cyg and V4641 Sgr. Finally, we summarize our findings in Section 5.

## 2. Observations and Data Reduction

### 2.1. NuSTAR

J1858 was observed with NuSTAR (Harrison et al. 2013) on 2018 November 3 (MJD 58425.28, obsID 80401317002) for ~52 ks after correcting for dead time. The data were reduced using the NuSTAR Data Analysis Software (NuSTARDAS) package version 1.8.0 and the 20181022 version of the Calibration database (CALDB). First, the photon arrival times were corrected to the solar system barycenter using *nupipe-line*, which also includes a clock correction from the CALDB to account for NuSTAR’s clock drift.<fnref n="10" /> Then, the energy spectra and light curves of the source were extracted from both the focal plane detector modules (FPM) A and B detectors using an *r* = 90" circle centered on the source. The corresponding background spectra and light curves were extracted from a source-free circular region (*r* ≈ 50") on the same detector chip as the source. Prior to fitting, the energy spectra were grouped to have a signal-to-noise ratio of at least five in each energy bin. We fit all X-ray energy spectra with XSPEC version 12.10.1 (Arnaud 1996), accounting for interstellar absorption with v2.3.2 of the Tuebingen-Boulder interstellar matter absorption model, *tbnew*,<fnref n="11" /> with solar abundances adopted from Wilms et al. (2000). All uncertainties in this paper are reported at the 90% confidence level unless otherwise noted.

<footnote n="9">⁹ See https://swift.gsfc.nasa.gov/results/transients/weak/SWIFTJ1858.6-0814/.</footnote>

<footnote n="10">¹⁰ See http://www.srl.caltech.edu/NuSTAR_Public/NuSTAROperationSite/clockfile.php.</footnote>

<footnote n="11">¹¹ See http://pulsar.sternwarte.uni-erlangen.de/wilms/research/tbabs/.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** FPMA light curve of Swift J1858 with 100 s time bins. We show only the FPMA light curve for clarity. The inset shows a zoom in with 1 s time bins over a ~1 ks span (denoted by the vertical red lines) that includes the two large flares that occurred during the simultaneous *NuSTAR* and *NICER* observations (see Sections 2.2 and 3.2.2). *NICER*'s count rates in the 0.5–7 keV energy range, shown as red points in the inset, are roughly a factor of 10 larger than *NuSTAR*'s count rates. The green lines in the inset show the time intervals in which each flare's energy spectra were extracted (see Section 3.2.2).</figcaption></figure>

## 2.2. NICER

NICER (Gendreau et al. 2012) observed J1858 twice simultaneously with the NuSTAR observation reported here. The two NICER observations had exposure times of ~5.7 ks and 5.2 ks, for obsIDs 1200400103 and 1200400104, respectively. The NICER data were reduced following the standard processing and filtering procedures using the NICER-DAS (V005) software package and the xti20190520 version of the CALDB. We also excluded data from the two detectors that are known to exhibit increased detector noise.<fnref n="12" /> We use the latest response functions (i.e., nixtiaveonax-is20170601v002.arf and nixtiref20170601v001. rmf) for the spectral analyses of the NICER data performed in this paper. To minimize the effects of residuals that still remain in the NICER response functions, and to further minimize the effects of detector noise, we restrict our analysis to the 0.5–7 keV energy range.<fnref n="13" />

The NICER data are primarily used in this paper to constrain the soft part of J1858’s X-ray spectrum. Since J1858 shows significant spectral evolution during its flares, we only used bright flares simultaneously observed by NICER and NuSTAR. During the first observation, NICER unfortunately observed J1858 while the source was occulted by the Earth for NuSTAR, so there is only ~1400 s of strictly simultaneous data, none of which contains any particularly bright flares. Therefore, we do not use any of the data from the first observation. However, the second NICER observation overlapped with the NuSTAR

observation for ~4.5 ks, and caught two of the brightest flares observed by NuSTAR, occurring ~570 s apart (see the inset in Figure 2). These flares are referred to as flare 1 and flare 2, hereafter, with flare 1 occurring first, and flare 2 being brighter.

# 3. Results

## 3.1. *Variability and Timing*

To study the variability of Swift J1858 on different timescales, we have produced light curves using a number of different time bins (i.e., 1 s, 10 s, 100 s, 1 ks, 5 ks), which have also been corrected for NuSTAR’s various detector effects (e.g., dead time, point-spread function, vignetting). These light curves reveal that the source was highly variable throughout the NuSTAR observation, showing large amplitude flares (see Figure 2). These flares typically lasted between ~10–100 s with the largest having a peak count rate ~50 times higher than the source’s average count rate (see the inset in Figure 2). Throughout the observation, the source also showed changes in its hardness ratio, defined here as the 10–79 keV count rate divided by the 3–10 keV count rate (i.e., $R_{10-79\text{keV}}/R_{3-10\text{keV}}$). The FPMA and B averaged hardness-intensity diagram (HID) shows that the source is softer during the flaring episodes, while spanning a range of hardness ratios during the nonflaring periods (see Figure 3). The flares are also observed across NuSTAR’s entire bandpass, but are most strongly observed in the 3–10 keV energy band (see Figure 4).

To better understand the timing properties of the source, and to look for differences between the flaring and nonflaring periods, we use the HID to define two distinct modes of the source, which we designate as either “flaring” or “nonflaring.” The source is considered to be in the flaring mode (black points

<footnote n="12">¹² See https://heasarc.gsfc.nasa.gov/docs/nicer/data_analysis/nicer_analysis_tips.html for additional details.</footnote>

<footnote n="13">¹³ See, for example, https://heasarc.gsfc.nasa.gov/docs/nicer/data_analysis/nicer_analysis_tips.html.</footnote>
<figure><img src="image_4.png" /></figure>

**Figure 3.** NuSTAR hardness-intensity diagram (HID) with 1 ks time bins averaged over the FPMA and B detectors. The solid black line shows the distinction between flaring (black points above the line) and nonflaring (red and blue points below the line) time intervals. The dashed vertical black line shows the further distinction between nonflaring soft (NF<sub>S</sub>; red points to the left of the line) and nonflaring hard (NF<sub>H</sub>; blue points to the right of the line) time intervals (see Sections 3.1 and 3.2.1). The gray cross shows where the interval containing the two flares simultaneously observed with NuSTAR and NICER and discussed in Section 3.2.2 lands on the HID.

<figure><img src="image_5.png" /></figure>

**Figure 4.** NuSTAR FPMA energy resolved (3–10 keV, red; 10–20 keV, black; 20–79 keV, blue) light curves of flare 2 (see Sections 2.2 and 3.2.2) with 1 s time bins. The flares are detected across all energy bands, but have the largest amplitudes at soft (i.e., 3–10 keV) energies.

in Figure 3) when the 3–79 keV energy band count rate and hardness ratio, averaged over 1 ks time intervals, are above the line defined as $R_{3\text{–}79\text{keV}} = 4 \times R_{10\text{–}79\text{keV}}/R_{3\text{–}10\text{keV}}$ (the black solid line shown in Figure 3), while the nonflaring mode includes all of the points below this line (i.e., red and blue points shown in Figure 3). The flaring and nonflaring light curves, after being split, consisted of exposure times of ~25 ks and ~26 ks, respectively.

To characterize the observed variability and to search for possible quasi-periodic oscillations (QPOs), we constructed power density spectra<fnref n="14" /> (PDS) from the NuSTAR event lists<fnref n="15" /> using the Stingray python package (Huppenkothen et al. 2019).

First, light curves with a 4 ms binning were produced from the barycenter corrected NuSTAR event files in the 3–79 keV energy range. We also removed 100 s from the beginning and 200 s from the end of each good time interval (GTI) to eliminate any possible effects from an increased background that may occur near the borders of GTIs (see, e.g., Section 5 in Bachetti et al. 2015).

The PDS were produced for the full, flaring, and nonflaring time intervals spanning a 0.001–125 Hz frequency range and were averaged over 1 ks time segments. The PDS were geometrically rebinned by a factor of 1.08 (see Figure 5). The error bars in the PDS become very large at frequencies ≳1 Hz, so these points were excluded from our analysis after verifying no significant peak is detected in the PDS at these high frequencies. We then simultaneously fit a single, zero frequency centered Lorentzian model to the FPMA and FPMB PDS. The single-Lorentzian model fits the data reasonably well, with the exception of the flaring data (see Figure 5, right column), and we do not find evidence of any remarkable features (e.g., QPO, orbital periodicity) in the PDS. The best-fit widths of the Lorentzian models, and their $1\sigma$ uncertainties, are $(1.27 \pm 0.06) \times 10^{-2}$ Hz, $(1.8 \pm 0.1) \times 10^{-2}$ Hz, and $(1.4 \pm 0.1) \times 10^{-2}$ Hz, for the full, flaring, and nonflaring time intervals. We derived the fractional rms amplitudes and their $1\sigma$ uncertainties by integrating the best-fit Lorentzian models, which give 147% ± 3%, 135% ± 4%, and 129% ± 4% for the full, flaring, and nonflaring time intervals, respectively. These large rms fractional values are indicative of the large flux fluctuations exhibited by the source.

The single-Lorentzian model is a relatively poor fit to the flaring PDS ($\chi^2_{\text{red}} = 1.54$), so we also fit a model including a second zero frequency centered Lorentzian. This model provides a substantially better fit ($\chi^2_{\text{red}} = 0.98$), reducing the chi-squared from $\chi^2 = 172.3$ to $\chi^2 = 107.3$ (or $\Delta\chi^2 = 65$) for two fewer degrees of freedom (dof). The best-fit widths for this two-Lorentzian model and their $1\sigma$ uncertainties are $(1.51_{-0.09}^{+0.10}) \times 10^{-2}$ Hz and $(69_{-16}^{+12}) \times 10^{-2}$ Hz, for the low (dominating between $10^{-3}$ and $10^{-1}$ Hz) and high (dominating between $10^{-1}$ and 1 Hz) frequency Lorentzians, respectively. The fractional rms amplitudes and their $1\sigma$ uncertainties are derived in the same way as described in the previous paragraph and are 135% ± 3% and 27% ± 2% for the low- and high-frequency Lorentzians, respectively.

## 3.2. X-Ray Spectra

### 3.2.1. Relativistic Reflection

To characterize the spectral differences between J1858's flaring and nonflaring intervals, we extract the energy spectra from three different modes. The first mode is the same as the "flaring" mode defined above in Section 3.1 and encompasses the black data points above the solid black line shown in Figure 3. Since the nonflaring mode, defined as points below the solid black line in Figure 3, spans a broad range of hardness ratios, we further divide this mode by the dashed black vertical line shown in Figure 3 into nonflaring soft ($R_{10\text{–}79\text{keV}}/R_{3\text{–}10\text{keV}}<1.5$) and nonflaring hard ($R_{10\text{–}79\text{keV}}/R_{3\text{–}10\text{keV}}>1.5$) modes (i.e., shown as red and blue points in Figure 3, respectively). The HID bin size of 1 ks was chosen to ensure that the hardness ratio error bars were small enough to confidently differentiate the nonflaring data points

<footnote n="14">¹⁴ This observation was not significantly affected by dead time, having a dead time fraction of <10%. Therefore, we used the typical PDS and not the cross-spectrum (see, e.g., Bachetti et al. 2015).</footnote>

<footnote n="15">¹⁵ We note that Stingray constructs light curves using the event lists, and therefore does not correct the light curves for the various detector effects mentioned above.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Left: NuSTAR power spectra of Swift J1858 for the full, flaring, and nonflaring time intervals. The uncertainties in the rms are reported at the $1\sigma$ level. Right: $\Delta\chi^2$ residuals for the best-fit single, zero frequency centered Lorentzian models.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Left: unfolded NuSTAR energy spectra of J1858 from the full duration of the observation. Right: data to model ratio for the flaring (red and black), nonflaring hard (green and dark blue), and nonflaring soft (cyan and purple) spectra of J1858 fitted with an exponentially cutoff power law in the 3–4, 8–10, and 30–79 keV energy ranges and then plotted after reintroducing the data in the full 3–79 keV energy range. All spectra have been rebinned for easier visualization.</figcaption></figure>
into the hard and soft modes.<fnref n="16" /> After dividing the data in this way, exposures of ~25 ks, ~11 ks, and ~15 ks remained for the flaring, nonflaring hard, and nonflaring soft modes, respectively. For the remainder of the paper, we denote these three modes as flare, $NF_S$, and $NF_H$ for the flaring, nonflaring soft, and nonflaring hard modes, respectively.

The NuSTAR energy spectra in the 3–79 keV range spanning the entire observation (i.e., not split into the three modes) are shown in Figure 6 and exhibit a number of features typical of accreting BH systems. These features include an excess of emission at energies between 5 and 7 keV, typical of iron K features, an absorption edge around 7 keV, and a broad Compton hump above 10–15 keV. To highlight the differences between the X-ray spectra during the three intervals, we fit them with an exponentially cutoff power-law model in the 3–4, 8–10, and 30–79 keV energy bands (i.e., excluding the 4–8 and 10–30 keV energy ranges). These energy bands are chosen as they provide relatively unbiased access to the underlying continuum by avoiding the strongest spectral features previously mentioned (see, e.g., Walton et al. 2017). We then plot the ratio of the data to the folded model, after re-including the data in the full 3–79 keV energy range, in the right panel of Figure 6. Clear differences can be seen between the spectra from the three time intervals. For instance, the iron lines observed in the $NF_S$ and $NF_H$ spectra are more strongly peaked around 6.4 keV, while also having a more pronounced absorption edge at ~7 keV and Compton hump above ~10 keV when compared to the flaring spectrum. All three spectra show evidence of a red wing, i.e., a broadening of the iron line extending to lower energies. The red wing provides strong evidence that this emission is coming from a relativistically broadened reflection of photons off the inner-most regions of the accretion disk (see, e.g., Fabian et al. 1989; Laor 1991). Lastly, the Compton hump emission is most pronounced in the $NF_H$ spectra.

Given the broad reflection features evident in the X-ray spectra, we use the collection of RELXILL models (version 1.2.0; García et al. 2014) to fit them. The RELXILL models combine the XILLVER (García & Kallman 2010) reflection model with the RELCONV relativistic convolution model, which captures the relativistic effects due to the emitting material’s close proximity to the BH (Dauser et al. 2010). In particular, we use the RELXILLPCP model to fit the X-ray spectra of J1858. This model uses the thermally Comptonized continuum model nthComp (Zdziarski et al. 1996; Życki et al. 1999) for the input continuum spectrum and assumes a lamp-post geometry (i.e., a point source directly above the spin axis of the BH) for the illuminating X-ray source.

RELXILLPCP is characterized by several physical parameters intrinsic to the BH binary. This includes the inclination of the inner accretion disk, *i*, the iron abundance of the accreted material, $A_{Fe}$, and the spin of the BH, $a^*$. Since these parameters are not expected to change during the duration of our observation, they are linked between all spectra during the fitting procedure. This model also contains parameters that can change between the different time intervals, including the photon index of the power-law emission incident on the

<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) NuSTAR energy spectra of J1858 from the flaring (red and black), nonflaring hard (green and blue), and nonflaring soft (cyan and magenta) time intervals. The best-fit model components from the *relxill_lp cp +xillver cp +diskbb* model are plotted as dashed lines in each corresponding color. (b) Residuals of the data to the model for the best-fit model. (c) Residuals of the data to the model for the *relxill_lp cp +xillver cp* model. Note that large residuals appear at soft X-ray energies. (d) Residuals of the data to the model for the *relxill_lp cp* model. Large residuals appear at soft X-ray energies and around the iron line complex. All spectra have been rebinned for easier visualization.</figcaption></figure>

accretion disk, $\Gamma$, the temperature of the electrons in the corona, $kT_e$, the height of the source above the BH that is irradiating the accretion disk, $h$, and the ionization state of the iron in the accretion disk, $\log\zeta$.<fnref n="17" /> These parameters are allowed to vary for the flaring, $NF_S$, and $NF_H$ spectra during our fits. The reflection fraction, $R_{\text{ref}}$, is another parameter in this model and is defined as the ratio of the amount of light from the primary source that is emitted toward the disk versus the amount that escapes to infinity (Dauser et al. 2016). For the lamp-post geometry, the RELXILL package offers the option to calculate the reflection fraction in a self-consistent way using relativistic ray-tracing, which we take advantage of for our fits. The model also contains the inner radius of the accretion disk, $R_{\text{in}}$, which we allow to vary across all three spectral modes.

The best-fit RELXILLPCP model to the data has a reduced chi-squared of $\chi^2/\nu = 3243/2933$ and still shows large residuals both at low energies and near the iron complex around 6.4 keV (see Figure 7(d)). The iron line features are more strongly peaked in the residuals of the $NF_S$ and $NF_H$ spectra, and are possibly due to reflection from distant cold material. This suggests that additional model components are

<footnote n="16">¹⁶ While it is difficult to accurately separate the nonflaring soft and hard modes with 100 s binning, it is still possible to separate the flaring mode data. We have carried out this exercise and found that the best-fit model is consistent with the best-fit flaring model shown in Table 1 (within uncertainties), with the exception of slightly higher normalizations.</footnote>

<footnote n="17">¹⁷ Here $\xi = 4\pi F_X/n$, where $F_X$ is the ionizing flux incident on the accretion disk, and $n$ is the density of the material in the disk.</footnote>
needed to adequately fit the spectra. To account for the excess
near the iron complex, we add a neutral (i.e., logξ = 0) XILLVERCP component to the model. The XILLVERCP model was chosen because its continuum emission model (i.e., nthComp) is the same as that used in the RELXILLPCP model. Therefore, we tie the parameters shared by both models together as the same source should be illuminating both the accretion disk and distant reflector.<fnref n="18" /> Additionally, we fix the ionization of the distant reflector to log ξ = 0 because the narrow part of the residual is peaked near the 6.4 keV Fe Kα line, implying that it is likely coming from neutral iron. We also assume that all of the emission coming from this component of the model is reflected (i.e., denoted in the XILLVERCP model by setting the reflection fraction to -1). We allow the normalization of the XILLVERCP component to vary between the flaring, NF<sub>H</sub>, and NF<sub>S</sub> spectra. Following this addition, the reduced chi-squared, χ²/ν = 3051/2930, improved by >190 for three fewer dof. However, the excess at soft X-ray energies still remains (see Figure 7(c)), so we add a multitemperature blackbody (diskbb; Mitsuda et al. 1984) to account for thermal emission from the accretion disk, which is likely causing this excess. We allow the accretion disk temperature and normalization of the DISKBB component to vary between the flaring, NF<sub>H</sub>, and NF<sub>S</sub> spectra. Hence, the complete model is `const*tbnew*(diskbb+relxill.lp cp+xillver cp)`. This additional component further improved the fit, leading to an improvement of ~120 in the reduced chi-squared (χ²/ν = 2928/2924) for six fewer dof. It is interesting to note that about half of the reduction in the chi-squared comes from the addition of the thermal component to the flaring spectrum. The best-fit model parameters for J1858 can be found in Table 1, while the best-fit spectra and their residuals are shown in Figures 7(a) and (b), respectively.

Using the best-fit model, we have also tried to require the absorbing column density to be the same across the spectra of all three source modes. This led to a best-fit absorbing column density of $N_H = 18 \times 10^{22} \text{ cm}^{-2}$ and a reduced chi-squared of $\chi^2/\nu = 2958/2926$. This fit has an additional $\Delta\chi^2 = 30$ for two additional dof. We disfavor this model because Figure 6 shows clear differences in the strength of the absorption edge at ~7 keV, strongly suggesting that there is additional intrinsic absorption during the nonflaring mode compared to the flaring mode (see, e.g., Walton et al. 2017).

3.2.2. *Joint NICER and NuSTAR Fits*

The relativistic reflection fits to the *NuSTAR* energy spectra
place constraints on several spectral features that are more
prominently observed in the soft X-ray band (e.g., large
intrinsic absorbing column density, *diskbb* temperature, and
normalization). Therefore, to verify that the values of these
features derived from the *NuSTAR* spectra are reasonable, we
use the simultaneous *NICER* observations.<fnref n="19" /> To minimize the
effects of *NICER*'s background contribution to the X-ray

energy spectra we chose to focus on short (~100 s) flares
simultaneously observed by *NuSTAR* and *NICER*.

The energy spectra from the first flare, which reached peak
count rates of ~770 cts s⁻¹ and ~85 cts s⁻¹ for NICER and
NuSTAR, respectively, were extracted from a 60 s time interval
around the flare. These spectra contain a total of ~950 counts
(cts) in each NuSTAR focal plane module and ~11,100 cts in
NICER. The energy spectra of the second flare, reaching peak
NuSTAR and NICER count rates of ~135 cts s⁻¹ and
~860 cts s⁻¹, were extracted from a 100 s window containing
the flare (see the inset in Figure 2). This flare's spectrum
contained ~2300 cts in each NuSTAR focal plane module and
~17,000 cts in NICER. These flares have a 1 ks average
NuSTAR count rate of ~10 cts s⁻¹ and a hardness ratio of ~0.6
(see the black cross in Figure 3). Prior to fitting the spectra we
grouped them to have 100 cts per energy bin for the NICER
spectra and 50 cts per energy bin for the NuSTAR spectra.
Unfortunately, there appears to be strong systematic residuals
in the NICER spectra around 0.5 keV, so we avoid energies
below 0.6 keV in our fits.

To check that the large intrinsic absorption and thermal
component observed by *NuSTAR* is consistent with the *NICER*
spectra, we jointly fit the *NICER* and *NuSTAR* spectra with the
best-fit flaring model shown in Table 1. Due to the relatively
small number of counts in the *NuSTAR* energy spectra for
these short duration flares we freeze all of the parameters
of the model except for the absorbing column density, the
normalizations for the three model components (i.e.,
*rellxill.lp cp*, *xillver.lp cp*, *diskbb*), and the
cross-calibration constants.<fnref n="20" /> This model produces a very
poor fit to the spectra of flares 1 and 2, having reduced chi-
squared values of $\chi^2/\nu = 200/133$ and $\chi^2/\nu = 675/224$,
respectively, and giving small absorbing column densities of
$N_H \approx 4 \times 10^{21} \text{ cm}^{-2}$. However, the fits substantially improve
if a partially covering absorber (pcfabs) is added to the
model and allowed to vary ($\chi^2/\nu = 129/131$, and
$\chi^2/\nu = 209/222$, for flares 1 and 2, respectively). The best-
fit parameters of the model to the two flares are shown in
Table 2, while the best-fit spectra and residuals for flare 2 are
shown in Figure 8. Most notably, the *diskbb* normalization is
consistent with the *NuSTAR*-only fits for both flares 1 and 2,
suggesting that the thermal component required by the fits to
the *NuSTAR*-only data is confirmed by the *NICER* data.

4. Discussion

4.1. *Spectral and Timing Features*

The best-fit model parameters of the physical values inherent to the binary are an iron abundance of the accreting material, $A_{\text{Fe}} = 1.00_{-0.2}^{+0.4}$, and a 90% upper limit on the inclination of the inner accretion disk $i < 23^{\circ}$. We also calculate the $3\sigma$ upper limit on the inclination, finding $i < 29^{\circ}$. This suggests that the disk is viewed almost face on. We also find 90% upper limits for the inner radius of the accretion disk, $R_{\text{in}}$, of $<8 \text{ r}_{\text{ISCO}}$, $<6 \text{ r}_{\text{ISCO}}$, and $<5 \text{ r}_{\text{ISCO}}$ for the flaring, $\text{NF}_{\text{H}}$, and $\text{NF}_{\text{S}}$ spectra, respectively, which are consistent with $R_{\text{in}} = R_{\text{ISCO}}$. It is important to note that the spin and $R_{\text{in}}$ are degenerate with each other. Therefore, to calculate a lower limit on the spin of the BH, we set $R_{\text{in}}$ to be at $r_{\text{ISCO}}$, refit the model, and calculate the lower limit on the spin. After fitting, we find that all best-fit

<footnote n="18">¹⁸ Realistically, the distant reflector will see the illuminating source gravitationally redshifted. However, due to the relatively large illuminating source heights, the gravitational redshift is small (≤0.25), making this effect negligible.</footnote>

<footnote n="19">¹⁹ The Neil Gehrels Swift Observatory's X-ray Telescope (XRT) also observed J1858 simultaneously with NuSTAR for ~400 s (obsID 00010955002). However, due to the short duration of the observation, no bright flares from J1858 were detected by Swift-XRT leaving only ~250 total counts.</footnote>

<footnote n="20">²⁰ For simplicity, the inner radius of the accretion disk is frozen at the ISCO (i.e., $R_{\text{in}} = R_{\text{ISCO}}$) for these fits.</footnote>
<table><caption>**Table 1**

Best-fit Parameters of the `const*tbnew(diskbb+relxill.lp cp+xillver cp)` Model Simultaneously Fit to the Flaring, Nonflaring Hard ($NF_H$), and Nonflaring Soft ($NF_S$) *NuSTAR* Spectra of Swift J1858</caption><thead><tr><th>Model Component</th><th>Parameter</th><th>Units</th><th>Flare</th><th>NF<sub>H</sub></th><th>NF<sub>S</sub></th></tr></thead><tbody><tr><td>Constant</td><td>...</td><td>FPMB/FPMA</td><td>1.020 ± 0.007</td><td>1.02 ± 0.02</td><td>1.02 ± 0.02</td></tr><tr><td>TBNEW</td><td>N<sub>H</sub></td><td>10<sup>22</sup> cm<sup>-2</sup></td><td>14 ± 2</td><td>42<sup>+8</sup><sub>-9</sub></td><td>30<sup>+8</sup><sub>-7</sub></td></tr><tr><td>DISKBB</td><td>T<sub>in</sub><br>Norm</td><td>keV<br>10<sup>3</sup> × (R<sub>in,km</sub><sup>a</sup>/D<sub>10</sub><sup>b</sup>)<sup>2</sup>cosi</td><td>0.36<sup>+0.03</sup><sub>-0.05</sub><br>7<sup>+20</sup><sub>-4</sub></td><td>0.31 ± 0.03<br>62<sup>+97</sup><sub>-51</sub></td><td>0.35 ± 0.04<br>8<sup>+13</sup><sub>-6</sub></td></tr><tr><td rowspan="10">RELXILLPCP</td><td>i<sup>c</sup></td><td>degrees</td><td>&lt;23°</td><td>...</td><td>...</td></tr><tr><td>A<sub>Fe</sub><sup>c</sup></td><td>solar</td><td>1.0<sup>+0.4</sup><sub>-0.2</sub></td><td>...</td><td>...</td></tr><tr><td>Γ</td><td>...</td><td>1.50 ± 0.03</td><td>1.41<sup>+0.10</sup><sub>-0.08</sub></td><td>1.49<sup>+0.1</sup><sub>-0.03</sub></td></tr><tr><td>kT<sub>e</sub></td><td>keV</td><td>15.0<sup>+0.9</sup><sub>-0.7</sub></td><td>14 ± 2</td><td>16 ± 2</td></tr><tr><td>h</td><td>r<sub>g</sub></td><td>13<sup>+7</sup><sub>-5</sub></td><td>6<sup>+6</sup><sub>-4</sub></td><td>5<sup>+3</sup><sub>-2</sub></td></tr><tr><td>a<sup>*c</sup></td><td>...</td><td>&gt;0.0</td><td>...</td><td>...</td></tr><tr><td>log ξ</td><td>log (erg cm s<sup>-1</sup>)</td><td>3.54<sup>+0.10</sup><sub>-0.13</sub></td><td>3.0 ± 0.3</td><td>3.2<sup>+0.1</sup><sub>-0.2</sub></td></tr><tr><td>r<sub>in</sub></td><td>r<sub>ISCO</sub></td><td>&lt;8</td><td>&lt;6</td><td>&lt;5</td></tr><tr><td>R<sub>refl</sub></td><td>...</td><td>1.3<sup>d</sup></td><td>1.6<sup>d</sup></td><td>1.5<sup>d</sup></td></tr><tr><td>Norm</td><td>10<sup>-4</sup></td><td>6.6<sup>+1.2</sup><sub>-0.7</sub></td><td>2.5<sup>+3.3</sup><sub>-0.9</sub></td><td>2.6<sup>+1.7</sup><sub>-0.8</sub></td></tr><tr><td rowspan="3">XILLVERPCP</td><td>log ξ</td><td>log (erg cm s<sup>-1</sup>)</td><td>0.0<sup>e</sup></td><td>0.0<sup>e</sup></td><td>0.0<sup>e</sup></td></tr><tr><td>R<sub>refl</sub></td><td>...</td><td>-1<sup>e</sup></td><td>-1<sup>e</sup></td><td>-1<sup>e</sup></td></tr><tr><td>Norm</td><td>10<sup>-4</sup></td><td>2.0 ± 0.7</td><td>2.0 ± 0.9</td><td>1.8<sup>+0.7</sup><sub>-0.5</sub></td></tr><tr><td>Observed Flux</td><td>3.0–79 keV</td><td>10<sup>-10</sup> erg cm<sup>-2</sup> s<sup>-1</sup></td><td>6.79 ± 0.02</td><td>2.10 ± 0.02</td><td>1.88 ± 0.02</td></tr><tr><td>Unabsorbed Flux</td><td>3.0–79 keV</td><td>10<sup>-10</sup> erg cm<sup>-2</sup> s<sup>-1</sup></td><td>7.40 ± 0.02</td><td>2.56 ± 0.03</td><td>2.19 ± 0.02</td></tr><tr><td>χ<sup>2</sup>/dof</td><td>2928/2924</td><td>1.001</td><td></td><td></td><td></td></tr></tbody></table>

**Notes.**

<sup>a</sup> The apparent inner disk radius in units of km.

<sup>b</sup> Distance to the source in units of 10 kpc.

<sup>c</sup> Parameter is tied across Flare, $NF_S$, and $NF_H$ spectral models. The best-fit value is given in the Flare column.

<sup>d</sup> Calculated self-consistently by the RELXILLPCP model.

<sup>e</sup> Fixed value.

parameters in this updated model are consistent (within the 90% uncertainties) with those found previously and shown in Table 1. This model gives a 90% confidence lower limit on the spin of $a^* > 0.0$ and a $3\sigma$ lower limit of $a^* > -0.8$ (see Figure 9).

One of the most striking features of the best-fit model to the *NuSTAR* energy spectra is the absorbing column density, which is very large even in the flaring spectrum (i.e., $N_H = (14 \pm 2) \times 10^{22} \text{ cm}^{-2}$). This absorption is a factor of ~2–3 larger for the $NF_H$ ($N_H = 42_{-9}^{+8} \times 10^{22} \text{ cm}^{-2}$) and $NF_S$ ($N_H = 30_{-7}^{+8} \times 10^{22} \text{ cm}^{-2}$) spectra. However, this absorption cannot fully account for the change in the source’s flux (see unabsorbed fluxes in Table 1), implying that the source must also be intrinsically variable. Furthermore, the absorption appears to be anticorrelated with the source’s intrinsic (i.e., unabsorbed) flux between the flaring, $NF_H$, and $NF_S$ modes. On the other hand, the intrinsic absorption found in our joint NICER and *NuSTAR* fits is about a factor of two smaller than the intrinsic absorption found from the fits to the averaged flaring *NuSTAR* spectra. Assuming that the intrinsic absorption is caused by dense clouds of material intersecting the observer’s line of sight (LOS), similar to V404 Cyg (see, e.g., Motta et al. 2017a; and Section 4.3 for more details), this difference is somewhat expected. This is because flaring modes are defined in the *NuSTAR* data using the 1 ks binned light curves. However, given the rapid variability of the source, this binning will inevitably include some short periods when the

source is not flaring.<fnref n="21" /> Therefore, since these nonflaring periods are correlated with a larger intrinsic absorption, then the 1 ks “averaged” absorbing column density would be expected to be larger than the absorbing column density observed during a shorter ~100 s flare. Furthermore, in a later Swift-XRT observation (i.e., not the one coincident with our *NuSTAR* observation) it was found that energy spectrum at soft X-ray energies could also be adequately fit by a similar partially covered thermal plasma plus a power-law model, with the partially covering absorber having a significant absorbing column density of $N_H = 17_{-0.5}^{+0.6} \times 10^{22} \text{ cm}^{-2}$ (Reynolds et al. 2018), which is consistent with the value found from the fits to the flaring *NuSTAR* energy spectra.<fnref n="22" /> Additional support for the fact that most of this absorption must be intrinsic to the binary itself, is that the total Galactic HI absorbing column density in this direction is only $N_H \approx 1.8 \times 10^{21} \text{ cm}^{-2}$ (HI4PI Collaboration et al. 2016). We also note that the Galactic HI absorbing column density is consistent (within a factor of about two) with the value found for the interstellar absorption in our joint NICER+$NuSTAR$ spectral fits (i.e., $N_H \approx 4 \times 10^{21} \text{ cm}^{-2}$).

Although no QPOs or orbital modulation were found in our X-ray light curves or PDS they did show large amounts of variability. The full PDS had a large fractional rms amplitude

<footnote n="21">²¹ See the inset in Figure 2 for an example of a typical 1 ks flaring bin viewed on shorter 1 s timescales.</footnote>

<footnote n="22">²² The *Swift* spectra were not divided into flaring and nonflaring modes for this spectral fit.</footnote>
<table><caption>**Table 2**

Best-fit Parameters of the `const*tbnew*pcfabs(diskbb+relxill.lp cp+xillver cp)` Model Fit to Flares 1 and 2, which Were Simultaneously Observed by *NuSTAR* and *NICER*</caption><thead><tr><th>Model Component</th><th>Parameter</th><th>Units</th><th>Flare 1</th><th>Flare 2</th></tr></thead><tbody><tr><td>Constant</td><td>...</td><td>FPMB/FPMA</td><td>1.12 ± 0.08</td><td>1.00 ± 0.05</td></tr><tr><td>Constant</td><td>...</td><td>NICER/FPMA</td><td>1.18 ± 0.08</td><td>1.07 ± 0.05</td></tr><tr><td>TBNEW</td><td>N<sub>H</sub></td><td>10<sup>22</sup> cm<sup>-2</sup></td><td>0.41 ± 0.03</td><td>0.39 ± 0.02</td></tr><tr><td>PCFABS</td><td>N<sub>H</sub><br>Covering Fraction</td><td>10<sup>22</sup> cm<sup>-2</sup></td><td>5.3<sup>+1.1</sup><sub>-0.9</sub><br>0.78<sup>+0.09</sup><sub>-0.21</sub></td><td>6.2<sup>+0.7</sup><sub>-0.5</sub><br>0.84<sup>+0.05</sup><sub>-0.09</sub></td></tr><tr><td>DISKBB</td><td>Norm</td><td>10<sup>3</sup> × (R<sub>in,km</sub>/D<sub>10</sub>)<sup>2</sup> cos i</td><td>4<sup>+4</sup><sub>-3</sub></td><td>5<sup>+3</sup><sub>-2</sub></td></tr><tr><td>RELXILLPCP</td><td>Norm</td><td>10<sup>-4</sup></td><td>24 ± 1</td><td>32 ± 2</td></tr><tr><td>XILLVERLPCP</td><td>Norm</td><td>10<sup>-4</sup></td><td>&lt;2.5</td><td>&lt;9</td></tr><tr><td>Observed Flux</td><td>0.6–79 keV</td><td>10<sup>-9</sup> erg cm<sup>-2</sup> s<sup>-1</sup></td><td>2.52 ± 0.04</td><td>3.46 ± 0.04</td></tr><tr><td>Unabsorbed Flux</td><td>0.6–79 keV</td><td>10<sup>-9</sup> erg cm<sup>-2</sup> s<sup>-1</sup></td><td>3.51 ± 0.05</td><td>4.78 ± 0.05</td></tr><tr><td>χ<sup>2</sup>/dof</td><td></td><td></td><td>129/130</td><td>209/222</td></tr></tbody></table>

Note. Model parameters not listed here were fixed to the best-fit values shown in Table 1.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Best-fit `const*tbnew*pcfabs(diskbb+relxill.lp cp+xillver cp)` NICER (green) and NuSTAR (black and red for FPMA and FPMB, respectively) energy spectra. The NICER data were limited to the 0.6–7 keV band (see Section 3.2.2). The best-fit model parameters are shown in Table 2.</figcaption></figure>

of 1.47 (see Section 3.1), capturing the large flares observed in the light curves. It should be noted that there has been a report of a relatively low frequency QPO in a later NuSTAR observation at a frequency of $2.7 \times 10^{-3}$ Hz (i.e., ~364 s; Hare et al. 2019), suggesting that this QPO unfortunately formed while the source was Sun constrained for NuSTAR. The analysis of this later NuSTAR observation and potential QPO will be presented in a future publication.

## 4.2. Optical Counterpart and Distance

An optical/UV counterpart coincident with J1858's Swift-XRT position was first detected by Swift-UVOT, having a UVW2 Vega magnitude of $17.38 \pm 0.08$ (Kennea & Krimm 2018). Additional follow-up observations found that this optical source was highly variable, varying between $r'$ magnitudes of ~18.4–16.3, with variations as large as 1 magnitude on timescales of 2 minutes to as short as <5 s

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Constraints on the lower limit of the spin parameter of J1858 using XSPEC's `steppar` command. The dashed black lines show the 90% and 3σ confidence intervals for a single parameter. The solid black line shows the $\Delta \chi^2$ of the spin parameter if the inner radius of the accretion disk is frozen at the ISCO.</figcaption></figure>

(Baglio et al. 2018; Paice et al. 2018; Rajwade et al. 2018; Vasilopoulos et al. 2018). Further, the slower, but brighter flares were found to be more blue, while the fast flares were found to be more red (Paice et al. 2018). Interestingly, the optical counterpart was bright enough to be detected by Pan-STARRs, with an $r'$ magnitude of $19.97 \pm 0.03$ prior to the source going into outburst$^{23}$ (Chambers et al. 2016; Kennea & Krimm 2018). This suggests that the optical source will still likely be detectable after J1858 returns to quiescence, allowing for follow-up optical/near-infrared spectroscopy to determine the spectral type of the companion star. It may also be possible to constrain the binary orbital period and inclination once the source returns to quiescence.

The optical counterpart has also shown strong P-Cygni profiles in multiple emission lines, suggesting the system contains a strong optical wind (Munoz-Darias et al. 2019). The terminal velocity of the wind was found to be ~2500 km s$^{-1}$.

<footnote>$^{23}$ We note that the Pan-STARRs observations were obtained in 2012, long before the source was observed to be in outburst (Chambers et al. 2016).</footnote>
The features observed in these emission lines were also found to vary on ~5 minute timescales (Munoz-Darias et al. 2019). There is also evidence that the terminal wind velocity has changed in magnitude over time, with earlier observations having smaller wind velocities of ~500–1500 km s⁻¹ (Munoz-Darias et al. 2019).

It is interesting to note that the source was detected in optical prior to its outburst, and in UV during the outburst, suggesting that the source is likely to be relatively nearby. There appears to be a low amount of extinction (i.e., $E(B-V) \approx 0.3$ out to a distance of ~4 kpc) in this direction from the Bovy et al. (2016) extinction maps, but unfortunately the maps do not extend beyond ~4 kpc. Therefore, it is difficult to set any constraints on the source distance at this time. If, however, we assume a fiducial distance of ~5 kpc, then the unabsorbed X-ray luminosity of flare 2 (see Table 2) reaches $L_X = 1.4 \times 10^{37} (d/5 \text{ kpc})^2 \text{ erg s}^{-1}$ in the 0.6–79 keV energy band. We find an unabsorbed peak X-ray luminosity for flare 2 of $L_X \approx 7 \times 10^{37} (d/5 \text{ kpc})^2 \text{ erg s}^{-1}$ in the 0.6–79 keV energy band as the peak count rate for flare 2 is about a factor of five larger than the average flare count rate. This luminosity is about 5% of the Eddington luminosity assuming a 10 solar mass BH. Therefore, this system would need to be at a distance of $\gtrsim 20$ kpc to be accreting near the Eddington limit. Moreover, optically bright BH transients have been found to lie at a median distance of ~2 kpc (Gandhi et al. 2019), suggesting that distances $>20$ kpc are somewhat unlikely.

### 4.3. Comparison with Similar Systems

J1858’s flaring behavior in the X-ray through optical bands,<fnref n="24" /> lack of a well-defined spectral state, large intrinsic absorbing column density, and winds detected in the optical band is reminiscent of the well-known sources V404 Cyg and V4641 Sgr. V404 Cyg, similar to J1858, showed rapid and large amplitude flares in its X-ray light curves during its outbursts in 1989 and more recently in 2015 (see, e.g., Kitamoto et al. 1989; Życki et al. 1999; Walton et al. 2017). In fact, during V404 Cyg’s 2015 outburst, the hardness ratios of a majority of the flares were found to be around $R_{10-79 \text{ keV}}/R_{3-10 \text{ keV}} \approx 0.5$, similar to what is observed for J1858 (see Figure 3 in Walton et al. 2017). Additionally, V404 Cyg showed dramatic changes in the intrinsic absorbing column density, reaching values as large as $N_H \approx 10^{25} \text{ cm}^{-2}$, which were anticorrelated with the flux of the source (Motta et al. 2017a, 2017b). At optical wavelengths, slow blue flares and fast red flares, like those observed in J1858, were also exhibited by V404 Cyg (Gandhi et al. 2016). The detection of varying P-Cygni profiles in the optical spectra of J1858 also suggests that it has a high-velocity (~2500 km s⁻¹) outer accretion disk wind, similar to those observed during the 2015 outburst of V404 Cyg (Munoz-Darias et al. 2016; Munoz-Darias et al. 2019). Lastly, low-frequency QPOs have been detected in both V404 Cyg and J1858 (Huppenkothen et al. 2017; Hare et al. 2019).

V4641 Sgr had a major outburst in 1999, which showed short duration large amplitude flares similar to those observed in V404 Cyg and J1858 (Wijnands & van der Klis 2000; Revnivtsev et al. 2002a). Additionally, it was suggested that an enshrouding envelope, which caused variable absorption,

surrounded the inner accretion disk of V4641 Sgr (Revnivtsev et al. 2002b). This source has also exhibited a fast wind, with velocities up to ~3000 km s⁻¹, and also showed P-Cygni profiles in its optical spectra, suggesting that it is coming from the outer accretion disk (Lindstrøm et al. 2005; Muñoz-Darias et al. 2018).

While both V404 Cyg and V4641 Sgr are similar to J1858 in many ways, there still remain two distinct differences between the first two systems and J1858. The first difference is that the orbital inclination angles for both V404 Cyg and V4641 Sgr are relatively large, ~67° and ~72°, respectively (Khargharia et al. 2010; MacDonald et al. 2014), while the $3\sigma$ upper limit on the disk inclination angle for J1858 derived from our reflection fits is relatively low, $i < 29°$. These large inclinations in the first two systems imply that the system is being viewed close to edge on, suggesting that the large amounts of variable obscuring material can be explained by a flared disk (possibly with clumps of material) which intersects the LOS between the observer and the inner accretion disk (see, e.g., Figure 8 in Motta et al. 2017b). Given that the inclination of the inner accretion disk of J1858 appears to be low, this explanation appears less likely to be applicable to J1858. However, misalignments of ~15° between the orbit of the system and the inner accretion disk have been observed in Cygnus X-1 (see, e.g., Tomsick et al. 2014; Walton et al. 2016) and even larger misalignments (possibly ~30°–50°) have been observed in the two systems V4641 Sgr and V404 Cyg (Maccarone 2002; Miller-Jones et al. 2019), which are very similar to J1858. Thus, the possibility of a large misalignment between the inclination of the orbit and inner accretion disk, leading to the obscuration of the inner regions of the accretion disk by a flared disk, cannot be entirely ruled out.

An alternative possibility is that there may be some systematic effects on the derived inclination of the inner accretion disk if there is a complicated source geometry (e.g., a thick disk; Taylor & Reynolds 2018) near the inner accretion flow, which is currently unaccounted for in the simplified RELXILL reflection models. For instance, using a reflection model similar to the one used here (i.e., *relxill_lp+xillver*), Walton et al. (2017) found a range of inclinations for V404 Cyg, spanning $i = 27°–52°$. Furthermore, Connors et al. (2019) found an inclination of the inner accretion disk (~40°) that also largely differed from the well-determined binary inclination (~75°) in XTE J1550–564. We also note that lamp-post geometry is an idealized, point-source geometry, so if the corona is vertically or horizontally extended it may also impact the inferred inclination. Additionally, if the illuminating X-ray source is associated with the base of the jet, it may be moving with a mildly relativistic velocity, which could also affect the inferred inclination (see, e.g., Section 5.3 in Connors et al. 2019). In any case, since the binary inclination for J1858 still remains unknown, the disk inclination cannot be compared to the binary inclination to look for a possible misalignment or discrepancy. However, since it appears that the source was detected in optical by Pan-STARRs prior to its outburst, we reiterate that it may be possible to constrain the binary inclination once the source returns to quiescence.

The second way in which J1858 differs from V404 Cyg and V4641 Sgr is that it appears to be accreting at only a few percent of the Eddington luminosity, whereas V404 Cyg and V4641 Sgr were accreting at, or possibly even above, the Eddington limit during their outbursts (Revnivtsev et al. 2002a;

<footnote n="24">²⁴ We note that rapid X-ray flaring has also been observed in several other BHBs, such as GRS 1915+105 and IGR J1709-3624 (see, e.g., Kimura et al. 2016).</footnote>
Motta et al. 2017b). The distance is not well known for J1858, so once better constraints are placed on its distance a clearer picture of its accretion rate will emerge. However, even if the distance to J1858 is as large as 10 kpc, the source luminosity would still only be ~10%-20% of the Eddington luminosity. One other possible explanation for this difference in luminosities is that for large enough absorbing column densities (i.e., $N_H \gtrsim 10^{24} \text{ cm}^{-2}$), scattering processes become important. For example, during one of V404 Cyg's plateau states, Motta et al. (2017a) found a luminosity of only a few percent of the Eddington luminosity when using a reflection model having a large best-fit absorbing column density (i.e., $N_H \approx (1-3) \times 10^{24} \text{ cm}^{-2}$), but after using a model that also accounted for scattering, they found an increase in intrinsic flux by a factor of ~30, pushing the true luminosity to the Eddington luminosity. This explanation is also somewhat unlikely though, considering that during flaring episodes (e.g., during flares 1 and 2), J1858's intrinsic absorption is about a factor of 20 lower than in V404 Cyg during this plateau state.

## 5. Summary

We have reported on the first *NuSTAR* observation of Swift J1858.6–0814, which is a newly discovered BH binary candidate. The main findings of this work can be summarized as follows:

1. The source exhibits large amplitude flares, showing an increase in count rate by a factor of ~100 on timescales of 10–100 s, but we find no evidence of any periodicity or QPOs in the *NuSTAR* light curves.

2. The observed flares are accompanied by large changes in the source's hardness ratio, suggesting that the source's spectrum also significantly changes during the flares. We split the source's energy spectra into three different modes based on where the source is located in the HID during a given interval. We then fit these spectra with a relativistic reflection model, allowing us to constrain a number of the source's physical parameters.

3. Interestingly, we find a large and variable partially covering absorbing column density ($N_H = (14-42) \times 10^{22} \text{ cm}^{-2}$) dependent on the mode of the source. We also find that a thermal component is required to adequately fit the spectra. These results are supported by the joint fits to the energy spectra from two flares simultaneously observed by *NuSTAR* and *NICER*.

4. We constrain the BH spin to be $a^* > 0.0$ at the 90% level and $a^* > -0.8$ at the 3σ level, assuming the inner radius of the accretion disk is at the ISCO.

5. The inclination of the inner accretion disk derived from our fits appears to be relatively low ($i < 29^\circ$ 3σ upper limit), making the origin of the large amount of obscuring material unclear. Future comparisons of the derived inner accretion disk inclination with the binary inclination angle (once it is known) can help further our understanding of this system.

6. The source shows many similarities to the well-known Galactic BH binaries V404 Cyg and V4641 Sgr. However, J1858's low inclination and low luminosity differ greatly from these two sources, making a direct comparison difficult at this point in time.

This source was followed up by a large multiwavelength campaign, including five additional *NuSTAR* observations.

Therefore, our understanding of this interesting source will continue to grow as more results are released.

This work made use of data from the *NuSTAR* mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by the National Aeronautics and Space Administration. We thank the *NuSTAR* Operations, Software and Calibration teams for support with the execution and analysis of these observations. This research has made use of the *NuSTAR* Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (USA). J.H. and J.A.T. acknowledge partial support from *NuSTAR* Guest Observer grant 80NSSC19K0404. J.H. acknowledges support from an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center, administered by the USRA through a contract with NASA. D.J. W. acknowledges support from STFC in the form of an Ernest Rutherford Fellowship. M.C. acknowledges support from the Centre National d'Etudes Spatiales (CNES). J.A.G. acknowledges support from NASA grant 80NSSC19K1020 and from the Alexander von Humboldt Foundation. We thank the anonymous referee for providing useful comments which have improved the quality of this paper.

**Software:** XSPEC (v12.10.1; Arnaud 1996), NUSTARDAS (v1.8.0), NICERDAS (V005), Stingray (Huppenkothen et al. 2017), Xselect (v2.4e), Matplotlib (Hunter 2007), HEASOFT (v6.25), MWDust (Bovy et al. 2016), XILLVER (García & Kallman 2010; García et al. 2013), RELXILL (v1.2.0; Dauser et al. 2014; García et al. 2014).

## ORCID iDs

Jeremy Hare https://orcid.org/0000-0002-8548-482X

John A. Tomsick https://orcid.org/0000-0001-5506-9855

Maïca Clavel https://orcid.org/0000-0003-0724-2742

Poshak Gandhi https://orcid.org/0000-0003-3105-2615

Javier A. García https://orcid.org/0000-0003-3828-2448

Brian W. Grefenstette https://orcid.org/0000-0002-1984-2932

Dominic J. Walton https://orcid.org/0000-0001-5819-3552

Yanjun Xu https://orcid.org/0000-0003-2443-3698

## References

Arnaud, K. A. 1996, in ASP Conf. Ser 101, Astronomical Data Analysis Software and Systems V (San Francisco, CA: ASP), 17

Bachetti, M., Harrison, F. A., Cook, R., et al. 2015, ApJ, 800, 109

Baglio, M. C., Russell, D. M., Pirbhoy, S., et al. 2018, ATel, 12180

Belloni, T. M., & Motta, S. E. 2016, in Astrophysics of Black Holes: From Fundamental Aspects to Latest Developments, ed. C. Bambi (Berlin: Springer), 61

Beri, A., Tetarenko, B. E., Bahramian, A., et al. 2019, MNRAS, 485, 3064

Bhalerao, V., Romano, P., Tomsick, J., et al. 2015, MNRAS, 447, 2274

Bovy, J., Rix, H.-W., Green, G. M., et al. 2016, ApJ, 818, 130

Bright, J., Fender, R., Motta, S., et al. 2018, ATel, 12184, 1

Buisson, D. J. K., Fabian, A. C., Barret, D., et al. 2019, MNRAS, 490, 1350

Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv:1612.05560

Connors, R. T., García, J., Steiner, J., et al. 2019, ApJ, 882, 179

Dauser, T., García, J., Parker, M. L., Fabian, A. C., & Wilms, J. 2014, MNRAS, 444, L100

Dauser, T., García, J., Walton, D. J., et al. 2016, A&A, 590, A76

Dauser, T., Wilms, J., Reynolds, C. S., et al. 2010, MNRAS, 409, 1534

Fabian, A. C., Rees, M. J., Stella, L., et al. 1989, MNRAS, 238, 729

Ferrigno, C., Bozzo, E., Sanna, A., et al. 2019, A&A, 624, A142

Gandhi, P., Bachetti, M., Dhillon, V. S., et al. 2017, NatAs, 1, 859

Gandhi, P., Littlefair, S. P., Hardy, L. K., et al. 2016, MNRAS, 459, 554
Gandhi, P., Rao, A., Johnson, M. A. C., et al. 2019, MNRAS, 485, 2642

García, J., Dauser, T., Lohfink, A., et al. 2014, ApJ, 782, 76

García, J., Dauser, T., Reynolds, C. S., et al. 2013, ApJ, 768, 146

García, J., & Kallman, T. R. 2010, ApJ, 718, 695

Gendreau, K. C., Arzoumanian, Z., & Okajima, T. 2012, Proc. SPIE, 8443, 844313

Hare, J., Gandhi, P., Paice, J. A., et al. 2019, ATel, 12512

Hare, J., Halpern, J. P., Clavel, M., et al. 2019, ApJ, 878, 15

Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103

HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116

Homan, J., Steiner, J. F., Lin, D., et al. 2018, ApJ, 853, 157

Hunter, J. D. 2007, CSE, 9, 90

Huppenkothen, D., Bachetti, M., Stevens, A. L., et al. 2019, ApJ, 881, 39

Huppenkothen, D., Younes, G., Ingram, A., et al. 2017, ApJ, 834, 90

Jaisawal, G. K., Naik, S., & Chenevez, J. 2018, MNRAS, 474, 4432

Kennea, J. A., & Krimm, H. A. 2018, ATel, 12160

Khargharia, J., Froning, C. S., & Robinson, E. L. 2010, ApJ, 716, 1105

Kimura, M., Isogai, K., Kato, T., et al. 2016, Natur, 529, 54

Kitamoto, S., Tsunemi, H., Miyamoto, S., et al. 1989, Natur, 342, 518

Krimm, H. A., Barthelmy, S. D., Cummings, J. R., et al. 2018, ATel, 12151

Krimm, H. A., Holland, S. T., Corbet, R. H. D., et al. 2013, ApJS, 209, 14

Laor, A. 1991, ApJ, 376, 90

Lindström, C., Griffin, J., Kiss, L. L., et al. 2005, MNRAS, 363, 882

Ludlam, R. M., Miller, J. M., Arzoumanian, Z., et al. 2018, ATel, 12158

Maccarone, T. J. 2002, MNRAS, 336, 1371

MacDonald, R. K. D., Bailyn, C. D., Buxton, M., et al. 2014, ApJ, 784, 2

Miller-Jones, J. C. A., Tetarenko, A. J., Sivakoff, G. R., et al. 2019, Natur, 569, 374

Mitsuda, K., Inoue, H., Koyama, K., et al. 1984, PASJ, 36, 741

Mori, K., Gotthelf, E. V., Zhang, S., et al. 2013, ApJL, 770, L23

Motta, S. E., Kajava, J. J. E., Sánchez-Fernández, C., et al. 2017a, MNRAS, 468, 981

Motta, S. E., Kajava, J. J. E., Sánchez-Fernández, C., et al. 2017b, MNRAS, 471, 1797

Muñoz-Darias, T., Casares, J., Mata Sánchez, D., et al. 2016, Natur, 534, 75

Munoz-Darias, T., Jimenez-Ibarra, F., Armas Padilla, M., et al. 2019, ATel, 12881

Muñoz-Darias, T., Torres, M. A. P., & Garcia, M. R. 2018, MNRAS, 479, 3987

Paice, J. A., Gandhi, P., Dhillon, V. S., et al. 2018, ATel, 12197

Rajwade, K., Kennedy, M., Breton, R., et al. 2018, ATel, 12186

Reis, R. C., Fabian, A. C., & Miller, J. M. 2010, MNRAS, 402, 836

Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49

Revnitsev, M., Gilfanov, M., Churazov, E., et al. 2002a, A&A, 391, 1013

Revnitsev, M., Sunyaev, R., Gilfanov, M., et al. 2002b, A&A, 385, 904

Reynolds, M. T., & Miller, J. M. 2013, ApJ, 769, 16

Reynolds, M. T., Miller, J. M., Ludlam, R. M., et al. 2018, ATel, 12220

Taylor, C., & Reynolds, C. S. 2018, ApJ, 855, 120

Tomsick, J. A., Nowak, M. A., Parker, M., et al. 2014, ApJ, 780, 78

Vasilopoulos, G., Bailyn, C., & Milburn, J. 2018, ATel, 12164

Walton, D. J., Mooley, K., King, A. L., et al. 2017, ApJ, 839, 110

Walton, D. J., Tomsick, J. A., Madsen, K. K., et al. 2016, ApJ, 826, 87

Wijnands, R., & van der Klis, M. 2000, ApJL, 528, L93

Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914

Xu, Y., García, J. A., Fürst, F., et al. 2017, ApJ, 851, 103

Xu, Y., Harrison, F. A., García, J. A., et al. 2018a, ApJL, 852, L34

Xu, Y., Harrison, F. A., Kennea, J. A., et al. 2018b, ApJ, 865, 18

Zdziarski, A. A., Johnson, W. N., & Magdziarz, P. 1996, MNRAS, 283, 193

ZYcki, P. T., Done, C., & Smith, D. A. 1999, MNRAS, 309, 561