Chandra reveals a possible ultra-fast outflow in the super-Eddington Be/X-ray binary Swift J0243.6+6124

Accretion at super-Eddington rates is expected to be accompanied by strong outflows. Such outflows are observed in Galactic X-ray binaries and extragalactic Ultra-luminous X-ray sources (ULXs). However, due to their large source distances, ULX outflows are challenging to detect and study in detail. Galactic neutron stars accreting from a Be-star companion at super-Eddington rates show many similarities to ULX pulsars, and therefore offer an alternative approach to study outflows in this accretion regime. Here, we present Chandra high-resolution spectroscopy of such a super-Eddington accreting neutron star, Swift J0243.6+6124, to search for wind outflow signatures during the peak of its 2017/2018 giant outburst. We detect narrow emission features at rest from Ne, Mg, S, Si, and Fe. In addition, we detect a collection of absorption features which can be identified in two ways: either as all Fe transitions at rest (with a possible contribution from Mg), or a combination of three blue-shifted Ne and Mg lines at $\sim 0.22c$, while the remaining lines are at rest. The second scenario would imply an outflow with a velocity similar to those seen in ULXs, including the ULX pulsar NGC 300 ULX-1. This result would also imply that Swift J0243.6+6124 launches both a jet, detected in radio and reported previously, and an ultra-fast wind outflow simultaneously at super-Eddington accretion rates.


INTRODUCTION
The accretion and subsequent ejection of matter is a ubiquitous process in the Universe, occurring in objects ranging from young stellar objects and planet-forming systems (e.g.Kuiper et al. 2015;Beltrán & de Wit 2016) to X-ray binaries (see e.g.Fender et al. 2004;Migliari & Fender 2006, for overviews) and active galactic nuclei (e.g.Merloni et al. 2003;Falcke et al. 2004).In X-ray binary systems (XRBs) a stellar-mass compact object, either a black hole or a neutron star, accretes from an companion star in a close orbit.X-ray binaries come in numerous classes, with different combina-E-mail: a.j.vandeneijnden@uva.nl(JvdE) tions of compact object type, and mass and type of companion.Additionally, the accretion can take place through various channels (Frank et al. 1992), triggered by for instance Roche-lobe overflow of the donor (Kuiper 1941;Paczyński 1971), a stellar wind (e.g.Reig 2011) or the movement of the compact object through the circumstellar disk of the donor star (Okazaki & Negueruela 2001).
Accretion in X-ray binaries is often accompanied by the ejection of matter, either through disk winds or via jets.The latter are strongly-collimated outflows traveling near the speed of light, launched from the inner accretion flow, while winds are launched further out from the accreting object at lower velocities (ranging from hundreds of km/s to ∼ 0.3c) and with wider opening angles.In X-ray binaries, winds can carry away a large fraction of the mass from the accretion flow (Neilsen & Lee 2009;Ponti et al. 2012), possibly triggering instabilities in the flow (Begelman et al. 1983;Muñoz-Darias et al. 2016) and potentially affecting the outburst profiles of transient sources (Tetarenko et al. 2018).Similarly, jets can remove accretion power from the X-ray binary and deposit large amounts of energy in the surrounding interstellar medium (Fabrika 2004;Fender et al. 2005;Gallo et al. 2005;Pakull et al. 2010).
Outflows can alternatively be studied in X-ray binaries accreting around or above the Eddington luminosity; in such sources, strong outflows are expected due to the enhanced radiation pressure exerted on the accretion flow (Shakura & Sunyaev 1973;Ohsuga & Mineshige 2011;McKinney et al. 2014McKinney et al. , 2015;;Hashizume et al. 2015;King & Muldrew 2016).Famous examples of X-ray binaries launching strong outflows in this regime are the black-hole LMXBs GRS 1915+105 (Mirabel & Rodriguez 1994;Neilsen & Lee 2009) and V404 Cygni (Muñoz-Darias et al. 2016;Tetarenko et al. 2017), and the accreting neutron star Cir X-1 (Brandt & Schulz 2000).During super-Eddington X-ray binary states, the apparent dichotomy between winds and jets can also break down; for instance, black holes in XRBs and Z-sources -a subset of neutron star LMXBs categorized based on their X-ray color-color diagram (Hasinger & van der Klis 1989) and accreting around the Eddington limit -are thought to launch a wind and jet simultaneously at such accretion rates (Homan et al. 2016;Allen et al. 2018).
In this regard, Ultra-luminous X-ray sources (ULXs) are particularly interesting (see Kaaret et al. 2017, for a recent review).These extragalactic X-ray emitters have X-ray luminosities (greatly) exceeding ∼ 10 39 erg/s, or the Eddington luminosity of a 10M black hole.Recently, a handful of ULXs has been identified as accreting neutron stars through the detection of pulsations (Bachetti et al. 2014;Fürst et al. 2016;Israel et al. 2017a,b;Carpano et al. 2018), with several additional candidates found through possible cyclotron resonance scattering features (Brightman et al. 2018;Koliopanos et al. 2019;Walton et al. 2018c).This confirms the super-Eddington nature of at least a fraction of ULXs.While it is unclear what fraction of ULXs contains a pulsar, both theoretical (King et al. 2001;King & Lasota 2016) and observational studies (Koliopanos et al. 2017;Walton et al. 2018a) suggest it could be substantial.
Theoretically, outflows have often been suggested to explain the soft spectra of (some) ULXs (King 2001;Begelman 2002;King & Pounds 2003;Gladstone et al. 2009;Feng & Soria 2011;Urquhart & Soria 2016).From the observational side, jets have been observed indirectly from these sources through their impact on the surrounding medium (Middleton et al. 2013;Cseh et al. 2014Cseh et al. , 2015a,b;,b;Mezcua et al. 2015) and inferred from radio detections (e.g.Kaaret et al. 2003;Webb et al. 2012;Mezcua et al. 2013).Winds have been seen through X-ray (Pinto et al. 2016) and optical (Zepf et al. 2008) spectroscopy of several targets, including one ULX pulsar (NGC 300 ULX-1; Kosec et al. 2018a).However, the extragalactic nature of ULXs complicates the study of their outflows.Scaling for instance typical radio luminosities of compact jets launched by black holes accreting around the Eddington luminosity (e.g.Gallo et al. 2018) to Mpc distances, yields flux densities at best around the detection limit for current generation radio arrays.Similarly, the detection of winds through X-ray spectroscopy is limited by the low number of counts in the X-ray grating spectra.As a result, most high-resolution X-ray spectra of ULXs currently available in the archive are not sensitive enough to reveal any wind signatures (Kosec et al. 2018b).
With their smaller distances, X-ray binaries in the Milky Way or the Small Magellanic Cloud may offer a valuable alternative avenue to study these super-Eddington accretion states at higher signal-to-noise ratio (although the sample of such sources is limited by the smaller volume and extreme count rates can introduce calibration issues, as discussed later).In particular, neutron star Be/X-ray binaries, wherein the donor is a Be-star, show many similarities to the known ULX pulsars (Mushtukov et al. 2015;Koliopanos et al. 2017): strong (≥ 10 12 G) magnetic fields and slow spins (i.e.periods on the order of seconds).Importantly, Be/X-ray binaries can also show super-Eddington accretion rates during the peaks of their giant outbursts (e.g Reig 2011).
In September 2017, the Swift satellite discovered the new neutron star Be/X-ray binary Swift J0243.6+6124(Kennea et al. 2017, hereafter Sw J0243;): a strongly-magnetized neutron star (Tsygankov et al. 2018, e.g.B > 10 12 G;) with a ∼ 9.8 second spin period (Kennea et al. 2017).It reached super-Eddington X-ray luminosities during the peak of its outburst (van den Eijnden et al. 2018b) and has been referred to as the first Galactic ULX pulsar by both Tsygankov et al. (2018) and Wilson-Hodge et al. (2018).Through Very Large Array (VLA) radio and Swift X-ray monitoring, Sw J0243 was found to launch a relativistic jet during its super-Eddington state (van den Eijnden et al. 2018b), as well as at lower accreting rates (van den Eijnden et al. 2019).This jet detection constituted the first from a strongly-magnetized neutron star, contrary to the predictions of neutron star jet formation theory (Blandford & Payne 1982;Massi & Kaufman Bernadó 2008).
Here, we present Chandra high-resolution gratings Xray spectroscopy of Sw J0243 in its super-Eddington state.We find evidence for a wind with a velocity of ∼ 0.22c through the detection of blue-shifted absorption features (Section 4.3.2),similar to those detected in ULXs (Kosec et al. 2018a).If these features indeed arise from a wind, this would imply both a jet and a wind are launched simultaneously by Sw J0243 during its super-Eddington state.

OBSERVATIONS AND DATA REDUCTION
In Figure 1, we show the X-ray and optical light curves of the 2017/2018 outburst of Sw J0243.Chandra observed the target around the peak of this giant outburst.While the distance to Sw J0243 is not precisely known, the Gaia DR2 implies a minimum of 5 kpc at 99% confidence (van den Eijnden et al. 2018b).Given this minimum distance, the Xray luminosity during the stage of the outburst around the Chandra epoch exceeded 10 39 erg/s in the 0.5 − 10 keV band (Wilson-Hodge et al. 2018;van den Eijnden et al. 2018b).As the theoretical Eddington luminosity for a neutron star is 2× 10 38 erg/s, this luminosity implies a firmly super-Eddington accreting rate.
Chandra performed Director's Discretionary Observations of Sw J0243 on 11 November 2017 (MJD 58068) for ∼ 25 ks of exposure (ObsID 20859; PI Degenaar) with the high energy transmission grating spectrometer (HETGS).Given the extreme flux of the source, the observation setup had to mitigate photon pile-up as well as minimize telemetry saturation.The telescope aimpoint was moved to the CCD readout where the zero order dithered between the framestore and a few CCD rows and thus only partially covered the active CCD area.In this configuration only two grating dispersion arms are recorded, the medium energy grating (MEG) +1st order and the high energy grating (HEG) -1st order plus their higher order dispersions.The observation was recorded in continuous clocking mode (CC-mode) with a fast readout time of 3.85 msec to effectively mitigate photon pile-up in the dispersed spectra.The data were transmitted via GRADED mode which includes onboard event grading and grade summing.This results in some loss of data information and aspects of calibration become more approximate.While this mostly preserves the detection of discrete line features and edges, it does affect the calibration of the spectral continua in the first order spectra see also Schulz et al. 2009 for another example of CC-mode observations where mainly the discrete line features are preserved, and Miller et al. 2003 for the analysis of CC-mode spectra of the accreting black hole XTE J1550-564).
The observation data were re-processed via CIAO 4.9 using the latest calibration product at the time (CALDBv4.7.7).Changes in more recent versions did not have any impact on the analysis at the time of submission.We used the run pipe thread within the tgcat package in ISIS1 .Under normal circumstances this determines the wavelength scale to the accuracy of a quarter of a resolution element, i.e. 0.005 Å for MEG and 0.002 Å for HEG.However due to the fact that the zero order is piled as well as that it dithers on and off the chip likely adds another quarter in systematic uncertainty.In gratings the dispersion scale is linear in wavelength and all line analysis will be done in wavelength space.We used standard wavelength redistribution matrices (RMFs) and generated ancillerary response files (ARFs) applying provided aspect solutions, bad pixel maps, and CCD window filters2 .
Given the possible continuum issues due to the extreme count rates, we also extracted two Swift XRT spectra to compare with the Chandra HEG and MEG spectra.We used the Swift XRT data products generator (Evans et al. 2007) 3 to extract the Swift spectra taken on 10 and 13 November 2017 (MJDs 58067 and 58070, with ObsIDs 10336023 and 10336024, respectively).Both observations were taken in WT-mode, which can deal with high count rates, and the data products generator automatically corrects for any pile-up issues for very bright sources.

Continuum analysis and spline modeling
In Figure 2 (top), we show both the Chandra MEG and HEG spectra with the two Swift XRT spectra taken before and after the Chandra epoch.As described in Section 2, the unconventional HETGS observing setup required by the extreme flux can affect the calibration of the continuum in the first order spectra.This is clearly at play in our observations of Sw J0243: inaccuracies in the HETGS continua are obvious from the disparity between the MEG and HEG detectors, the large jump in the HEG spectrum around ∼ 7 Å (∼ 1.77 keV), and both the offset and difference in spectral shape between the Swift and Chandra spectra.The latter is highlighted by the bottom panel of Figure 2, where we show the ratio of the MEG and HEG data to a simple tbabs*(bbodyrad+po) model fit jointly to the XRT spectra with xspec (Arnaud 1996).
However, while the calibration of the continuum shape is affected by the observational setup, discrete line features and edges are preserved (Schulz et al. 2009).Indeed, individual narrow features remain.In the HEG spectrum, for instance, a clear Fe K complex around ∼ 6.4-7 keV is visible, as discussed in detail in Section 4.2.1.Therefore, one can still search for narrow emission and absorption features from, for instance, outflows or donor star material.As the offset between the Chandra and Swift spectra demonstrates, the HEG and MEG data will not provide an accurate measurement of the flux.As a result, conventional line strength measures such as the equivalent width or normalization of any such narrow features will not be accurate.Furthermore, without an accurate continuum measurement, physically motivated modelling is challenging for any model that contains a (pseudo)-continuum component.But despite these restrictions, the presence of narrow emission and absorption features can still be tested and give valuable physical information about the system.
Any search algorithm for individual narrow X-ray spectral lines requires an accurate description of the continuum.This holds especially for the approach that we adopt for Sw J0243, which was originally developed for the detection of faint features in ULX spectra by Pinto et al. (2016) and is introduced in Section 3.2 (for an illustration of the effects of a poorly modeled continuum on the detection and significances of narrow line features, see also van den Eijnden et al. 2018a).However, this line search approach only requires an accurate shape of the continuum model to find narrow deviations from; this continuum model does not necessarily have to be physically motivated.A example of this can be found in Grinberg et al. (2017), where a Chandra spectrum of the HMXB Vela X-1 is modeled; in that work, the continuum consists of four independent power law models, which are not physically motivated, fitted over a limited wavelength range.However, these models do provide an accurate description of the underlying continuum shape and allow for the search and identification of narrow line features.A similar mathematical approach to model the continuum can be found in Yao et al. (2009), where a combination of broad Gaussians makes up the continuum model.Therefore, instead of fitting the Chandra Sw J0243 continuum with physical modelswhich is not possible for the full spectral range, as shown in Figure 2 -we apply a spline interpolation as the continuum instead.
To calculate the spline interpolations of the HEG and MEG spectra, we first used xspec to write out the flux as a function of wavelength.We then choose the step size of the interpolation -as we aim to search for deviations from the spline continuum, we should not interpolate every spectral bin but instead bins separated by a fixed wavelength -and defined a wavelength grid with such steps on the considered wavelength range (note that, therefore, this grid is not the same as wavelength bins of the HEG and MEG detectors).Simply calculating a spline between the fluxes on this grid has the risk of accidentally using either a statistical outlier or a spectral bin inside a narrow line feature as part of the continuum.To prevent this effect, we instead calculated a first continuum estimate with the third degree spline between the fluxes F i at each grid point in wavelength λ i .Then, to obtain the final spline continuum model, we fitted the spline function to the entire spectrum with the F i values as free parameters, recalculating the spline for each updated combination of F i and minimizing the χ 2 value between spline and data.Finally, the resulting best-fit continuum spline model was saved as an additive xspec fits table model.
We calculated separate spline continuum models for the MEG (2.07-13.78Å) and HEG spectra.We used two individual splines for the HEG data, covering 1.55-7 Å and 7.1-12.4Å, in order to account for the steep jump at 7.06 Å.As the presence of a narrow feature in both the HEG and MEG data is important to conclude it is not a mere statistical fluctuation, this implies we exclude the 7-7.1 Å range from our entire analysis.We tried different combinations of step sizes (0. 5, 1, 1.5 and 2 Å) and spectral  rebinning, and rebinning to a S/N of 10 and 50 per spectral bin) for the calculation of the continuum splines.After a combination of visual inspection and comparison of the continuum χ 2 values, we concluded that a 0.5 Å stepsize and no rebinning provided the most accurate continuum for both the HEG and the MEG data: this combination systematically resulted in the lowest χ 2 values (of the order of χ 2 ν ≈ 1.2), while the other combinations (especially with stepsize ≥ 1 Å) introduced significant residual structure between the gridpoints interpolated by the spline.

Line search
We adopt the line search algorithm developed by Pinto et al. (2016) and refer the reader to that paper and to Kosec et al. (2018a), Kosec et al. (2018b), andvan den Eijnden et al. (2018a) for an extensive description of its details.The basic rationale is as follows: after setting a continuum model -the spline interpolations in the case of Sw J0243 -we define a grid in wavelength and choose a fixed velocity line width.We then step through the wavelength grid, fitting a single Gaussian function with the fixed velocity width and a free normalization, centered at the grid point.The significance of an emission (i.e.positive normalization) or absorption (negative normalization) line at that wavelength is then recorded as the fitted normalization divided by its one-sigma error.Alternatively, the line significance can also be probed by the improvement in fit statistic (either χ 2 or C-statistic).
For the correct interpretation of the results of this line search method, several caveats are important to keep in mind.Firstly, despite fitting the spline continuum to cancel the effect of outliers, this continuum does not necessarily describe the entire spectrum accurately.In extreme cases, such as the 7 Å jump in the HEG spectrum, this requires the calculation of multiple splines.However, for less extreme cases, it also implies that care should be taken when considering the physical origin of any suggested lines, and one should carefully inspect the spectra and splines around the possible line features.Furthermore, the returned significances are single-trial estimates, while estimating the number of independent trials is not straightforward.Therefore, line search results of a single spectrum should be treated with caution.
To reflect this, we do not quote the single-trial significances of any detected lines as actual significances.
For Sw J0243, we have two simultaneous but independent spectra from two different detectors with different instrument responses.This latter point is important, as any imperfections in the response modeling might appear as deviations from the continuum and resemble a narrow spectral line.However, such features are not expected to appear in both spectra, unless the same response feature is present in both detectors.
To take these caveats into accounts, we require that any possible spectral lines possess the following properties before considering them as real spectral features of the X-ray binary: (i) the line should be 3σ significant (single trial) in both the HEG and the MEG spectrum, (ii) the line should not be located on top of a shared response feature of both detectors, (iii) the continuum model should look accurate around the central wavelength of the line in both detectors and the presence of a spectral feature should hold up to visual inspection of the spectrum, and (iv) the centroid energies, where the significance peaks, of the line in the two detectors should be close: to account for slight statistical deviations between the peak wavelengths and possible small inaccuries in calibration, we require those centroids to be within ∼ 0.01 keV.Any combination of spectral features adhering to these requirements should of course also fit within a consistent, physically realistic picture of the X-ray binary system and its state.In addition, we also performed a careful visual inspection of the spectra and the line-search results to identify possible features in low S/N parts of the spectra, where lines are less likely to be picked up as significant by the line search.

Robustness of the spline continuum
As we did not model the continuum with a physical model, but instead with a spline interpolation, we performed several checks of our approach; specifically, we tested whether the line search results, and our inferences, were directly affected by the choice of continuum.For this purpose, we designed two tests: comparing the interpolated continuum with a physical continuum, and slightly varying the step size of the spline grid points.
Figure 2 shows that the Chandra and Swift spectral shapes do not generally match.However, above ∼ 1.8 keV, the (blue) HEG spectrum and both Swift XRT spectra appear to have a similar shape.Therefore, we fitted two continuum models to these three spectra, using the HEG data between ∼ 1.77 and ∼ 8 keV (7.0 -1.55 Å), and both XRT spectra between 1 and 10 keV.Using xspec, and assuming abundances from Wilms et al. (2000) and cross-sections from Verner et al. (1996), we fitted both a tbabs*po and a tbabs*(bbodyrad+po) model as simple phenomenological continuum models.In both cases, we included a multiplicative constant to account for offsets between the spectra, while keeping all other parameters tied.Using both these continuum models, we re-applied our line search pipeline with a 500 and 2000 km/s line width.This re-analysis finds the same narrow features in the line search, although some residual trends in the line significances remain when using the physical continuum models.These trends suggest that while the HEG and XRT spectra appear similar above 1.8 keV, small deviations in shape are present.We show these results in more detail in Figure A1 in Appendix A.
Secondly, we re-performed our analysis using a slightly smaller step size for the calculation of the spline continuum -0.48 Å instead of 0.5 Åtherefore smoothly connecting different spectral bins with the spline.This check should therefore reveal any imperfections due to the spline by chance connecting statistical outliers and/or narrow lines, instead of probing the continuum.As shown in Figure A2 in Appendix A in more detail, the line search results of this test are consistent with our first analysis and do not imply changes in the detected narrow features.

Line search results
In Figure 3, we show an overview of the analysis and the results of the line search.In the top panel, we show the HEG and MEG spectra.The large jump in the HEG spectrum around 7 Å signals the complications in measuring the continuum shape discussed extensively in the previous section.The solid, black lines show the spline continuum models used in the line search.To deal with the 7 Å jump, two different splines are used for the HEG spectrum, while the range between 7 and 7.1 Å is removed from the analysis.The top panel also shows the effective area shape for both detectors in arbitrary units, to indicate the instrument response.This can be used to test whether any narrow features identified in the line search might be instrumental, and shows that the 7 Å jump is directly on top of the strongest HEG detector feature.We note that using the gain command in xspec to investigate the response feature did not provide a simple gain shift solution to reduce the large response residuals at for instance ∼ 7 Å and ∼ 6.75 Å.
The middle panel shows the ratio between the spectra and the spline continuum models, providing a visual aid in searching for and confirming the physical nature of any narrow features.The results of the line search are shown in the bottom panel: we plot the single-trial significance (N/σ N ) of a Gaussian line of fixed width, added at the given energy.The red and black lines shows the 2000 km/s velocity width search, while the black area shows the results for 500 km/s.The 3σ and 5σ single-trial significance thresholds are shown to guide the eye.We re-emphasize that we analyzed the HEG and MEG spectra separately, to obtain independent search results that can be compared to distinguish physical lines from statistical fluctuations or instrumental features.Note that a negative significance signals an absorption feature.
Finally, in all panels, the grey dotted lines indicate the narrow lines identified following the requirements set out in Section 3.2.Note that a clear Fe K complex is visible in the HEG spectrum below 2 Å (∼ 6.4-7 keV), which we do not indicate with grey lines for clarity of the figure.While this region is only covered in the HEG spectrum, the shape and centroid energies -matching the Fe fluorescence lines expected in BeXRBs (Torrejón et al. 2010, see the next section) -of the three lines clearly show that this feature is real.All identified lines are listed in Tables 1 (emission), 2, and 3 (both absorption).

Significance simulations
An alternative significance estimator for narrow features is the change in fit statistic after the addition of a narrow line at a certain wavelength, ∆C(λ).This estimator offers the options of combining the results from two independently analysed spectra, by linearly adding them as ∆C(λ) = ∆C MEG (λ) + ∆C HEG (λ).We plot the combined ∆C values from the MEG and HEG detectors, for the 500 and 2000 km/s velocity width separately, in  B2).
To assess the significance of narrow features beyond the single-trial estimates shown in Figure 3, we perform Monte-Carlo simulations of the continuum spline models.For each combination of detector (HEG or MEG) and velocity width (500 or 2000 km/s), we use XSPEC to simulate 3000 spectra based on the continuum spline model.For each simulated spectrum, we then repeat the line search, with a reduced resolution of 0.05 Å to optimize computational time.At each trial wavelength, we calculate the 2 and 3σ confidence levels by calculating the 95.4 th and 99.7 th percentile of the simulated ∆C values, respectively.For the combined HEG and MEG results, we first add the simulated fit improvements, and calculate the above percentiles for this summed ∆C.The simulated confidence levels are shown in Figures 4, B1, and B2, as the blue dotted and dash-dotted lines, respectively.
The majority of lines identified in Figure 4 and listed in Tables 1, 2, and 3, shown by the grey dotted lines, show up as ≥ 3σ significant in either the 500 or 2000 km/s line search (or in both).The exception is the 12.15 Å emission line, which is only ∼ 2σ significant at 500 km/s.Unsurprisingly, this line is also not formally significant in Figure 4, which can be explained by the poor S/N in this region of the spectra.However, given the prominence of this feature in other BeXRBs (see Section 4.2.1 and 5), we include it in our further analysis.Finally, several other features appear significant above 3σ.However, these all originate from a strong feature in only one of the two detectors and we therefore do not further analyse these (see also Appendix B).

Line identification
High-resolution X-ray observations of both super-Eddington X-ray binaries (Pinto et al. 2016;Koliopanos & Vasilopoulos 2018) and BeXRBs at lower accretion rates (La Palombara et al. 2016;Grinberg et al. 2017) typically reveal emission lines at rest .Indeed, six out of seven emission lines detected by our line-search algorithm can be straightforwardly identified as such.All emission line identifications are summarized in Table 1.The observed wavelengths λ obs of four of these features are consistent with strong Lyα emission from Ne X (λ obs = 12.15 Å), Mg XII (λ obs = 8.4 Å), Si XIV (λ obs = 6.18Å), and S XVI (λ obs = 4.75 Å).As shown in the above references, these ions are often observed in rest-emission in X-ray binaries accreting above the Eddington limit.
This leaves three lines to be identified, at 5.04 Å, 5.96 Å, and 9.50 Å.The first wavelength, 5.04 Å, coincides with the resonance line of Heα-like S VI, which fits with the detection of the Lyα emission line of S XVI.The 9.50 Å emission line fits with Lyδ emission of Ne X, which we know is present from the Lyα line.Finally, no emission lines appear to be present within 0.1 Å of the 5.96 Å feature and we do not identify this line.While a feature is visible in the HEG spectrum at this wavelength, the MEG feature is very close to a large instrumental residual associated with a large feature in the instrument response.It might therefore be a spurious detection.
In addition to these seven emission lines present in both detectors, the iron fluorescence complex below 2 Å is clearly visible in the HEG spectrum.We measured centroid wavelengths of 1.93 Å, 1.85 Å, and 1.77 Å, consistent with Fe Kα, Heα-like Fe XXV, and Fe Kβ, all at rest.Such emission lines at rest are also seen in all Chandra gratings spectra of the ten HMXBs analysed in the overview work by Torrejón et al. (2010).Analysing a NuSTAR observation early on in the outburst of Sw J0243 (during the sub-Eddington phase), Bahramian et al. (2017) report the presence of a broad (σ = 0.3 ± 0.1 keV) Gaussian iron line at 6.42 ± 0.07 keV (∼ 1.931 Å).This could either be the same feature as present in the HEG spectrum, only not resolved into the three individual lines.Alternatively, the NuSTAR feature might be a relativistically broadened reflection feature, which transitioned into the three narrow lines we observe as the mass accretion rate increased.Finally, the observed HEG structure could arise from two absorption features on top of a broad emission feature.However, we could not find an satisfactory fit of the HEG spectrum with such a combination of emission and absorption.Combined with the accurate match between the centroid wavelengths and the expected iron line energies, we conclude that the final option is unlikely.

Emission line modelling
To further analyse the emission lines, we performed spectral fits with an emission line model added to the spline continuum.Such physical modelling can provide insights in the properties of the emitting gas, such as temperature and ionisation state.By comparing different models, the origin of the ionisation can also be constrained.We performed line modelling using two models in xspec: BAPEC, which models a velocity-broadened, shock-ionised gas, and PHOTEMIS, describing emission from a photo-ionised plasma.As we find no evidence for red or blueshifts in the identified emission lines, we freeze the redshift parameter in both models to zero.In both models, we assume Solar abundance ratios.We fit the MEG and HEG spectra simultaneously, both with their own spline continuum model, and keep the line model parameters tied between the spectra.
The BAPEC model provides the best fit for a velocity broadening of v = 1100 +200 −340 km/s and a temperature kT = 0.68 ± 0.03 keV, with an improvement in C-statistic of ∆C = 47.3 for three extra free parameters (including normalisation).Visual inspection of the residuals reveals that this improvement largely arises from fitting the 12.15 Å Ne X line, while no other emission lines are fitted.An issue with the BAPEC model is the presence of a significant pseudo-continuum of lines, which cannot be fitted to the non-physical continuum of the Chandra observations.There- fore, this systematic effect prevents a more accurate fit of the spectra.
The PHOTEMIS model, however, provides a formally better description of the emission lines with a ∆C = 146.8for three additional free parameters, for an ionisation parameter r log ξ = 2.77±0.05and a turbulent velocity v = (2.2±0.2)×10 2 km/s.The five strongest lines in the model are located at 12.15 Å, 8.4 Å, 6.18 Å, 5.04 Å, and 4.75 Å, fitting the Ne X, Mg XII, SI XIV, S XV, and S XVI features, respectively.The only detected narrow features not described in the model (see Table 1) are the unidentified 5.96 Å feature, and the three iron lines below 2 Å, which are located outside the fitted wavelength range.While this suggests the emitting gas could be photo-ionised rather than shock-ionised, the comparison with the BAPEC model is complicated by the systematic pseudo-continuum issues in the latter.

No outflow scenario
The identification of the detected absorption lines is more ambiguous than that of the emission lines, as the possible presence of blue-shifted lines greatly increases the feasible line identifications.Here, we will first focus on an identification scenario where no outflow was present and all lines are at rest.In this scenario, the absorption lines are either dominated by only Fe lines or by a combination of Fe and Mg lines.
Out of the seven detected absorption lines, listed in Table 2, four can be identified with the same ions in both the only-Fe and Fe+Mg interpretations.Iron absorption can account for the features observed at 8.21 Å, 9.45 Å, and 9.80 Å, where several transitions of respectively Fe XXI-XXIIV, Fe XX-XXII, and Fe XIX-XXII are located.The 3.00 Å absorption line could be associated with the Lyα transition of Ca XX at rest.However, while iron is often invoked to explain observed narrow lines in X-ray binaries (e.g Pinto et al. 2016) and the clear iron fluorescence lines show that iron is present in Sw J0243 (c.f.Section 4.2.1),Ca XX is not typically observed in these systems.Additionally, the 3.00 Å feature appears predominantly present in the MEG spectrum.Therefore, there is a possibility that this line is merely a spurious detection.
Assuming no blue-shifted absorption, we can link the other three detected lines with either iron (the Fe interpretation in Table 2) or magnesium (the Fe+Mg interpretation).For the former interpretation, the 7.35 Å and the 7.90 Å lines can be associated with Fe XXII-XXIV and Fe XXII-XXIII, respectively.This would leave the 6.50 Å feature unidentified.For the Fe+Mg interpretation, this 6.50 Å feature could arise from Mg XII, while the 7.35 Å and the 7.90 Å lines would be Mg XI and the Hβ-like resonance of Mg XI, respectively.However, this Fe+Mg interpretation has several caveats: while the wavelengths match up and the 8.4 Å Mg XII Lyα emission line shows that Mg is present, it is unexpected that Mg XII would be observed at rest in both emission and absorption (Grinberg et al. 2017).Furthermore, the 6.50 Å Mg XII and 7.35 Å Mg XI lines are relatively weak transitions, and one would therefore expect to see other or a larger number of Mg absorption lines instead.Finally, the rest and observed wavelengths do not match up perfectly in this interpretation.
For the above two interpretations, we can calculate the cumulative improvement in fit statistic by adding the combined ∆C values for each identified line.The six identified lines in the Fe interpretation yield a cumulative ∆C of 127.1 (124.5) for 500 (2000) km/s, while the seven identified lines in the Fe+Mg interpretation sum up to ∆C = 141.7 (147.0) for 500 (2000) km/s.

Outflow scenario
Alternatively, we consider a scenario where a selection of the absorption features are identified as blue-shifted transitions with the same outflow velocity.The firm detection of Ne, Mg, Si and S in emission aids in this approach, as it provides a starting point to identify lines that might be expected in absorption.In fact, in the discovery of ultra-fast outflows (UFOs) in ULXs, Pinto et al. (2016) observed many of the observed rest emission features in absorption with the same blue-shift.For instance, in these ULXs, the Ne X Lyα rest emission line that is also present in Sw J0243, is also observed in absorption with a ∼ −0.2c velocity shift.
To test for a similar scenario in Sw J0243, we calculated the blue-shifts required to explain every combination of an observed absorption line and an observed higher-wavelength emission line (excluding the unidentified 5.96 Å emission feature).In the case of an outflow, we would expect a similar blue-shift to appear for a number of such pairings.Indeed, for an outflow velocity of ∼ 0.22c, the 6.50 Å, 7.35 Å and 9.45 Å absorption lines can be linked to, respectively, the observed 8.4 Å Mg XII Lyα line, 9.50 Å Ne X Lyδ line, and 12.15 Å Ne X Lyα line (see Table 3).This scenario provides a seemingly more feasible explanation of the 6.50 Å feature than that in Section 4.3.1, while the required velocity is similar to that of UFOs in ULXs with an unknown accretor (Pinto et al. 2016) and the ULX pulsar NGC 300 ULX-1 (Kosec et al. 2018a).The cumulative ∆C for these three blueshifted lines is 57.7 (56.8) for a velocity width of 500 (2000) km/s.
If Sw J0243 launches an outflow with a velocity of 0.22c, we can ask two more questions.Firstly, why do we only observe the Ne X and Mg XII emission lines in absorption as well?Shifting the two S XVI and the Si XIV emission lines by the same velocity returns wavelengths of 5.45 Å, 4.15 Å, and 4.44 Å.Out of these, hints for an absorption feature can only be seen around ∼ 4.15 Å in the HEG spectrum.However, this is not a convincing feature, and no hints of a line are present at the other two wavelengths.Given the strength of the Ne X and Mg XII lines in HMXBs in general, and in Sw J0243 specifically, it is however not surprising that these ions are most easily detected in blue-shifted absorption.
Secondly, we consider whether any of the other absorption features might be associated with 0.22c blue-shifted lines from species not observed in emission.Shifting these four remaining absorption lines, only the 8.21 Å feature yields a possible match; its shifted wavelength of 9.33 Å is similar to the 9.31 Å forbidden transition of Heα-like Mg XI.However, it appears unlikely that only this forbidden line is observed, while other stronger transitions are not seen.Therefore, the only direct evidence for the outflow consists of the three absorption lines listed in 3, and we interpret the remaining absorption features as in Section 4.3.1.

DISCUSSION AND CONCLUSIONS
We have reported high-resolution Chandra X-ray spectroscopy of the super-Eddington outburst of Sw J0243.A search for narrow emission and absorption features reveals a number of both, present in both the HEG and MEG spectrum.The emission lines can be identified with Fe, S, Si, Mg, and Ne ions at rest.The absorption features can either be interpreted to be all at rest (from Fe and possibly Ca and Mg), or a combination of some lines at rest and three blue-shifted Mg and Ne absorption lines at v ≈ −0.22c.Here, we briefly review our method, discuss the possible outflow in the context of ULXs and close-by BeXRBs, and finally present future improvements for the study of outflows from BeXBRs during super-Eddington phases.

Line-search method
As shown in Figure 2, the continuum of the Chandra spectrum differs greatly between detectors and deviates from the shape measured by Swift.However, the iron fluorescence complex below 2 Å, i.e. around 6.5 keV that is observed in all HMXBs with Chandra gratings observations (Torrejón et al. 2010), is clearly detected.Similarly, the Ne X and Mg XII Lyα emission lines at respectively 12.15 Å and 8.4 Å, often observed in neutron star HMXBs (e.g.Grinberg et al. 2017;Koliopanos & Vasilopoulos 2018;La Palombara et al. 2016), are visible in the spectra even by eye (c.f. the top panel of Figure 3).This suggests that indeed, while the continuum is affected by the observing setup, narrow features remain detectable (Schulz et al. 2009).In addition, the consistency checks of the spline continuum model (e.g.Appendix A) show that none of the detected lines or conclusions are due to the non-physical nature of this model.Our adopted line-search method follows the rationale first used by Pinto et al. (2016), and later applied to both XMM-Newton RGS and Chandra observations by Degenaar et al. (2017), Kosec et al. (2018a), Kosec et al. (2018b), andvan den Eijnden et al. (2018a).However, as discussed in more detail in van den Eijnden et al. (2018a), estimating formal significances of detected features is challenging.The significances shown in Figure 3 are single-trial values, while estimating the number of independent trials is difficult: fitting a Gaussian line at neighbouring wavelength gridpoints is not independent, as the Gaussian width exceeds the resolution of the grid.For a velocity width of 500 km/s, a Gaussian line covers between 4 and 28 wavelength bins of 0.01 Å width in its one-sigma range, where the number varies since a constant velocity width translates to a variable width in wavelength-space.For the 2000 km/s search, these numbers are multiplied by four.Therefore, while we fit the normalization of a Gaussian line at 1161 wavelength bins, a much smaller -but difficult to estimate precisely -number of those trials are truly independent.
Therefore, we also performed Monte-Carlo simulations of the continuum shape to estimate how likely random fluctuations can reproduce the observed excesses (e.g van den Eijnden et al. 2017;Kosec et al. 2018a).Secondly, we opted for an independent analysis and subsequent comparison of the HEG and MEG spectrum, since statistical fluctuations or response-effects are less likely to show up in both detectors at the same wavelength.Thirdly, searching with two different velocity widths decreases the probability of statistical fluctuations being identified as a line: while a small number of bins fluctuating either above or below the continuum by chance might mimick a narrow line, it would not be identified as such when searching with a broader velocity width.Finally, it is important that any possible line can be identified within a coherent physical picture of the system.For this reason, we do not identify for example the possible line at 5.96 Å.

An ultra-fast outflow from Sw J0243?
While the combination of detected absorption lines can be interpreted as simply Fe (and Mg) ions at rest, a more interesting possibility is the presence of an outflow suggested by the presence of absorption features at a ∼ 0.22c blue-shift from Mg Lyα, Ne X Lyα, and Ne X Lyδ -which are all observed in emission.The presence of such an outflow during the super-Eddington regime fits both theoretical predictions and simulations (e.g.Shakura & Sunyaev 1973;Ohsuga & Mineshige 2011;McKinney et al. 2014McKinney et al. , 2015;;Hashizume et al. 2015;King & Muldrew 2016) and observational work (e.g. Lee et al. 2002;Pinto et al. 2016;Kosec et al. 2018a;Allen et al. 2018).
In the outflow scenario, the blueshifted absorption lines originate from the same ions as several rest-emission lines.This combination of rest emission and blue-shifted absorption from the same ions might arise from the outflow, seen in absorption, shocking with the surrounding medium, seen in emission.Such a scenario is also invoked in Pinto et al. (2016), where the emission can be modeled as a shockionised gas.Alternatively, the emission might arise from a different region in the system, instead of the outflow but with the same ions present (Grinberg et al. 2017), such as the accretion flow.
The Sw J0243 high-resolution X-ray spectrum is similar to that of other Galactic and SMC BeXRBs in several aspects.For instance, the XMM-Newton RGS spectrum of SMC X-3 during its super-Eddington state also shows restemission lines of Ne X, Mg XII, and Fe XXIII-XXIV (Koliopanos & Vasilopoulos 2018).Additionally, a possible blueshifted Mg XII absorption line is detected, which would imply an outflow with a ∼ 0.07c velocity.The presence of Mg XII in both rest emission and blue-shifted absorption mirrors our outflow interpretation for Sw J0243.At slightly sub-Eddington X-ray luminosity (L X ∼ 10 38 erg/s), the BeXRB SMC X-2 shows both rest emission lines of Ne X and Si XIV (La Palombara et al. 2016), as we also identify in Sw J0243 (in addition to several other rest emission lines not seen in Sw J0243 due to the difference between the HEG/MEG and RGS bandpasses).However, no hints for absorbtion features or an outflow are present.Even further down in the sub-Eddington regime, Grinberg et al. (2017) report a plethora of rest emission lines in Vela X-1, including the Fe, Si, Mg, and Ne species identified in Sw J0243.But again, no outflow is detected, despite the high-quality observations which would likely reveal an outflow similar to the one we possibly detect in Sw J0243.
Given the super-Eddington accretion rate of Sw J0243 during the Chandra observation, ULXs form a second interesting source class for comparison.Kosec et al. (2018a) reported the detection of a possible outflow from the ULX pulsar NGC 300 ULX-1.This outflow was observed through the identification of blueshifted O XVII and O XVIII absorption lines in RGS spectra, which fall outside the Chandra bandpass.The inferred velocity of 0.22c is however consistent with the velocity of the possible wind in Sw J0243.In addition, Pinto et al. (2016) present the discovery of ∼ 0.2c outflows from the unclassified ULXs NGC 1313 X-1 and NGC 5408 X-1; interestingly, next to the similarity of the wind velocity, both sources show a combination of emission lines at rest, absorption lines from the same species at a blue-shift, and additional rest absorption lines.This mimicks exactly our outflow scenario for Sw J0243.Finally, an outflow with a higher velocity of ∼ 0.34c was tentatively claimed in NGC 5204 X-1 by Kosec et al. (2018b), but this result awaits confirmation.
Around the time of the Chandra observation of Sw J0243, the source also launched a radio jet (van den Eijnden et al. 2018b).Radio observations taken four days later show an optically-thin radio spectrum, implying that during this super-Eddington state, the jet consisted of discrete ejecta (e.g.Fender 2006) 4 .While there is no simultaneous radio and high-resolution X-ray coverage, Sw J0243 remained in a very similar state between the jet and possible wind detection.Therefore, we deem it likely that both an ultra-fast outflow and a jet are launched at the same time.Such behaviour is observed more often during the super-Eddington regime in other sources: black holes and Z-sources also show winds and jets during the same super-Eddington accretion states (Homan et al. 2007(Homan et al. , 2016;;Allen et al. 2018).However, while both winds and jets have been inferred in ULXs (e.g.Middleton et al. 2013;Cseh et al. 2014;Kaaret et al. 2017), these have never been observed in the same target, let alone at the same time.Given the difficulty to find these outflows (e.g.Kosec et al. 2018b), due to the large distances to ULXs, Galactic BeXRBs offer a new avenue to explore the connection between winds and jets in the super-Eddington regime.
One particularly puzzling aspect of the jet launched by Sw J0243 is its faintness compared to fast-spinning, weakly-magnetized accreting neutron stars at similar super-Eddington accretion rates (van den Eijnden et al. 2018b).This can naively be explained by the slow spin of Sw J0243 (e.g.Parfrey et al. 2016); however, at lower accretion rates, the radio brightness of Sw J0243 is consistent with fasterspinning neutron stars, which implies a more complicated picture (van den Eijnden et al. 2019).Possibly, the presence of an ultra-fast disk wind during the super-Eddington phase of the outburst can regulate the jet power; depending on its launch radius, a wind might decrease the mass accretion rate in the inner accretion flow and reduce the matter available to form the jet.Alternatively, it might carry away excess angular momentum and reduce the jet power.A similar interplay between the wind and the jet has earlier been proposed to explain the jet-wind dichotomy in GRS 1915+105 (Neilsen & Lee 2009, see also Díaz Trigo & Boirin 2013).Since the jet in Sw J0243 is consistent with the population of other NS jets at lower X-ray luminosity, this scenario assumes that the ultra-fast outflow was driven by the super-Eddington accretion rate and disappeared as the outburst decayed.
A scenario where the presence of a strong wind outflow regulates the jet power does not occur in Z-sources: these sources launch powerful jets with the highest radio brightness of any type of accreting neutron star.However, while the winds in these systems can carry away significant amounts of mass (Ponti et al. 2012;Allen et al. 2018), they generally do not reach velocities similar to those in ULXs and inferred here for Sw J0243 (i.e.maximally one per cent of the speed of light; Díaz Trigo & Boirin 2013).In addition, the inner accretion flow and jet launching regions differ greatly between weakly-magnetized Z-sources and the more strongly-magnetized ULX pulsars and BeXRBs (e.g Mushtukov et al. 2017;Walton et al. 2018b).Finally, while Z-sources accrete close to or above the Eddington limit, our Sw J0243 Chandra observation was taken at one order of magnitude higher X-ray luminosity.Therefore, any coupling between the (super)-Eddington winds and jets would not necessarily be the same in Sw J0243, ULX pulsars, and Zsources.
out the entire outburst show that this radio emission can not originate from either a stellar or disk wind.

Future observations
While the presence of an ultra-fast outflow from Sw J0243 fits with the Chandra high-resolution X-ray spectrum, the continuum issues complicate a full analysis and more detailed physical modelling than described in Section 4.2.2.The possible wind detection does however showcase the power of studying Galactic BeXRBs for understanding pulsating ULXs and their outflows.Therefore, future observational campaigns combining radio, X-ray and UV observations would be highly valuable: dense radio monitoring can track the jets, while high-resolution X-ray and UV spectra (from i.e. the Hubble Space Telescope) taken at different phases in the outburst can track the onset and evolution of any wind outflow.Through such detailed monitoring, the relation between the jet and wind can be studied as well, for instance aiming to understand if and how the (presence of the) wind might influence the jet power.These future observations can also reveal how commonly super-Eddington BeXRBs launch a wind and jet simultaneously, to better understand the expected outflow properties of ULX pulsars.

Figure 1 .
Figure 1.Multi-band light curves of the 2017/2018 outburst of Sw J0243.The top panel shows the Swift/BAT 15-50 keV count rate, the middle panel shows the MAXI /GSC 4-10 keV count rate, and the bottom panel shows the ASAS-SN V-band magnitude (Vega Mag).The Chandra epoch is indicated with the red line in the top panel.

Figure 2 .
Photons cm 2 s 1 keV 1 Figure 4.The grey dotted lines indicate the same identified features as in Figure 3.For visual clarity, we removed three wavelength ranges (shown by the grey bands) where large instrumental features in a single detector yield extreme ∆C values.In Appendix B, we also show the uncombined ∆C results for each individual detector (Figures B1 and

Figure 4 .
Figure 4.The summed improvement in HEG and MEG fit statistic ∆C for the 500 km/s (top panel) and the 2000 km/s (bottom panel) velocity width line search.The grey dotted lines show the identified emission and absorption features, while the grey bands show regions containing instrument features that are excluded to improve the clarity of the figure.Finally, the dashed and dash-dotted blue lines show the simulated 2 and 3σ confidence levels.For the uncombined fit statistic improvement of the individual detectors, see Appendix B

Figure A1 .Figure A2 .Figure B1 .
Figure A1.The first consistency check of the spline continuum, shown the HEG spectrum below 7 Å: both panels show the results of the line search algorithm for different velocity widths (top: 500 km/s, bottom: 2000 km/s).The black curves use the spline continuum model, while the red and blue curves use a power law and a power law + blackbody continuum model, respectively.In both cases, and especially in the bottom panel, residual trends remain when using the physical continuum, which can artificially enhance single-trial significances.However, the individual narrow features appear for all continuum models.

Figure B2 .
Figure B2.Same as Figure B1 for the 2000 km/s line width search.

Table 1 .
Identification of the detected emission lines in Sw J0243.See Section 4.2.1 for details.

Table 2 .
Identification of absorption lines in Sw J0243 in the no-outflow scenario.See Section 4.3.1 for details.Note: for the Fe ions, multiple transitions fall close to the observed wavelength.Therefore we do not list a single rest wavelength.

Table 3 .
Identification of a selection of absorption lines in Sw J0243 in the outflow scenario.See Section 4.3.2 for details.The remaining absorption lines are interpreted as in Table 4.3.1 and 4.3.1.