© ESO 2020

1 Introduction

With a column density on the order of 10²⁴ cm⁻², IGR J16318−4848 is one of the most strongly intrinsically absorbed X-ray sources in our Galaxy. It was discovered during routine scans of the Galactic plane by the IBIS/ISGRI instrument on board the INTEGRAL satellite in 2003 January (Courvoisier et al. 2003; Walter et al. 2003). A follow-up XMM-Newton observation in 2003 February revealed the very strong absorption and strong K-shell emission lines of Fe and Ni (Matt & Guainazzi 2003). These authors presented a detailed study of the absorbing and emitting material in the vicinity of the source, and detected a very weak Compton shoulder of the emission lines from which they concluded an anisotropy of the absorber with an average column density of a few 10²³ cm⁻².

The high absorption column and the strong fluorescence lines, which dominate the source flux below 10 keV, were conf- irmed in theanalysis of 2006 August observations with Suzaku by Barragán et al. (2009). The upper limit on the strength of the Compton shoulder of the emission lines in these observations again strongly argues for an inhomogeneous and anisotropic absorber. A re-analysis of archival ASCA data revealed the source to be persistent (Murakami et al. 2003). Later INTEGRAL monitoring revealed IGR J16318−4848 to be a strongly variable hard X-ray source (Barragán et al. 2010a,b). A Hitomi observation in 2016 March constrained the ionization state of iron to be less than Fe IV and revealed line width broadening corresponding to 160 km s⁻¹ (Hitomi Collaboration 2018).

A possible optical counterpart for IGR J16318−4848 was proposed by Foschini et al. (2003). Dedicated optical and near-infrared (NIR) observations confirmed and refined the detection of the optical companion, which Chaty & Filliatre (2004) and Filliatre & Chaty (2004) concluded was a sgB[e] star with a distance estimate ranging from 0.9 to 6.2 kpc. The infrared observations indicate that the star is probably embedded in dense and absorbing circumstellar matter. Moon et al. (2007) and Chaty & Rahoui (2012) suggest that this source-intrinsic absorber contains cold and warm dust. The interstellar absorption is two orders of magnitude lower than the circumstellar absorption (Chaty & Filliatre 2004). More recently, Fortin et al. (2020) performed stellar atmosphere and wind modeling using VLT/X-shooter observations and concluded an inclination of the system higher than 76° from broadband NIR and optical data. They also observed optical forbidden lines originating from an optically thin medium, possibly a disk wind. The nature of the compact object is still under debate (e.g., Walter et al. 2003; Barragán et al. 2009). The compact X-ray source and the sgB[e] companion place this source among the high mass X-ray binaries (HMXBs) and Iyer & Paul (2017) propose a ~80 d orbit based on intensity variations in the Swift/BAT light curve.

In this paper we discuss the analysis of simultaneous observations of IGR J16318−4848 made on 2014 August 22–23 with the Nuclear Spectroscopic Telescope Array (NuSTAR) and with the X-ray Multi-Mirror Mission (XMM-Newton). This is the highest signal-to-noise ratio (S/N) observation of the source made to date, covering the 3–60 keV band. In Sect. 2 we give an overview of the observations and our data reduction approach. Section 3 describes our modeling of the time-averaged spectrum with empirical and physical fluorescence line models, including models for a dusty absorber. We study the time variability of the spectral shape in Sect. 4. We discuss our results and conclude in Sects. 5 and 6.

2 Observations and data reduction

2.1 NuSTAR

NuSTAR (Harrison et al. 2013) was launched on 2012 June 13 as a NASA Small Explorer mission. NuSTAR has two co-aligned grazing incidence X-ray telescopes, which focus X-rays on two focal plane modules, FPMA and FPMB. NuSTAR provides a usable energy range of 3–79 keV with a field of view of 13′ × 13′ at 10 keV.

NuSTAR observed IGR J16318−4848 on 2014 August 22 and 23 (ObsID 30001006002). We processed these data with the standard NUSTARDAS pipeline version 1.8.0 with CalDB version 20191219. The source regions were circles of 80″ radius, while the background regions were three circles of 140″ radius each, placed at the remaining corners of the field of view. We included data taken in SCIENCE_SC mode, where the attitude reconstruction is performed only by the spacecraft star trackers and not by the star trackers on the optical module. The position reconstruction in this observing mode is less precise and it is triggered, for example when the optical module star trackers are blinded by the moon or sun. We carefully checked that the spectral shape of the SCIENCE_SC mode data is in very good agreement with the standard SCIENCE mode data. We therefore combined the SCIENCE and SCIENCE_SC data, adding ~ 8 ks more exposure per FPM to yield a total net exposure of ~64 ks per FPM. All times were barycentered with the FTOOL barycorr. We used the Interactive Spectral Interpretation System (ISIS v. 1.6.2; Houck & Denicola 2000) for all spectral and timing analyses. Uncertainties are given at the 90% confidence level for one parameter of interest unless stated otherwise.

2.2 XMM-Newton

The XMM-Newton satellite (Jansen et al. 2001) was launched on 1999 December 10 by the European Space Agency. It carries three main instruments: the European Photon Imaging Camera (EPIC; Strüder et al. 2001), the Reflection Grating Spectrometer (RGS; den Herder et al. 2001), and the Optical Monitor (OM; Mason et al. 2001).

XMM-Newton observed IGR J16318−4848 on 2014 August 22 and 23 for ~64 ks (ObsID 0742270201), simultaneous with NuSTAR. Using the Science Analysis System (SAS) version 18.0.0, we extracted the EPIC-pn, -MOS1, and -MOS2 spectra and light curves following standard reprocessing and screening criteria and employing circular source and background regions of 40′′ and 60′′ radius, respectively.

We use single and double events for EPIC-pn and single events for EPIC-MOS for our spectral analysis to ensure the most accurate charge transfer efficiency and response calibration at low energies as required for constraining the continuum absorption and emission line profile. The loss of hard X-ray events due to the pattern restriction is tolerable because of the overlap with the high quality NuSTAR data. We exclude spectral channels below 3 keV where the data are strongly absorbed and the spectrum has a very low S/N. We also ignore the RGS as the energy range accessible to the RGS is also too strongly absorbed to obtain a significant S/N value.

3 Time-averaged spectroscopy

For our spectroscopic analyses we used the 3–10 keV and 5–60 keV spectra of the XMM-Newton/EPIC-pn and NuSTAR/ FPMs, respectively. The NuSTAR/FPMA and -B spectra were jointly rebinned requiring a minimum S/N of 15 and adding at least 2 and 4 channels below and above 40 keV, respectively. The XMM-Newton/EPIC-pn spectrum was rebinned adding at least four channels and additionally requiring a minimum S/N of 5 and 10 below and above 6 keV, respectively. We introduce a constant gainshift per instrument to account for possible remaining calibration uncertainties of the energy scale (see, e.g., Briel et al. 2005; Kitaguchi et al. 2011, for details on the charge transfer modeling in XMM-Newton and NuSTAR, respectively). This approach has proven necessary because in physically motivated photoionization models, line and edge energies are not free parameters but reference values from some corresponding atomic database. Therefore, instrument calibration uncertainties cannot be accounted for by the spectral model. Furthermore, we introduced flux cross-calibration constants with respect to NuSTAR/FPMA.

3.1 Empirical emission line modeling

Following previous analyses (e.g., Ibarra et al. 2007; Barragán et al. 2009), we first describe the data with a power law with an exponential cut off (called cutoffpl in ISIS/XSPEC) with photon index Γ and folding energy E_fold, accounting for photoelectric absorption with the model tbabs v.2¹, parameterizing absorption using the equivalent hydrogen column density N_H. We used abundances and cross sections according to Asplund et al. (2009) and Verner et al. (1996), respectively. We left the iron abundance free to model the iron K-edge. The absorber is further modified by an optically thin Compton scattering component (cabs) with the same column density as the neutral absorber. Similar to all previous observations, the spectrum shows very strong fluorescence lines of Fe Kα, Fe Kβ, and Ni Kα (e.g., Matt & Guainazzi 2003). Based on the Hitomi results (Hitomi Collaboration 2018), we modeled these with Gaussian emission lines and fixed the centroids at 6.404, 7.058, and 7.478 keV, respectively. Since the line broadening seen with Hitomi cannot be resolved with XMM-Newton or NuSTAR, we fixed all line widths to 3.4 eV. In XSPEC-like notation, our model is defined as $$\begin{array}{l}S(E)=gainshift\times detconst\times tbnew\times cabs\\ \times \left(cutoffpl\times ComptonShoulder+3\times gauss\right).\end{array}$$(1)

The time-averaged spectrum and best-fit model are shown in Fig. 1. Table 1 lists the best-fit parameters. The fit confirms the general results of earlier papers that the spectrum is dominated by a very strong fluorescence Fe Kα line and a strong Fe K-edge that overlays an exponentially cut off power law. The flux below 10 keV is fully dominated by the emission lines. Because the different components are treated as independent of each other, however, the physical interpretation of these parameters with respect to the absorber geometry and physics is difficult. In the next sections we therefore attempt to model the spectrum with more physics-based models. We start in Sect. 3.2 with a study of the strength of the Compton shoulder and then model the spectral parameters with a simple photoionization model (Sect. 3.3) and with a dusty absorber motivated by the optical behavior of the system (Sects. 3.4 and 3.5).

Fig. 1 Panel a: time-averaged counts spectrum of IGR J16318−4848 with best-fit cutoff power law model for NuSTAR/FPMA (gray), FPMB (gold), XMM-Newton/EPIC-pn (blue), and XMM-Newton/EPIC-MOS1 and -MOS2 (dark and light green, respectively). Panels b and c: residuals and ratio for the best-fit model.

3.2 Compton shoulder

The Compton shoulder is a spectral feature that originates from the down-scattering in energy of the fluorescence photons off electrons in the ambient material. It is an important diagnostic tool to learn about the structure and geometry of the absorber. Its possible presence in the X-ray spectrum of IGR J16318−4848 has already been discussed in previous analyses (Matt & Guainazzi 2003; Barragán et al. 2009).

In order to constrain the strength of a possible Compton shoulder, we add a Gaussian emission feature at 6.3 keV with a width of 50 eV to the model and then determined the statistical significance and flux limits of this component. Adding the Compton shoulder component only slightly decreases χ² from 1069.2 to 1036.4, while introducing one additional free parameter. Following Protassov et al. (2002), we estimate the significance of this putative Compton shoulder component in a Monte Carlo simulation where we simulate the X-ray spectrum without a Compton shoulder and then estimate how often a false positive detection of a Compton shoulder is seen with a Δχ² that is larger than that found in the data. We find no false positive detections of a Compton shoulder in 10⁵ Monte Carlo runs (i.e., the existence of a Compton shoulder is required at least at 4.4σ confidence by the joint NuSTAR and XMM-Newton data).

We find single-parameter 90% uncertainty of the ratio of the flux of the Compton shoulder with respect to that of the Fe Kα line (f_CS) of 11 ± 4%. This value is to be compared with the prediction of the strength of the Compton shoulder for the transmission through a centrally illuminated spherical absorber. Specifically, for column densities of ~ 2 × 10²⁴ cm⁻² Matt (2002) estimate that the Compton shoulder ratio should be ~45% (i.e., a factor of ~3–7 higher than observed in our data). This discrepancy is in agreement with earlier measurements. To explain this discrepancy with the data, previous authors have suggested low covering fractions (Matt & Guainazzi 2003) or highly structured inhomogeneous absorbers (e.g., Barragán et al. 2009) as a possible explanation. As we show in the next section, however, both explanations are difficult to reconcile with the absorption seen in the continuum.

3.3 Constraints set by photoionization

The main inconsistency of the empirical models that we applied so far and that were also used in earlier modeling of the spectrum of IGR J16318−4848 is that absorption and fluorescence emission are treated separately; this means that the absorption model, for which N_H and some element abundances are free to vary, provides no feedback on the emission lines, which are independent, additive components. This approach is mainly due to a lack of available models for XSPEC. One of the few models that aims at self-consistent modeling of photoabsorption and fluorescence is the XSPEC model PEXMON (Nandra et al. 2007). This model combines in an empirical way reflection from cold matter (Magdziarz & Zdziarski 1995) with line strengths by George & Fabian (1991) and Compton shoulder simulations by Matt (2002). However, PEXMON assumes a pure power law as incident spectrum with photon index 1.1 < Γ < 2.5, and is therefore mostly applicable to Seyfert 1 galaxies, but not necessarily to HMXBs such as IGR J16318–4848. We therefore have to perform additional checks using more physical modeling to see whether our empirical model is self-consistent.

XSTAR (Kallman & Bautista 2001), used here in version 2.53, is a software package for the self-consistent modeling of the ionization structure and emission spectra of photoionized gas clouds. The code assumes a cloud of gas irradiated by a central source. The physical parameters that can be varied by the user are the chemical composition of the gas, its temperature, density, pressure, and the strength and energy distribution of the irradiating source. In the determination of the ionization structure, all relevant physical processes such as photoionization and collisional ionization are considered. A covering fraction parameter allows for the realization of both spherical and slab-like geometries.

In order to model the structure and emerging spectrum of a putative photoionized plasma around IGR J16318–4848, we ran a grid of XSTAR simulations in the parameter range of interest for IGR J16318–4848. We used an assumed distance of 5 kpc, well within the distance estimate interval and also supported by Fortin et al. (2020) with a corresponding source luminosity of 1.33 × 10³⁶ erg s⁻¹ in the 13.6 eV to 13.6 keV energy range. The cloud that surrounds the source is irradiated with a power law with an exponential cut off derived from the empirical fits as the input spectrum. We caution that the derivation of the shape of the unabsorbed incident continuum depends on the choice of the iron abundance in our empirical fit. The iron abundance in our XSTAR simulation is therefore not allowed to vary in order to avoid inconsistencies in the model assumptions. We cannot rule out that a systematic exploration of different iron abundances along with continuum shapes could result in a better fit to the data. This is, however, beyond the scope of this work. We also note that the observational data do not completely span the energy range required by XSTAR, so our spectra and luminosity rely on the assumption that the spectral shape does not deviate significantly from an extrapolation below 3 keV which is, however, likely to lead to a significant ultraviolet excess. Elements included in the simulation are H, He, C, N, O, F, Ne, Mg, Si, S, Ar, Ca, Cr, Mn, Fe, Co, and Ni with solar abundances (based on Grevesse et al. 1996). Furthermore, we set the hydrogen particle number density to 6.5 × 10¹⁰ cm⁻³, which we estimated from the column density of the empirical fits, and the estimated size of the emission region of 3 × 10¹³ cm from Matt & Guainazzi (2003).

We first determine the rough model parameter range by visually comparing example XSTAR runs with the data, and then run a finer grid for spectral modeling. As computational constraints force us to leave some of the XSTAR parameters fixed in the preparation of the grid, we vary only those parameters that have the largest effect on the resulting spectral shape. Specifically, we produced table models for XSPEC/ISIS on a model grid where we varied the covering fraction from 0.9–1.0, N_H from 5 × 10²³ cm⁻² to 5 × 10²⁴ cm⁻², and the logarithm of the ionization parameter (logξ) from 0 to 5. Here, ξ = L∕(nR²) (Tarter et al. 1969), where L is the sourceluminosity, n the particle density, and R is the absorber’s distance from the source.

The choice of the limited range of covering fractions in our table model is motivated by preliminary XSTAR calculations where the coveringfraction was varied over a range from 0.1–1.0. Figure 2 shows some of these spectra as an example to illustrate how the line flux with respect to the continuum depends on the covering fraction. At low covering fractions, a significant part of the incident radiation is transmitted, which dilutes the line. These early calculations show that the covering fraction has to be very close to 100% to reproduce the overall strong absorption and the strong emission lines with respect to the observed continuum. This is in contradiction to the conclusions made in earlier analyses of the X-ray spectrum of IGR J16318–4848 (e.g., Matt & Guainazzi 2003), where a lower covering fraction was inferred from the observed Compton shoulder flux and equivalent width of the Fe Kα line.

To compare these models with the data and explore the full parameter space, we employed Markov chain Monte Carlo (MCMC) calculations with affine-invariant ensemble sampling (see Goodman & Weare 2010; Foreman-Mackey et al. 2013, for details), employing 300 walkers for six free parameters and 30 000 iterations. This approach is motivated by the complexity of physical photoionization models and their parameter space, where commonly used χ² -minimization algorithms are prone to get stuck in local minima. When we call a set of parameters “best fit” in the following, we refer to the most probable parameter combination unless noted otherwise. To reduce the number of free parameters, we fixed the detector flux and gainshift calibration constants to the values obtained from the empirical fits (Table 1). This is legitimate because these parameters are model independent and can be accurately determined from the earlier fits. Since fluorescence emission from ionization stages less than Ni IX is not included in XSTAR, the neutral Ni Kα fluorescence line was modeled by an additional Gaussian emission component with width fixed to 3.4 eV and energy fixed to 7.47 keV (Bearden 1967). In XSPEC-like notation, our fit model is given by $$S(E)=gainshift\times detconst\times \left(XSTAR+gauss\right).$$(2)

Despite the large volume of parameter space covered, pure single-zone photoionization modeling fails to result in a satisfactory description of the data. Figure 3 shows the best fit from the MCMC. The low quality of the fit renders listing the corresponding parameters unnecessary. Although the spectral shape of the model looks similar to the data (a strong Fe line and approximately the correct continuum curvature), neither the Fe Kα/Kβ line ratio nor the depth of the iron K-edge (and possibly also at the nickel K-edge at ~8.34 keV) are described correctly. We note that the mismatch of the hard continuum in this fit is merely a normalization issue since the fit is dominated by the fluorescence lines. We note that even though numerical issues prohibit a further exploration of the parameter space, the best fit is of such low quality that it is very unlikely that opening up further parameters would yield a satisfactory description of the data.

We note, however, that there is a very clear indication from the empirical fit and from the photoionization modeling that the medium responsible for the emission of the Fe Kα emission line and for the formation of the Fe K-edge is for some reason overabundant in Fe.

Fig. 2 Simulated XSTAR spectra for a wider range of covering fractions. Other model parameters are NH = 1024 cm−2 and log ξ = 1. Vertical arrows indicate the maximum line flux for clarity.

Fig. 3 Best fit of the single XSTAR model. Colors and panels are defined in Fig. 1.

3.4 Dust absorption modeling

The large mismatch between the Fe Kα line strength and the Fe K-edge depth, in combination with a very weak Compton shoulder that implies a low electron column, is difficult to reconcile with a gaseous absorber of solar abundances. In the previous section we also showed that within the constraints set by a strong Fe Kα line from neutral iron it is not possible to reconcile these observations with ionization effects. In this section, we study another possibility to explain the X-ray spectral shape: absorption of X-rays in a dust layer in the source. This approach is motivated by broadband NIR observations which suggest the presence of cold and warm dust within the binary system (Moon et al. 2007; Chaty & Rahoui 2012; Fortin et al. 2020), and by the general observation that interstellar dust is strongly enriched in refractory elements (Dwek 2016).

Since the optical observations indicate that the dusty absorber is source-intrinsic, we ignore dust scattering which gives rise, for example, to dust scattering halos when the dust is at a significant distance from the source, and treat absorption in dust only (Corrales et al. 2016).

Specifically, our dust-absorption model is a multiplicative factor of the form $$\exp \left(-{\tau }_{\text{dust}}\right)=\exp \left(-N_{\text{H}}{\sigma }_{\text{dust}}\right),$$(3)

where σ_dust is the (energy-dependent) absorption cross section per H-nucleon.

We describe the dust absorption cross section in the vicinity of the Fe K-edge using laboratory measurements by Rogantini et al. (2018) for various Fe containing minerals such as crystalline olivine (Mg_1.56Fe_0.4Si_0.91O₄), troilite (FeS), and pyrrhotite (Fe_0.875S). We present our results using the Rogantini et al. (2018) cross section for olivine for two reasons. First, the knowledge of grain sizeand abundance models are most reliable for silicates and graphites (see, e.g., Draine 2003a for a review, and Draine 2003b; Clayton et al. 2003; Nozawa & Fukugita 2013 for details). Second, due to the limited resolution of NuSTAR and XMM-Newton, we found that our results are insensitive to the specific choice of dust composition. Rogantini et al. (2018) assume a power law grain size distribution after Mathis et al. (1977, hereafter MRN) and solar iron abundance according to Lodders & Palme (2009). According to their cross section normalization, there is one iron atom for each 2.5 olivine unit cells. We extrapolate these cross-sections outside of the tabulated energy range of 6.7–8.0 keV for the element composition of olivine based on the gas-phase cross sections by Verner & Yakovlev (1995). This approach is justified since specific solid-state effects, i.e., X-ray absorption fine structure, mainly lead to deviations from the gas-phase cross section in areas close to the absorption edge (Lee & Ravel 2005; Lee et al. 2009).

We note that the total iron column along the line of sight is the most robust quantity to compare with other analyses, because it is constrained directly by the depth of the iron K-edge. For consistency with older studies, we map this quantity to the equivalent hydrogen column density assuming an iron abundance given by Asplund et al. (2009). We point out, however, that in a sgB[e] system the element abundances may depend strongly on the line of sight; in particular, heavier elements that are highly depleted into dust will accumulate in the circumstellar disk around the companion star, such that N_H determined from any fit with such a model is in all likelihood not a good estimate for the true hydrogen column.

Absorption by the iron atoms in the olivine grains will result in the emission of the Fe fluorescence line. We estimate the fluorescence line flux, N_ph,i for i = α, β, from the dust absorber with a simple empirical model following Nagase et al. (1986, their Eq. (11)), $$N_{\text{ph},i}={\gamma }_i{\omega }_{\text{K}}{\displaystyle {\int }_{E_{\text{K}}}^{\infty }\frac{{\sigma }_{\text{Fe}\text{K}}(E)}{{\sigma }_{\text{olivine}}(E)}}{N}_{\text{ph}}(E)\left[1-\exp (-\tau (E))\right]\text{d}E,$$(4)

where N_ph(E) is the incident photon spectrum, ω_K = 0.351 is the fluorescence yield of iron (Hubbell et al. 1994), and γ_i is the branching ratio of the Kα to Kβ line. For neutral iron, as appropriate for a constituent in minerals, the intensity ratio Kβ/Kα is 0.132 (Han & Demir 2009). This quantity is expected to be at most slightly affected by solid state effects. In the Fe-band, self-absorption of Fe fluorescence photons in grains is negligible (Wilms et al. 2000, and Sect. 3.5).

To summarize, in order to describe the data with the dust model, we consider a pure dust absorber, which we describe in pseudo-XSPEC notation by $$\begin{array}{l}S(E)=gainshift\times detconst\times (\eta (N_{\text{ph},\alpha }(cutoffpl)\\ +N_{\text{ph},\beta }(cutoffpl))+\exp \left(-{\tau }_{\text{dust}}\right)\times (cutoffpl+gauss)) ,\end{array}$$(5)

where τ_dust and N_ph,α,β were defined in Eqs. (3) and (4), respectively. We note that in this model N_ph,α,β are convolution models, since the fluorescence line strength is determined from the overall spectral shape above the Fe K-edge through Eq. (4). In addition to these physics-motivated components, our model includes an empirical correction factor η that accounts forthe physical effects excluded in Eq. (4), such as self-absorption of the fluorescence radiation in the absorber, any residual effects of self-shielding by the grains, and geometrical dilution effects due to the unknown absorber geometry.

Since our olivine model does not include absorption by nickel, and thus the fluorescence yield cannot be calculated self-consistently, similar to our photoionization modeling (Sect. 3.3) we account for the Ni Kα line by adding a Gaussian emission line. The other fit parameters are N_ph and the spectral shape of the continuum. In Sect. 3.5 we check whether the results of our fits are consistent with detailed Monte Carlo simulations in which these effects are explicitly taken into account.

Figure 4 shows the spectrum with the most likely realization of the model. Even though we now have a physically much simpler model than the XSTAR one, the fit residuals look remarkably better than in our earlier attempts to model the spectrum. Despite the simplification of the olivine absorption model, the present fit also reproduces the Kα-to-Kβ ratio correctly. Parameter probability distributions are displayed in Fig. 5, and corresponding parameter values are listed in Table 2. Uncertainties on these parameters are obtained by numerical integration of the probability density function and thus represent the parameter range where 90% of the walkers settle.

3.5 Numerical modeling of a dusty absorber: fluorescence, line strength, and Compton shoulder

The results presented in the previous section show that a model in which the absorber purely consists of dust (i.e., a medium that is strongly enhanced in iron such as the olivine absorber used here as an example) can reproduce the observed spectral features in a quantitative way. In order to obtain a comparatively simple model that allows spectral fitting, we ignored several potentially important effects, including self-absorption of the fluorescence line emission, Compton scattering, or radiative transfer effects due to the source geometry. While the fitted iron column density and the continuum shape obtained from these fits have a direct physical meaning, these effects were modeled by the introduction of the scaling parameter η in the fit model of Eq. (5).

Having shown that the dust absorber yields a good description of the data, in this section we refine the model further and perform direct Monte Carlo simulations of photon absorption and propagation in olivine dust grains. In order to do so we utilize the SIXTE software package (Dauser et al. 2019). Although this package is designed to simulate instrument effects in X-ray detectors, it contains tools that allow direct modeling of the physics of radiation propagation in any medium, including the effects of Compton scattering off bound electrons and of fluorescence.

In a first step we apply Eq. (A1) from Wilms et al. (2000) to calculate the optical depth of individual dust grains for iron Kα fluorescence photons in order to quantitatively confirm that we can ignore self-absorption inside a single grain. We assume a grain density ρ of 1 g cm⁻³, resulting in a number density of 4.02 × 10²¹ unit cells cm⁻³. For simplicity, we used partial photoionization cross sections by Verner et al. (1996) without fine structure. For typical grain sizes of 0.025–0.25 μm (Draine & Lee 1984) we find the optical depth of the individual grains to range from 2 × 10⁻⁴ to 2 × 10⁻³. Re-absorption of the fluorescence photon only becomes dominant for grain sizes larger than ~ 1 μm, while smaller grains are mostly transparent for photon energies below the iron K-edge. We therefore conclude that the grain size distribution has little effect on the fluorescence properties.

One significant difference of dust versus gas absorption is known as self-blanketing (Fireman 1974) where atoms in the interior of the dust grain are partially shielded by atoms on the surface. Effectively, this reduces the optical depth of a medium when its constituents are confined in dust grains compared to the same number of atoms in gas phase (assuming identical chemical composition). We calculate the self-blanketing factor f for olivine grains using Eqs. (A1) and (A5) of Wilms et al. (2000). For display purposes, the fractional reduction of the optical depth, 1 − f, due to self-blanketing instead of f is shown in Fig. 6.

Again, self-blanketing becomes more significant with increasing grain size, but the change in optical depth is generally below 1% for the grain sizes and photon energies considered here and therefore well below any systematic error that is introduced by our specific choice for the dust material and the unknown size distribution.

Because of the weak dependency of the results on the grain size distribution, in order to study radiative transfer effects we next simulate the emergence of fluorescence photons from spherical and slab-like absorbers of considerable optical depth. We perform simulations for a central point source embedded in a homogeneous, quasi-spherical olivine absorber, and for an infinite slab between a point source and the observer. We explicitly calculate the absorption, fluorescence, and Compton scattering of primary and fluorescence photons in this absorber. Our implementation of Compton scattering correctly takes into account incoherent scattering from bound electrons (Hubbell et al. 1975). To estimate the scattering factor for olivine, we took the weighted mean of the constituents and obtained the energy-dependent, total Compton scattering cross section by numerical integration over all scattering angles. The incident spectrum on these clouds is fixed to a power law with exponential cut off with the best-fit parameters obtained in the previous section (Table 2), while we vary the absorber’s column density. Leahy & Creighton (1993, and references therein) and Matt (2002) present results of earlier simulations of fluorescence calculations and the Compton shoulder for the case of spherical and slab geometries, respectively, which used approximate photoabsorption cross sections and assumed free electrons.

The Monte Carlo simulations clearly show that the emerging flux of the Fe Kα fluorescence line is strongly geometry dependent (see Fig. 7). The reason is that the most probable location of the initial fluorescence event is located in the vicinity of the absorber’s τ = 1 surface. In the case of a slab geometry the fluorescence photons that are emitted into the hemisphere towards the observer will have to penetrate a higher optical depth than those emitted opposite the observers side (i.e., toward the point source), where in contrast to the K-edge energy, the optical depth at 6.4 keV is significantly less than 1. For a sphere, however, this asymmetry does not exist, and as a result the observable line flux for the slab geometry is a factor of 2–3 lower than for the spherical geometry.

To connect these results to the best-fit results from Sect. 3.4, we derive the scaling factor η included in Eq. (5) for both simulations by comparing the simulated fluorescence photon flux with that predicted from Eq. (4). Figure 7 shows η as obtained from our simulations for spherical and slab geometry as a function of column density. The parameters obtained in Sect. 3.4 are very close to those obtained for the spherical geometry. Given the large difference in the predicted values for η between the two geometries, we conclude that a dusty geometry with a quasi-spherical absorber is qualitatively consistent with the observations, while a slab geometry is unlikely to explain the data, even if some of the assumptions of our model, such as the homogeneous distribution of dust, are relaxed.

Finally, as discussed above, in addition to our estimate of the primary Fe Kα fluorescence line flux, our Monte Carlo code considers Compton scattering off the electrons in the dust grains. This allows us to quantify the formation of the Compton shoulder in an absorber that consists of dust grains. Figure 8 shows the simulated energy distribution of the escaping photons around the Fe Kα fluorescence lines with Compton continuum for a quasi-spherical, homogeneous olivine absorber with a Fe column density 6 × 10¹⁹ cm⁻². For the range of iron column densities considered in our simulations (6 × 10¹⁸ cm⁻² ≤ N_Fe ≤ 1 × 10²⁰ cm⁻²) we find Compton shoulder flux ratios of ~3 × 10⁻⁴ to 2 × 10⁻², which is several orders of magnitude lower than for gas absorbers of the same iron column density with solar abundances. Among minor effects like the modification of the Compton scattering cross section for bound electrons, this is mainly because the effective electron number density per iron atom in our pure olivine dust model is roughly a factor 300 lower than the gas of solar or ISM abundances.

Fig. 4 Most probable olivine absorption and fluorescence model. Colors and panels are defined in Fig. 1.

Fig. 5 Parameter probability distribution for the olivine absorption model including fluorescence. Additional parameters are the photon index Γ and folding energy Efold of the continuum and the flux of the Ni Kα line (not accounted for by the dust absorption model). The column density $N_{\text{H}}^{\text{Dust}}$ refers to theequivalent hydrogen column density for the iron abundance of Asplund et al. (2009).

Fig. 6 Reduction of the optical depth due to self-blanketing effects.

Fig. 7 Correction factor η for sphericaland slab geometry (red and blue, respectively) as a function of column density obtained from Monte Carlo simulations. The black shaded region shows the 90% confidence contour from the MCMC spectral modeling.

Fig. 8 Simulated Fe Kα line profile in an olivine absorber of Fe column density of 6 × 1019 cm−2, including the Compton continuum.

4 Time-resolved spectral analysis

So far, we have studied the time-averaged spectrum of IGR J16318–4848. As shown in Fig. 9a, however, the source is strongly variable on kilosecond timescales, which has been associated with changes in the absorber (Matt & Guainazzi 2003; Ibarra et al. 2007; Barragán et al. 2009). Here, we investigate the spectral variability with time.

No periodicity is apparent in the light curve, which is also confirmed by a search for periodicities on the cleaned, combined event data from FPMA and -B standard and SCIENCE_SC mode using epoch folding (Leahy et al. 1983) for trial periods of 1–1000 s, as well as a power spectrum analysis. In order to trace spectral variability over time, we use the (H−S)/(H+S) hardness ratio. This quantity shows only moderate variability (Fig. 9b). To investigate their origin, we perform time-resolved spectroscopy on segments of the light curve over which the spectral shape does not vary significantly. Because the time intervals of the constant hardness ratio are likely to represent the same spectral state, we use the hardness ratio to define these time intervals using a Bayesian block analysis. We note, however, that for the given dataset the result of the Bayesian block analysis depends significantly on the ncp_prior parameter (for details, see Scargle et al. 2013), which controls the prior for the number of blocks and was set to 0.5. With this choice of prior, we caution that the detection of change points in the hardness ratio is only marginally significant.

Considering only time intervals of simultaneous NuSTAR and XMM-Newton coverage, we obtained 24 spectra of ~0.7 to ~5.3 ks net exposure, each. The EPIC-MOS spectra have insufficient S/N for splitting them up further so for the spectral analysis we used only NuSTAR/FPMA and FPMB and XMM-Newton/EPIC-pn spectra, which we rebinned to a minimum S/N of 5 using the same channel and energy restrictions as those used in Sect. 3 for the time-averaged spectrum.

The spectra are of insufficient quality to warrant a detailed analysis based on the dust modeling discussed above. We therefore describe thespectra using the empirical cutoffpl model (Sect. 3). An initial modeling run shows that several parameters, in particular the photon index and the folding energy of the X-ray continuum did not vary significantly over the observation (Figs. 10b and c). In subsequent modeling runs we therefore forced the detector cross-calibration constants, the emission line energies, the iron abundance, the photon index (red data points in Fig. 10), and finally also the folding energy of the continuum to be the same in all time intervals (blue data points in Fig. 10), leaving only the continuum normalization, the column density, N_H and the iron lineflux as free parameters. The evolution of the continuum parameters, continuum flux, and the ratio of the Fe Kα line flux to the hard continuum flux is shown in Fig. 10. This quantity represents information that is similar to the commonly used equivalent width, but is less affected by the continuum measurement at the line position, which is particularly difficult to constrain for this source. While the continuum flux varies by more than a factor of four throughout the observation,we find that the incident continuum shape stays constant within parameter uncertainties and that N_H and the Fe Kα flux ratio vary only moderately with time. We therefore consider it very unlikely that the results of our time-averaged analysis are affected by variability effects.

Finally, Fig. 11 shows the ratio of the Fe Kα flux to the continuum flux as a function of N_H and continuum flux. We find no correlation of the Fe Kα-to-continuum flux ratio with continuum flux, and only marginal variability of the Fe Kα flux with N_H. This behavior is in agreement with our simulations where at these optical depths of the absorber only a change in the column density by a factor of a few results in notable change in the observed fluorescence photon flux. The observed variability of N_H with time is, although large in absolute numbers, only ~30%.

Fig. 9 Panel a: NuSTAR/FPMA 3–80 keV (gray) and XMM-Newton/EPIC-pn 0.2–15 keV (blue) light curve with 500 s time resolution. Panel b: XMM-Newton/EPIC-pn light curve of the iron K band. Panel c: NuSTAR/FPMA (H−S)/(H+S) hardness ratio of the 8–15 keV and 15–80 keV band. The solid line shows the result of a Bayesian block search which defines the time intervals for the time-resolved spectroscopy.

Fig. 10 Time-resolved spectroscopy of IGR J16318–4848. The figure shows the parameter evolution with time for the cutoffpl model. Panel a: equivalent hydrogen column density in units of 1022 cm−2; panels b andc: photon index and the folding energy in keV, respectively; panel d: unabsorbed continuum flux in units of 10−9 erg cm−2 s−1; and panel e: line flux of the Fe Kα line divided by the continuum flux. Red and blue data points show fits where first Γ and then Γ and Efold have been fixed, respectively.

Fig. 11 Fe Kα flux ratio to the 10–50 keV continuum flux as function of NH (panel a) and continuum flux (panel b).

5 Discussion

We have presented a detailed analysis of the absorption and fluorescence features of the HMXB IGR J16318–4848. The main focus of this work is to investigate the possible presence of a dust absorber in order to resolve open questions raised by previous work using pure neutral gas absorption.

5.1 Neutral and ionized gas absorbers

The earliest XMM-Newton observations from 2003 showed already that IGR J16318–4848 is heavily absorbed with a column density on the order of 10²⁴ cm⁻², and shows avery strong fluorescence line complex of iron and nickel. Subsequent observations were mainly modeled with an empirical cut off power-law continuum with a neutral absorber and additional Gaussian emission lines. Recently, one of the few Hitomi observations confirmed that the fluorescence lines are consistent with a neutral medium (Hitomi Collaboration 2018).

The profile of the iron Kα line is narrow, and previous authors give very low flux or only upper limits on the Compton shoulder flux, a spectral feature emerging from the reprocessing of the fluorescence photons by the Compton-thick absorber. In our empirical spectral analysis we are able to confirm this low Compton shoulder flux, the column density of ~ 2 × 10²⁴ cm⁻², and the overall continuum shape. This empiric evidence led Matt & Guainazzi (2003) to the conclusion that the absorber is anisotropic and, more importantly, cannot have a covering fraction higher than 0.2–0.3. It is important to note that these and similar analyses by other authors treated the continuum absorption and fluorescence independently.

Aiming at a more self-consistent description of the absorbing medium we then compared the data to a grid of synthetic spectra using the photoionization code XSTAR for a wide range of covering fractions, ionization parameters, and column densities. Interestingly, a single gas absorber failed to provide an adequate description of the data. The Kα-to-Kβ flux ratio and the depth and structure of the absorption edge are very robust probes of the absorbing material, unlike the continuum absorption modeling process, which is systematically affected by the choice of the empirical incident continuum.

Although our fits of a single gas absorber are statistically not acceptable, they already indicate that covering fractions very close to 100% are necessary to reproduce the correct order of magnitude of continuum absorption and fluorescence line flux. This seems intuitively reasonable since at lower covering fractions significant leakage of the incident continuum will lower the line flux with respect to the continuum flux. The spectral shape below the iron K-edge is, however, completely dominated by the fluorescence lines and we do not see a significant contribution from the continuum. At the same time, high incident radiation is needed to produce the observed fluorescence flux in the first place. We therefore conclude that the covering fraction along the line of sight, despite the low Compton shoulder flux, has to be close to unity.

5.2 Dust absorption and fluorescence

To resolve the issue of the missing Compton shoulder and the very strong Fe Kα line and edge, we propose that a significant part of the X-ray absorption is due to dust and not to gas. This idea is also motivated by near IR observations of the system which suggest the presence of a significant amount of dust in the binary system (Moon et al. 2007; Chaty & Rahoui 2012). We further take into account that the dust not only contributes to the absorption, but that the physics of photo absorption also dictates that the absorber must be responsible for a part of the strong fluorescence emission line. We therefore implement a simple dust absorption model based on laboratory measurements of olivine absorption cross sections by Rogantini et al. (2018).

In our spectral modeling, we limit ourselves to an empirical estimate of the fluorescence yield of the dust absorber. Our dust absorption and fluorescence model assumes that the iron is confined in olivine dust grains for which both the K-edge depth and continuum absorption fit the data remarkably well, even though we note that fits with other dust models that are similarly enriched in iron might work equally well. In this fit, the total iron column is mainly constrained by the Fe K-edge depth, while the strength of the Fe Kα line encodes information about the geometry and about re-absorption probabilities of emitted photons.

In order to interpret this empirical estimate of the fluorescence yield in a physical way, we performed Monte Carlo simulations of the two extreme cases of spherical and slab absorber. We find that the spherical absorber reproduces the observed line strength well, whereas in the slab geometry too much fluorescence emission is scattered out of the line of sight to the observer.

In addition to the flux of the primary iron Kα fluorescence line, the Compton shoulder flux has raised questions about the structure and geometry of the absorber. We therefore include Compton scattering in our Monte Carlo simulation in order to compare the Compton shoulder flux ratio for a spherical dust absorber to the observational values. While the detailed value of the Compton shoulder flux ratio is still debated (e.g., Barragán et al. 2009 give only an upper limit, this work and Matt & Guainazzi 2003 find ~ 10%), all previous analyses agree that it is a factor of a few lower than expected from gaseous absorbers of the appropriate hydrogen column density of ~2 × 10²⁴ cm⁻² (see, e.g., simulations by Matt 2002). Our simulations for a pure dust absorber produce a very weak Compton shoulder flux ratio that is only consistent with the upper limit by Barragán et al. (2009), but lower than the value of 11 ± 4% observed in our spectral analysis. This indicates that the simple picture of a pure dust absorber, while able to reproduce the absorption and primary fluorescence line correctly, is still incomplete as it does not provide an explanation of the flux of the Compton shoulder. It is, however, very plausible that the absorber also contains a gas component that has a minor effect on the total absorption, but would give rise to a Compton shoulder feature, although much weaker than if all the aborption were through gas. Comparing our observed Compton shoulder flux ration value to the simulations for spherical geometries by Matt (2002) we conclude that the hydrogen column density is on the order of 10²³ cm⁻² assuming cold gas, similar to the average column density proposed by Matt & Guainazzi (2003). Such a low hydrogen column density cannot explain the strong continuum and iron edge absorption that can be resolved in the picture of an additional dusty absorber.

We note that these considerations assume a cold and neutral absorbing material. Scattering off a slightly warm plasma would significantly decrease the strength of the Compton shoulder (Odaka et al. 2016). Constraining the ionization state of the absorber is, however, not possible with the observations available to date since the existing Hitomi observation puts strong constraints only on the ionization of iron.

5.3 Time-resolved spectroscopy

Our analysis of time-resolved spectra of constant hardness showed that the source mainly varies in flux on kilosecond timescales while the incident continuum shape shows no significant variation and the absorption column and the relative strength of the iron fluorescence line with respect to the hard continuum flux only vary moderately on longer timescales. Interestingly, the ratio of iron fluorescence line to continuum flux shows no distinct correlation with the absorption column or the flux. While the time-averaged spectrum is very nicely described by a homogeneous spherical dust absorber, the origin of the time variability remains an open question. Possible explanations include a structured and possibly anisotropic absorber, as suggested in previous works (e.g., Matt & Guainazzi 2003; Ibarra et al. 2007; Barragán et al. 2009) or a reprocessing structure with light-traveling times different from those probed in this observation.

6 Conclusion

We presented a detailed study of the absorption and fluorescence properties, as well as the broadband spectrum and its time-variability, of the heavily absorbed HMXB IGR J16318–4848 with data taken by NuSTAR and XMM-Newton simultaneously. We showed that a single homogeneous gas absorber is not able to reproduce the observed absorption and fluorescence properties, independent of its ionization state. We also investigated the effect of a dust absorber on the observed spectrum, finding that a spherical homogeneous absorber of iron column density N_Fe = ~6–7 × 10¹⁹ cm⁻² is consistent with the observations. Our time-resolved spectroscopy, however, indicates the presence of structures in the absorbing material since the absorption appears to be more variable than the incident spectrum.

In order to improve our understanding of HMXBs containing a significant amount of dust like IGR J16318–4848 we hope that our results will foster the development of more sophisticated dust absorption models that also include fluorescence processes. There are currently a few very advanced dust absorption models (e.g., GGADT; Hoffman & Draine 2016) used to calculate differential and total X-ray scattering cross sections for a specified grain using Anomalous Diffraction Theory, or xscat (Smith et al. 2016), which computes scattering cross sections for a population of dust grains as a function of the distance from the source and extraction region. However, neither of these advanced models includes treatment of the fluorescence from dust grains.

We studied the formation of a Compton shoulder with Monte Carlo simulations that are based on the best fit of our dust absorption model. The Compton shoulder flux produced by a pure dust absorber is even lower than observed in many of the earlier studies. A future step in the absorption modeling of this source will be the inclusion of an additional gas and self-consistent Compton shoulder component into the dust absorption model which is, however, beyond the scope of this work.

From an observational point of view, the existing Hitomi observation of IGR J16318–4848 is of exceptional value because it resolved the Fe Kα line with unprecedented quality and put strong constraints on the ionization state of the absorber. We therefore expect high energy resolution instruments like the Resolve micro-calorimeter on board XRISM (Tashiro et al. 2018) to provide further insight into dust geometry and composition. Of particular interest will be any measurement of the X-ray absorption fine structure of the iron K-edge because it is a major diagnostic of the elemental composition and structure of the dust, which so far has been mainly studied through gratings observations in the L-edge band (e.g., Westphal et al. 2019) at an energy that is not accessible for IGR J16318–4848.

Acknowledgements

We thank the anonymous referee for constructive comments that helped to improve the manuscript. We thank Jakob Stierhof, Daniel Stern, and Stefan Licklederer for extensive and fruitful discussions. R.B. acknowledges partial funding by the Bundesministerium für Wirtschaft under Deutsches Zentrum für Luft- und Raumfahrt contract 50 OR 1410. V.G. is supported through the Margarete von Wrangell fellowship by the ESF and the Ministry of Science, Research and the Arts Baden-Württemberg. Work at LLNL is conducted under the auspices of the US DOE under Contract DE-AC52-07NA27344. This research is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA and has made use of data obtained with NuSTAR, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by NASA. We acknowledge the 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). This work has made use of a collection of ISIS functions (ISISscripts) provided by ECAP/Remeis observatory and MIT (http://www.sternwarte.uni-erlangen.de/isis/).

References

All Tables

All Figures

	Fig. 1 Panel a: time-averaged counts spectrum of IGR J16318−4848 with best-fit cutoff power law model for NuSTAR/FPMA (gray), FPMB (gold), XMM-Newton/EPIC-pn (blue), and XMM-Newton/EPIC-MOS1 and -MOS2 (dark and light green, respectively). Panels b and c: residuals and ratio for the best-fit model.
In the text

	Fig. 2 Simulated XSTAR spectra for a wider range of covering fractions. Other model parameters are NH = 1024 cm−2 and log ξ = 1. Vertical arrows indicate the maximum line flux for clarity.
In the text

	Fig. 3 Best fit of the single XSTAR model. Colors and panels are defined in Fig. 1.
In the text

	Fig. 4 Most probable olivine absorption and fluorescence model. Colors and panels are defined in Fig. 1.
In the text

Fig. 5 Parameter probability distribution for the olivine absorption model including fluorescence. Additional parameters are the photon index Γ and folding energy Efold of the continuum and the flux of the Ni Kα line (not accounted for by the dust absorption model). The column density $N_{\text{H}}^{\text{Dust}}$ refers to theequivalent hydrogen column density for the iron abundance of Asplund et al. (2009).

In the text

	Fig. 6 Reduction of the optical depth due to self-blanketing effects.
In the text

	Fig. 7 Correction factor η for sphericaland slab geometry (red and blue, respectively) as a function of column density obtained from Monte Carlo simulations. The black shaded region shows the 90% confidence contour from the MCMC spectral modeling.
In the text

	Fig. 8 Simulated Fe Kα line profile in an olivine absorber of Fe column density of 6 × 1019 cm−2, including the Compton continuum.
In the text

	Fig. 9 Panel a: NuSTAR/FPMA 3–80 keV (gray) and XMM-Newton/EPIC-pn 0.2–15 keV (blue) light curve with 500 s time resolution. Panel b: XMM-Newton/EPIC-pn light curve of the iron K band. Panel c: NuSTAR/FPMA (H−S)/(H+S) hardness ratio of the 8–15 keV and 15–80 keV band. The solid line shows the result of a Bayesian block search which defines the time intervals for the time-resolved spectroscopy.
In the text

Fig. 10 Time-resolved spectroscopy of IGR J16318–4848. The figure shows the parameter evolution with time for the cutoffpl model. Panel a: equivalent hydrogen column density in units of 1022 cm−2; panels b andc: photon index and the folding energy in keV, respectively; panel d: unabsorbed continuum flux in units of 10−9 erg cm−2 s−1; and panel e: line flux of the Fe Kα line divided by the continuum flux. Red and blue data points show fits where first Γ and then Γ and Efold have been fixed, respectively.

In the text

	Fig. 11 Fe Kα flux ratio to the 10–50 keV continuum flux as function of NH (panel a) and continuum flux (panel b).
In the text