Monte Carlo simulations of the detailed iron absorption line profiles from thermal winds in X-ray binaries

Monte Carlo simulations of the detailed iron absorption line profiles from thermal winds in X-ray... Abstract Blueshifted absorption lines from highly ionized iron are seen in some high inclination X-ray binary systems, indicating the presence of an equatorial disc wind. This launch mechanism is under debate, but thermal driving should be ubiquitous. X-ray irradiation from the central source heats disc surface, forming a wind from the outer disc where the local escape velocity is lower than the sound speed. The mass-loss rate from each part of the disc is determined by the luminosity and spectral shape of the central source. We use these together with an assumed density and velocity structure of the wind to predict the column density and ionization state, then combine this with a Monte Carlo radiation transfer to predict the detailed shape of the absorption (and emission) line profiles. We test this on the persistent wind seen in the bright neutron star binary GX 13+1, with luminosity L/LEdd ∼ 0.5. We approximately include the effect of radiation pressure because of high luminosity, and compute line features. We compare these to the highest resolution data, the Chandra third-order grating spectra, which we show here for the first time. This is the first physical model for the wind in this system, and it succeeds in reproducing many of the features seen in the data, showing that the wind in GX13+1 is most likely a thermal-radiation driven wind. This approach, combined with better streamline structures derived from full radiation hydrodynamic simulations, will allow future calorimeter data to explore the detail wind structure. accretion, accretion discs, black hole physics, X-rays: binaries, X-rays: individual: (GX 13+1) 1 INTRODUCTION Low-mass X-ray binaries (LMXBs) are systems where a companion star overflows its Roche lobe, so that material spirals down towards a compact object which can be either a black hole or neutron star. The observed X-ray emission is powered by the enormous gravitational potential energy released by this material as it falls inwards, lighting up the regions of intense space–time curvature and giving observational tests of strong gravity. This inflow also powers outflows. Accretion disc winds are seen via blueshifted absorption lines from highly ionized material in high inclination LMXBs (see e.g. the reviews by Ponti et al. 2012; Díaz Trigo & Boirin 2016) but the driving mechanism of these winds is not well understood. The potential candidates are acceleration of gas by the Lorentz force from magnetic fields threading discs (magnetic driving: Blandford & Payne 1982; Fukumura et al. 2014), radiation pressure on the electrons (continuum driving: Proga & Kallman 2002; Hashizume et al. 2015) and thermal expansion of the hot disc atmosphere heated by the central X-ray source, which makes a wind at radii which are large enough for the sound speed to exceed the local escape velocity (thermal driving: Begelman, McKee & Shields 1983; Woods et al. 1996, hereafter W96). Recent work has focused on magnetic driving, primarily because of a single observation of dramatic wind absorption seen from the black hole GRO J1655−40 at low luminosity, far below the Eddington limit where continuum driving becomes important, and with a derived launch radius which is far too small for thermal driving (Miller et al. 2006; Luketic et al. 2010; Higginbottom & Proga 2015). Magnetic wind models can fit this spectrum (Fukumura et al. 2017), but then it is difficult to understand why such high column and (relatively) low ionization winds are not seen from other observations of this object or any other high inclination systems with similar luminosities and spectra. Instead, the unique properties of this wind could potentially be explained if the outflow has become optically thick along the line of sight. This would suppress the observed flux from an intrinsically super-Eddington source (Shidatsu, Done & Ueda 2016; Neilsen et al. 2016). Alternatively, this singular wind may be a transient phenomena, not representative of the somewhat lower column/higher ionization winds which are normally seen. What then powers the more typical winds? Thermal driving qualitatively fits the observed properties, as these winds are preferentially seen in systems with larger discs (Díaz Trigo & Boirin 2016). X-rays from near the compact objects (the inner disc emission and any corona and/or boundary layer) irradiate the outer disc, and the balance between Compton heating and cooling heats the surface to the Compton temperature, which is a luminosity weighted mean energy given by $$T_{\mathrm{IC}} = \frac{1}{4}\int E L(E){\rm d}E/\int L(E) {\rm d}E$$ (Begelman, McKee & Shields 1983; Done, Tomaru & Takahashi 2018, hereafter D18). This temperature is constant with radius, as it depends only on the spectrum of the radiation from the central region, but gravity decreases with distance from the source. For a large enough disc, the isothermal sound speed from this Compton temperature is bigger than the escape velocity, so the heated material makes a transition from forming a bound atmosphere, to an outflowing wind. This defines the Compton radius (RIC = GMmp/kTIC = (6.4 × 104/TIC, 8)Rg, where Rg = GM/c2 and TIC, 8 = TIC/108 K) as the typical launch radius for a thermal wind (Begelman et al. 1983; D18). These thermal wind models give an analytic solution for the mass-loss rate from the disc (Begelman et al. 1983; W96). However, to calculate the observables such as column density and ionization structure requires an understanding of the velocity as a function of two-dimensional position (or equivalently, velocity as a function of length along a streamline, together with the shape of the streamlines). D18 assume a very simple two-dimensional density and velocity structure along radial streamlines, and show that this can match the column seen in the full hydrodynamic simulation of W96. This composite model was applied to multiple spectra of H1743−322, matching the wind seen in its disc dominated state, and predicting the disappearance of these wind absorption lines in a bright low/hard state, as observed (Shidatsu & Done, in preparation). This wind disappearance does not occur from simply the change in photoionization state from the changing illumination (Miller et al. 2012). The key difference is that thermal winds respond to changing illumination by changing their launch radius, density and velocity as well as responding to the changing photoionizing spectrum (Shidatsu & Done, in preparation). Here we explore the thermal wind predictions in more detail, using a three-dimensional Monte Carlo simulation code monaco (Odaka et al. 2011) to calculate the radiation transport through the material so as to predict the resulting emission and absorption line profiles. We calculate these for the very simple constant velocity structure assumed in D18, and then extend this to consider an accelerating wind along biconical streamlines, as a more appropriate disc wind geometry (Waters & Proga 2012). We show results from both these geometries for the specific case of L/LEdd = 0.3, TIC, 8 = 0.13, Rout = 5RIC to compare with W96. We apply the biconical wind model to the bright neutron star binary system GX13+1. This system is unique amongst all the black hole and neutron star binaries in showing persistently strong blueshifted absorption features in its spectrum (Ueda et al. 2004; D'Aì et al. 2014). The X-ray continuum is rather stable, with TIC, 8 ∼ 0.13, making it a good match to the simulations, but it has slightly higher luminosity at L/LEdd = 0.5. This means that radiation pressure should become important, decreasing the radius from which the wind can be launched and increasing its column density. We compare results from a hybrid thermal–radiative wind (calculated using the approximate radiation pressure correction from D18) to the detailed absorption line profiles seen in the third-order Chandra grating data from this source. This is the first quantitative, physical model for the wind, and the first exploration of the highest spectral resolution data from this source. While magnetic driving models are always possible, our results match the majority of the observed features, showing that the wind properties are broadly consistent with hybrid thermal–radiative driving. 2 RADIATIVE TRANSFER CODE We use the Monte Carlo simulation code monaco (Odaka et al. 2011) to calculate radiative transfer through the wind. monaco uses its original physics implementation of photon interactions (Watanabe et al. 2006; Odaka et al. 2011) while it employs the Geant4 toolkit library (Agostinelli et al. 2003) for photon tracking in an arbitrary three-dimensional geometry. We consider an azimuthally symmetric density and velocity field, then use the xstar photoionization code (Kallman & Bautista 2001) to calculate the equilibrium population of ions from each element assuming one-dimensional radiation transfer from the central source. We grid in radial distance, r, and θ and assume the illuminating flux is the transmitted part of the central spectrum along the line of sight to the element. This gives an ionization parameter   \begin{equation} \xi (r, \theta ) = \frac{L \exp (-\tau _{{\rm abs}})}{n (r, \theta ) r^2} \end{equation} (1)where τabs is the optical depth to absorption (but does not include electron scattering), and n(r, θ) is the density as a function of position. We then use monaco to track photons through this ionization structure, including their interaction with this material. Photons interacting with ions can be absorbed in photoionization or photoexcitation, and photons generated via recombination and atomic de-excitation are tracked. Doppler shifts of the absorption cross-sections from the velocity structure of the material are included, as is the Compton energy change on interaction with free electrons. Ideally, this calculated radiation field should then be used as input to xstar, the ion populations recalculated, and the process should be iterated until convergence. However, for our simulations here the wind is mostly optically thin, so we do not include this self-consistent iteration. This method of radiation transfer calculation is based on the Hagino et al. (2015), but the geometry and the velocity/density structure we use is reflected on the thermally driven winds in LMXBs whereas Hagino et al. (2015) is focused on the UV-line driven winds in active galactic nucleus. Also, as thermal winds are typically highly ionized, we consider only H-like and He-like ions of Fe and calculate the spectrum only over a restricted energy band of 6.5–7.2 keV. Table 1 details the transitions used. Table 1. Detailed parameters for each line included in these monaco simulations. Note that we list only lines which have larger oscillator strengths than 10−3. Line ID  Energy (keV)  Oscillator strength  Fe XXV Heαy  6.668  6.57 × 10−2  Fe XXV Heαw  6.700  7.26 × 10−1  Fe XXVI Lyα2  6.952  1.36 × 10−1  Fe XXVI Lyα1  6.973  2.73 × 10−1  Line ID  Energy (keV)  Oscillator strength  Fe XXV Heαy  6.668  6.57 × 10−2  Fe XXV Heαw  6.700  7.26 × 10−1  Fe XXVI Lyα2  6.952  1.36 × 10−1  Fe XXVI Lyα1  6.973  2.73 × 10−1  View Large 3 RADIAL STREAMLINES: D18 3.1 Geometry and parameters We first consider the radial streamline wind model of D18. This calculates the analytic mass-loss rate per unit area, $$\dot{m}(R)$$ where R denotes distance along the disc plane. Integrating over the whole disc gives the total mass-loss rate in the wind, $$\dot{M}$$. This is assumed to flow along radial (centred at the origin) streamlines from a launch radius which is 0.2RIC for high L/LEdd, with constant velocity set at the mass-loss weighted average escape velocity. The mass-loss rate along each radial streamline is weighted with angle such that $$\dot{M}(\theta )\propto \dot{M}(1-\cos \theta )$$, and then mass conservation gives n(r, θ) ∝ (1 − cos θ)/r2. D18 show that these assumptions lead to a total column density through the structure which matches to within a factor 2 of that in the hydrodynamic simulations of W96 (see also Section 4). We put this structure into monaco for L = 0.3LEdd with TIC = 1.3 × 107 K and Rout = 5RIC (mass-loss rate $$\dot{M}_w = 2.0 \times 10^{19}\,\mathrm{g\, s^{-1}}$$ for a 10 M⊙ black hole which means the ratio of mass-loss rate to mass accretion rate $$\dot{M}_w/\dot{M}_a = 3.9$$, where the $$\dot{M}_a = L/(0.1c^2)$$, launch radius of 0.2RIC ≈ 105Rg and weighted average vout = 420 km s−1). We include turbulence, assuming vturb = vout, and calculate the rotation velocity along each stream line assuming angular momentum conservation (see Appendix A). We make a grid which follows the symmetry of the assumed structure, i.e. centred on the origin, with 20 linearly spaced spherical shells from 0.2 to 5RIC, and 20 angles, linearly spaced in θ from 7° to 83° (see below). This density structure is shown in the left-hand panel of Fig. 1, while the right-hand panel shows the mean Fe ion state obtained from the xstar calculation. This is constant along each streamline because the ionization parameter ξ = L/(nr2), and the constant velocity radial streamlines mean that density decreases as 1/r2. Fe is almost completely ionized over the whole grid, with a small fraction of hydrogen-like iron remaining only for high inclination streamlines. Figure 1. View largeDownload slide Distribution of density (left) and mean Fe ionization state (right) for the radial streamline model. Figure 1. View largeDownload slide Distribution of density (left) and mean Fe ionization state (right) for the radial streamline model. 3.2 monaco output Fig. 2 shows the resulting spectra at three different inclination angles. These show that the emission lines are always similarly weak, and that the electron scattered continuum flux makes only a ∼0.5 per cent contribution to the total flux, but that the absorption lines strongly increase at higher inclination angles. We calculate the equivalent width (EW) of each emission and absorption line by fitting the continuum outside the emission and absorption regions with an arbitrary function [F(E) = aEbE + c; Odaka et al. 2016]. The EW of each emission and absorption line is then measured by numerical integration of the difference between the model and the simulation data. The left-hand panel of Fig. 3 shows the EW of the He-like (red) and H-like Lyα2 (green) and Lyα1 (blue) absorption lines. The corresponding emission lines always have EW lower than 0.1 eV so are not seen on this plot. The strong increase of the absorption line EW with inclination clearly shows that the wind is equatorial (by construction from the 1 − cos θ density dependence and constant velocity assumptions). At inclinations above 70°, the Doppler wings of the Kα1 and Kα2 absorption lines merge together due to the turbulent velocities, so Fig. 3 shows only a single EW for this blend. Figure 2. View largeDownload slide Spectra computed for the radial streamline model with 1 eV resolution for different lines of sight. Each panel shows the spectrum in a different inclination angle bin (the angular bin sizes are indicated the top of each panel). The total spectrum is shown in black (top), the spectrum direct photons in red and scattered/reprocessed spectrum is blue (bottom). Note that the vertical axis is plotted linearly in the top panels but logarithmically in the bottom panels. Lines are Fe XXV (6.668 keV for Heαy and 6.700 keV for Heαw) and Fe XXVI (6.952 keV for Lyα2 and 6.973 keV for Lyα1). The equatorial density structure of the wind means that the absorption is much stronger at high inclination angles. The emission is more isotropic, so it can clearly be seen at low inclination angles, but is absorbed by the wind at high inclinations. Figure 2. View largeDownload slide Spectra computed for the radial streamline model with 1 eV resolution for different lines of sight. Each panel shows the spectrum in a different inclination angle bin (the angular bin sizes are indicated the top of each panel). The total spectrum is shown in black (top), the spectrum direct photons in red and scattered/reprocessed spectrum is blue (bottom). Note that the vertical axis is plotted linearly in the top panels but logarithmically in the bottom panels. Lines are Fe XXV (6.668 keV for Heαy and 6.700 keV for Heαw) and Fe XXVI (6.952 keV for Lyα2 and 6.973 keV for Lyα1). The equatorial density structure of the wind means that the absorption is much stronger at high inclination angles. The emission is more isotropic, so it can clearly be seen at low inclination angles, but is absorbed by the wind at high inclinations. Figure 3. View largeDownload slide Left-hand panel: EW of the absorption lines as a function of inclination angle, Fe XXV (Heαw, red) and Fe XXVI (Lyα2, green and Lyα1, blue). The EW of all absorption lines increases strongly at higher inclination, showing the assumed equatorial disc wind geometry. The Doppler wings (with width set by turbulent velocity) of the two H-like absorption lines start to merge for inclinations above 70° so above this we show the total EW of the two lines. Right-hand panel: the blueshifted absorption line velocity for each ion species. This clearly shows the assumed constant velocity structure of the radial streamline. Figure 3. View largeDownload slide Left-hand panel: EW of the absorption lines as a function of inclination angle, Fe XXV (Heαw, red) and Fe XXVI (Lyα2, green and Lyα1, blue). The EW of all absorption lines increases strongly at higher inclination, showing the assumed equatorial disc wind geometry. The Doppler wings (with width set by turbulent velocity) of the two H-like absorption lines start to merge for inclinations above 70° so above this we show the total EW of the two lines. Right-hand panel: the blueshifted absorption line velocity for each ion species. This clearly shows the assumed constant velocity structure of the radial streamline. Fig. 3 (right) shows the outflow velocity, as measured from the energy of the deepest absorption lines (with error set by the resolution of the simulation to ±0.5 eV). These velocities are constant within 25  per cent as a function of inclination, again by construction due to the assumption of constant radial velocity along radial streamlines. 4 DIVERGING WIND In Section 3, we considered a wind model with constant velocity along radial streamlines. However, the expected thermal wind geometry is instead much more like an accelerating, diverging biconical wind Waters & Proga (2012). Full streamline structures which give the density and velocity of the wind at all points can only be found by hydrodynamic calculations (but see Clarke & Alexander 2016 for some analytic approximations). Since modern calculations only exist for the singular case of GRO J1655−40, we follow D18 and use the W96 simulation results. W96 does not give full density/velocity structures, but do give total column density through the wind at three different luminosities. We use these to match to our assumed streamline and velocity structure, which is the standard biconical diverging disc wind used in a variety of systems including cataclysmic variables (Knigge, Woods & Drew 1995; Long & Knigge 2002) and Active Galaxies (Sim et al. 2010; Hagino et al. 2015). 4.1 Geometry and parameters The geometry can be defined by three parameters (Fig. 4). Rin = 0.1RIC, the distance from the source to the inner edge of the wind Rout, the distance from the source to outer edge of the wind αmin, the angle from z-axis to the inner edge of the wind Figure 4. View largeDownload slide The geometry of the diverging biconical wind model. Figure 4. View largeDownload slide The geometry of the diverging biconical wind model. The disc wind is fan-shaped, with a focal point offset down from the centre by a distance d = 0.1RIC/tan αmin so that the wind fills the angles from αmax − αmin down to the disc surface. We use R to denote distance along the disc surface, and r, θ denote radial distance and polar angle from the origin, as before. αmin (or equivalently d) is a free parameter, which sets the wind geometry. Streamlines are assumed to be along lines of constant angle α (where αmin < α < αmax) from the focal point. Distance along a streamline which starts on the disc at radius R is l(R) (see Appendix. A). Velocity along the streamline is assumed to be of the form $$v (r, \theta )=f_v c_{{\rm ch}}(r)\sqrt{\frac{l(r, \theta )}{R(r)}}$$, i.e. this wind accelerates with distance along the streamline, with a terminal velocity which is related via a free parameter fv to the characteristic sound speed cch, given by the balance between heating and cooling in the time it takes the wind to reach a height H ∼ R (D18). The density structure is solved by the mass conservation continuity equation along streamlines (see Appendix B). We calculate the wind properties out to a distance which is twice that of the focal point of the wind to Rout. We set the free parameter values, fv and αmin, and calculate the total column along each line of sight, NH(θ), to the central source for parameters matching to the three W96 simulations. These are L/LEdd = 0.3, 0.08 with Rout = 5RIC and L/LEdd = 0.01 with Rout = 12RIC. We adjust fv and αmin to minimize the difference between our model and W96. We find αmin = 7° and fv = 0.25 matches within a factor 2 of the results from W96. Fig. 5 shows results with these parameters (filled circles), compared to the radial wind model of Section 3 (open circles) as well as the W96 results (solid line). This more physically realistic geometry and velocity gives a similarly good match to the simulations as the D18 radial wind. Figure 5. View largeDownload slide The solid lines show column density as a function of the cosine of the inclination angle through the wind resulting from the hydrodynamic simulations of W96 for L/LEdd = 0.01(green), 0.08 (red), 0.3 (black). The filled circles show that resulting from the diverging biconical wind (Section 4) while the open circles show the radial streamline model of D18 (Section 3). Figure 5. View largeDownload slide The solid lines show column density as a function of the cosine of the inclination angle through the wind resulting from the hydrodynamic simulations of W96 for L/LEdd = 0.01(green), 0.08 (red), 0.3 (black). The filled circles show that resulting from the diverging biconical wind (Section 4) while the open circles show the radial streamline model of D18 (Section 3). The resulting density structure from this different geometry and velocity are shown in the left-hand panel of Fig. 6 for L/LEdd = 0.3. Comparing this with the radial wind shows that the density is higher closer to the disc, and lower further away due to the material accelerating away from the disc rather than being at constant velocity. We run xstar as before, and the right-hand panel of Fig. 6 shows that this leads to a lower mean ionization state of Fe than before, and this is no longer constant along the radial sightline due to the different wind geometry. Figure 6. View largeDownload slide Distribution of density (left) and Fe ionization state (right) for the diverging wind geometry. The accelerating flow gives higher density material close to the disc compared to the constant velocity outflow model in Fig. 1, giving lower mean ionization state. Figure 6. View largeDownload slide Distribution of density (left) and Fe ionization state (right) for the diverging wind geometry. The accelerating flow gives higher density material close to the disc compared to the constant velocity outflow model in Fig. 1, giving lower mean ionization state. 4.2 monaco output We calculate the emission and absorption lines resulting from the different wind structure (Fig. 7). The diverging bipolar wind has higher density material closer to the source compared to the radial wind geometry, so it subtends a larger solid angle to scattering. This means that there is more emission line contribution, as well as a higher fraction of electron scattered continuum (around 2 per cent, see the lower panel of Fig. 7). The left-hand panel of Fig. 8 shows the emission line EW (dotted lines) for each ion species (red: Fe XXV w, green: Fe XXVI Lyα2, blue: Lyα1). These can now be of order 1 eV for face on inclinations, decreasing at higher inclination as they are significantly suppressed by line absorption. Figure 7. View largeDownload slide As in Fig 2, but for the diverging biconical wind geometry. Figure 7. View largeDownload slide As in Fig 2, but for the diverging biconical wind geometry. Figure 8. View largeDownload slide As in Fig. 3, but for the diverging wind model. The lower ionization state means that there is also a contribution from the intercombination line of Fe XXV Heαy (black) at the highest inclinations. Figure 8. View largeDownload slide As in Fig. 3, but for the diverging wind model. The lower ionization state means that there is also a contribution from the intercombination line of Fe XXV Heαy (black) at the highest inclinations. The corresponding absorption line EWs are shown as the solid lines (compare to Fig. 3). The lower mean ionization state leads to more He-like Fe, so there is more of this ion seen in absorption than in the radial streamline model. These absorption lines increase as a function of inclination angle as before, but now the Lyα1 and Lyα2 do not merge together at the highest inclination angles due to the different velocity structure (see right-hand panel of Fig. 8). The lines are formed preferentially in the higher density material close to the disc. The assumed acceleration law means that the typical velocities here are lower than in the constant velocity model, as the material has only just begun to accelerate. Thus, the turbulence is also lower, so the Doppler width of the absorption lines is smaller. This also means that the absorption line saturates to a constant EW at lower column density, so the EW of the absorption lines does not increase so strongly as before at the highest inclination angles. 5 COMPARISON WITH GX13+1 We now use the more physically motivated diverging biconical wind geometry to compare with observational data. An ideal source would be one which is not too different from the parameters simulated in the previous sections, as here we know the total column from W96 and know that our assumed velocity/density matches to this. Of the sources listed in Díaz Trigo & Boirin (2016), the neutron star LMXB GX13+1 is the source which has most similar L/LEdd and TIC to that assumed here, and it also has the advantage that it is a persistent source, with relatively constant luminosity and spectral shape, and it shows similarly strong absorption lines in multiple data sets. 5.1 Observational data GX13+1 was observed by the Chandra HighEnergy Transmission Grating (HETG) four times in 2 weeks in 2010 (Table 2). The first-order data are shown in D'Aì et al. (2014) and reveal multiple absorption features from highly ionized elements (see also Ueda et al. 2004 for similar features in an earlier observation). Higher order grating spectra give higher resolution, as demonstrated for the black hole binaries by (Miller et al. 2015). Here, we show for the first time the third-order Chandra data for GX13+1. We extract first- and third-order HEG spectra from these observations, using ciao version 4.9 and corresponding calibration files. We reprocess the event files with $$`\mathrm{chandra\_repro}$$’, and make response files using ‘mktgres’ to make the redistribution and ancillary response files. We run ‘tgsplit’ to get the HEG ±3 spectra, and run ‘combine_grating_spectra’ to combine HEG plus and minus orders for each observation to derive a single first-order spectrum (black), and a single third-order spectrum (red) as shown in Fig. 9. The first-order spectra can resolve the components of the He-like Fe triplet, with a clear dip to the low energy side at the resonance line energy of 6.7 keV, but the H-like Lyα1 and α2 are blended together. The higher resolution of the third-order spectra is able to clearly separate the He-like intercombination and resonance lines, and even the H-like Lyα1 and α2 (Miller et al. 2015). Figure 9. View largeDownload slide HEG spectra of GX 13+1 from first order (black) and third order (red). Figure 9. View largeDownload slide HEG spectra of GX 13+1 from first order (black) and third order (red). Table 2. List of the Chandra HETG observations. OBSID  MODE  Date  Exposure (ks)  11815  TE-F  24/07/2010  28  11816  TE-F  30/07/2010  28  11814  TE-F  01/08/2010  28  11817  TE-F  03/08/2010  28  OBSID  MODE  Date  Exposure (ks)  11815  TE-F  24/07/2010  28  11816  TE-F  30/07/2010  28  11814  TE-F  01/08/2010  28  11817  TE-F  03/08/2010  28  View Large 5.2 Model of GX 13+1 We fit the contemporaneous RXTE spectrum (ObsID 95338-01-01-05) with a model consisting of a disc, Comptonized boundary layer and its reflection. The resulting inverse Compton temperature of the continuum (disc plus Comptonization) is TIC ∼ 1.2 × 107 K, almost identical to the simulation (see also D'Aì et al. 2014). The luminosity is L = 0.5LEdd (Díaz Trigo, Migliari & Guainazzi 2014; D'Aì et al. 2014), similar to the maximum simulation value of L = 0.3LEdd in W96. The simulation also requires Rout, which can be calculated from the orbital period and mass of binary stars. GX 13+1 has 24 d orbital period, and the neutron star and companion have masses of 1.4 and 5 M⊙, respectively (Bandyopadhyay et al. 1999; Corbet et al. 2010). This gives a binary separation a = 4.6 × 1012 cm, for a Roche lobe radius RR/a = 0.27. The disc size is then Rout = 10RIC assuming that Rout = 0.8RR (Shahbaz, Charles & King 1998), double the value assumed in the simulations. D18 shows that this increase in disc size makes the predicted column slightly larger, but the effect is fairly small (Fig. 3: D18). Fig. 10 (blue line) shows the predicted column density through the wind as a function of inclination angle. This is very similar to the column predicted for the fiducial simulations (Fig. 5) Figure 10. View largeDownload slide The column density as a function of the cosine of the inclination angle for the diverging biconical wind calculated for the system parameters of GX13+1. The blue line shows the predictions for a purely thermal wind, while the red includes a very simple treatment of radiation pressure. The source has L/LEdd ∼ 0.5, so the thermal wind can be launched from closer in due to the lower effective gravity. This effect has a large impact on the predicted column, so the details of how this radiation pressure correction affects the velocity and density structure will be important in determining the line profiles. Figure 10. View largeDownload slide The column density as a function of the cosine of the inclination angle for the diverging biconical wind calculated for the system parameters of GX13+1. The blue line shows the predictions for a purely thermal wind, while the red includes a very simple treatment of radiation pressure. The source has L/LEdd ∼ 0.5, so the thermal wind can be launched from closer in due to the lower effective gravity. This effect has a large impact on the predicted column, so the details of how this radiation pressure correction affects the velocity and density structure will be important in determining the line profiles. However, D18 show that radiation pressure should make a rapidly increasing contribution to the wind as L/LEdd increases from 0.3 to 0.7. The GX13+1 luminosity is midway between these two, so radiation pressure should significantly lower the effective gravity, meaning that the wind can be launched from smaller radii. We follow D18 and estimate a radiation pressure correction to the launch radius of $$\bar{R}_{\mathrm{IC}} = (1.0-0.5L_{\mathrm{Edd}}/0.71L_{\mathrm{Edd}})R_{\mathrm{IC}}=0.30R_{\mathrm{IC}}$$, hence $$R_{\mathrm{out}}=33 \bar{R}_{\mathrm{IC}}$$, dramatically larger than assumed in the fiducial simulations. This correction predicts a density which is 11 times larger and column along any sightline which is 3.3 times larger assuming (as in D18) that the velocity structure is unchanged (red line, Fig. 10). This increase in Rout in terms of RIC means that more wind is produced (as in D18), so the wind efficiency increases to 4.0 (from 2.3). The column density goes close to 1024 cm−2 at high inclinations, so electron scattering becomes important. This effect reduces the illuminating ionizing flux by $$e^{-\tau _T}$$ from the central source along the line of sight to each wind element, and also increases the contribution of diffuse and scattered emission from the wind to the ionizing continuum. We include scattering, reducing the xstar illumination by $$e^{-\tau _{T}}$$ along each line of sight, but do not include the diffuse emission as the time-scale to integrate over the entire wind at each point is prohibitive. We run monaco on this wind structure to predict the detailed absorption line profiles for comparison to the third-order HEG data. Fig. 11 shows the result assuming an inclination angle of 80° (Díaz Trigo et al. 2012) which gives the best fit to the data. This gives a fairly good match to the overall absorption, except for the highest velocity material seen in the data. Lower inclination angles give higher blueshift, but lower absorption line EW, while higher inclination gives larger absorption line but lower blueshift (see Fig. 12). Thus, it is not possible to completely reproduce the observed lines in GX13+1 with our simple radiation pressure corrected thermal wind model. This is not surprising, as radiation pressure will almost certainly change the velocity law by radiative acceleration as well as changing the launch radius. Full radiation hydrodynamic simulations are required to predict the resulting velocity and density structure. None the less, our result demonstrates for the first time that hybrid thermal–radiative wind models can give a good overall match to the column and ionization state of the wind in GX13+1, and that current data can already give constraints on the velocity and density structure of this material. Figure 11. View largeDownload slide The model (red) and HEG third-order spectrum (black). The best-fitting inclination angle is i = 80°. This gives roughly the correct column of Fe XXV and XXVI at low velocity, but fails to match the observed higher velocity blue wing to the absorption features. Figure 11. View largeDownload slide The model (red) and HEG third-order spectrum (black). The best-fitting inclination angle is i = 80°. This gives roughly the correct column of Fe XXV and XXVI at low velocity, but fails to match the observed higher velocity blue wing to the absorption features. Figure 12. View largeDownload slide As in Fig. 8, but with the system parameters of GX 13+1 and the simplest radiation pressure correction to make a hybrid thermal–radiative wind. Figure 12. View largeDownload slide As in Fig. 8, but with the system parameters of GX 13+1 and the simplest radiation pressure correction to make a hybrid thermal–radiative wind. 6 DISCUSSION AND SUMMARY We construct a Monte Carlo code to calculate detailed spectra from any given density and velocity distribution of highly ionized material. We use this to explore the absorption and emission lines of H and He-like Fe for the mass-loss rates predicted from thermal wind models. We first use the radial streamline, constant velocity model of D18 which is able to reproduce the column derived from the hydrodynamic calculations of W96, but then extend this to a more realistic disc-wind geometry with gas accelerating along diverging streamlines, again reproducing the column from W96. The different assumed velocity and density structures for the thermal wind mass-loss rates give different predictions for the overall ionization state of the material, the resulting EW of emission and absorption lines, and their velocity shift. These show the potential of observations to test the detailed structure of the wind. We apply the biconical disc wind model to some of the best data on winds from an LMXB. The neutron star GX13+1 shows strong and persistent absorption features in Chandra first-order HETG spectra (Ueda et al. 2004; D'Aì et al. 2014), but here we show for the first time the higher resolution third-order data. We find that while the source is fairly well matched to the parameters of the brightest fiducial simulation in terms of TIC, the higher luminosity (L/LEdd = 0.5 compared to 0.3 for the simulation) makes a significant impact on the predicted wind properties as it puts the source firmly into the regime where radiation pressure driving should become important. We use the simple radiation pressure correction suggested by D18 and calculate the line profiles from a hybrid thermal–radiative wind. The additional radiation pressure driving means that the wind can be launched from much closer to the central source, and has higher mass-loss rate. This is the first detailed test of the absorption line profiles predicted by physical wind models on any source other than the singular wind seen in GRO J1655−40 (Luketic et al. 2010). Our simulations quantitatively match many of the observed features except for the highest velocity material. This is not surprising, given the simplistic assumptions about the effect of radiation pressure. In future work, we will use the wind velocity and density structure determined from full radiation hydrodynamics simulations in order to properly test the thermal–radiative wind models in GX13+1. None the less, our current simulations already show that the thermal–radiative winds can potentially explain all of the wind absorption features seen in GX13+1, so that there is very little room for any additional magnetically driven winds in this source. Acknowledgements We thank K. Hagino for help in setting up the monaco simulations. CD acknowledges STFC funding under grant ST/L00075X/1 and a JSPS long-term fellowship L16581. REFERENCES Agostinelli S. et al.  , 2003, Nucl. Inst. Meth. Phys. Res. A , 506, 250 https://doi.org/10.1016/S0168-9002(03)01368-8 CrossRef Search ADS   Bandyopadhyay R. M., Shahbaz T., Charles P. A., Naylor T., 1999, MNRAS , 306, 417 https://doi.org/10.1046/j.1365-8711.1999.02547.x CrossRef Search ADS   Begelman M. C., McKee C. F., Shields G. A., 1983, ApJ , 271, 70 https://doi.org/10.1086/161178 CrossRef Search ADS   Blandford R. D., Payne D. G., 1982, MNRAS , 199, 883 https://doi.org/10.1093/mnras/199.4.883 CrossRef Search ADS   Clarke C. J., Alexander R. D., 2016, MNRAS , 460, 3044 CrossRef Search ADS   Corbet R. H. D., Pearlman A. B., Buxton M., Levine A. M., 2010, ApJ , 719, 979 https://doi.org/10.1088/0004-637X/719/1/979 CrossRef Search ADS   D'Aì A., Iaria R., Di Salvo T., Riggio A., Burderi L., Robba N. R., 2014, A&A , 564, A62 CrossRef Search ADS   Díaz Trigo M., Boirin L., 2016, Astron. Nachr. , 337, 368 https://doi.org/10.1002/asna.201612315 CrossRef Search ADS   Díaz Trigo M., Sidoli L., Boirin L., Parmar A., 2012, A&A , 543, A50 CrossRef Search ADS   Díaz Trigo M., Migliari S., Guainazzi M., 2014, A&A , 76, 1 Done C., Tomaru R., Takahashi T., 2018, MNRAS , 473, 838 (D18) https://doi.org/10.1093/mnras/stx2400 (D18) CrossRef Search ADS   Fukumura K., Tombesi F., Kazanas D., Shrader C., Behar E., Contopoulos I., 2014, ApJ , 780, 120 https://doi.org/10.1088/0004-637X/780/2/120 CrossRef Search ADS   Fukumura K., Kazanas D., Shrader C., Behar E., Tombesi F., Contopoulos I., 2017, Nat. Astron. , 1, 0062 https://doi.org/10.1038/s41550-017-0062 CrossRef Search ADS   Hagino K., Odaka H., Done C., Gandhi P., Watanabe S., Sako M., Takahashi T., 2015, MNRAS , 446, 663 https://doi.org/10.1093/mnras/stu2095 CrossRef Search ADS   Hashizume K., Ohsuga K., Kawashima T., Tanaka M., 2015, PASJ , 67, 1 https://doi.org/10.1093/pasj/psu132 CrossRef Search ADS   Higginbottom N., Proga D., 2015, ApJ , 807, 107 https://doi.org/10.1088/0004-637X/807/1/107 CrossRef Search ADS   Kallman T., Bautista M., 2001, ApJS , 133, 221 https://doi.org/10.1086/319184 CrossRef Search ADS   Knigge C., Woods J. A., Drew J. E., 1995, MNRAS , 273, 225 https://doi.org/10.1093/mnras/273.2.225 CrossRef Search ADS   Long K. S., Knigge C., 2002, ApJ , 579, 725 https://doi.org/10.1086/342879 CrossRef Search ADS   Luketic S., Proga D., Kallman T. R., Raymond J. C., Miller J. M., 2010, ApJ , 719, 515 https://doi.org/10.1088/0004-637X/719/1/515 CrossRef Search ADS   Miller J. M., Raymond J., Fabian A., Steeghs D., Homan J., Reynolds C. S., van der Klis M., Wijnands R., 2006, Nature , 441, 953 https://doi.org/10.1038/nature04912 CrossRef Search ADS PubMed  Miller J. M. et al.  , 2012, ApJ , 759, L6 https://doi.org/10.1088/2041-8205/759/1/L6 CrossRef Search ADS   Miller J. M., Fabian A. C., Kaastra J., Kallman T., King A. L., Proga D., Raymond J., Reynolds C. S., 2015, ApJ , 814, 87 https://doi.org/10.1088/0004-637X/814/2/87 CrossRef Search ADS   Neilsen J., Rahoui F., Homan J., Buxton M., 2016, ApJ , 822, 20 https://doi.org/10.3847/0004-637X/822/1/20 CrossRef Search ADS   Odaka H., Aharonian F., Watanabe S., Tanaka Y., Khangulyan D., Takahashi T., 2011, ApJ , 740, 103 https://doi.org/10.1088/0004-637X/740/2/103 CrossRef Search ADS   Odaka H., Yoneda H., Takahashi T., Fabian A., 2016, MNRAS , 462, 2366 https://doi.org/10.1093/mnras/stw1764 CrossRef Search ADS   Ponti G., Fender R. P., Begelman M. C., Dunn R. J. H., Neilsen J., Coriat M., 2012, MNRAS , 422, L11 https://doi.org/10.1111/j.1745-3933.2012.01224.x CrossRef Search ADS   Proga D., Kallman T. R., 2002, ApJ , 20, 455 https://doi.org/10.1086/324534 CrossRef Search ADS   Shahbaz T., Charles P. A., King A. R., 1998, MNRAS , 301, 382 https://doi.org/10.1046/j.1365-8711.1998.01991.x CrossRef Search ADS   Shidatsu M., Done C., Ueda Y., 2016, ApJ , 823, 159 https://doi.org/10.3847/0004-637X/823/2/159 CrossRef Search ADS   Sim S. A., Proga D., Miller L., Long K. S., Turner T. J., Reeves J. N., 2010, MNRAS , 408, 1369 https://doi.org/10.1111/j.1365-2966.2010.17215.x CrossRef Search ADS   Ueda Y., Murakami H., Yamaoka K., Dotani T., Ebisawa K., 2004, ApJ , 609, 325 https://doi.org/10.1086/420973 CrossRef Search ADS   Watanabe S. et al.  , 2006, ApJ , 651, 421 https://doi.org/10.1086/507458 CrossRef Search ADS   Waters T. R., Proga D., 2012, MNRAS , 426, 2239 https://doi.org/10.1111/j.1365-2966.2012.21823.x CrossRef Search ADS   Woods D. T., Klein R. I., Castor J. I., McKee C. F., Bell J. B., 1996, ApJ , 461, 767 (W96) https://doi.org/10.1086/177101 (W96) CrossRef Search ADS   APPENDIX A: THE ROTATION VELOCITY FOR RADIAL STREAMLINES Here, we give details of how we calculate the rotation velocity of each element of the wind for radial streamlines (Section 3). We have a linear radial grid, with 20 points from Rin to Rout, so spaced by dR = (Rout − Rin)/20. The inner shell has midpoint R0 = Rin + dR/2. We inject all the mass-loss rate into this radial shell, distributed as (1 − cos θ), on a linear grid of 20 points in θ. Each point on this inner shell is at a horizontal distance of R0sin θ from the black hole. We assume the material has the Keplerian velocity at this horizontal distance i.e. $$v_\phi (R_{0}, \theta ) = \sqrt{GM /(R_{0}\,sin \theta )}$$. Angular momentum conservation along each stream line (of constant θ for these radial streamlines) then gives R0 sin θvϕ(R0, θ) = R sin θvϕ(R, θ) so vϕ(R, θ) = (R/R0)vϕ(R0, θ). APPENDIX B: DENSITY AND VELOCITY STRUCTURE FOR THE DIVERGING STREAMLINES The diverging wind streamlines originate from the focal point which is a distance d below the black hole (see Fig. 4). The innermost edge of the streamlines for the wind is at $$\alpha _{{\rm min}} = \arctan (R_{{\rm in}}/d)$$, and the outer edge is at $$\alpha _{{\rm max}} = \arctan (R_{{\rm out}}/d)$$. We make a linear grid so there are 40 angle elements in the wind, separated by dα = (αmax − αmin)/40 so that αi = α0 + idα for i = 0…40. We have α0 as a free parameter, set by comparison to the results of W96 (see Section 4). The maximum ‘streamline’ length below the disc is from αmax, where $$D=\sqrt{d^2+R_{{\rm out}}^2}$$. We follow this for the same length above the disc. This defines the outer radius of the simulation box which is $$R_{{\rm max}}=2\sqrt{d^2+R_{{\rm out}}^2}$$. We take the inner edge at Rin = 0.1RIC. We superpose a standard θ grid on this (measuring down from the z-axis to radial lines from the centre: Fig. B1). We set θ0 to the point where the innermost streamline edge (at angle α0) reaches Rmax from the origin, and take 41 angles from this to π/2, giving θj(j = 0, 1…40). We make shells using the crossing points of these angles θj with the initial angles αi (Fig. B1). We also define the midpoint angles Ai = (αi + αi + 1)/2 and $$\Theta _j =\frac{1}{2}(\theta _j+\theta _{j+1})$$. Figure B1. View largeDownload slide Details of the model geometry for the diverging wind streamlines. Figure B1. View largeDownload slide Details of the model geometry for the diverging wind streamlines. The velocity along each stream line at a distance lij from its launch point on the disc at radius Ri = d tan Ai is   \begin{equation} v_l (R_i,l_{ij})=f_v c_{{\rm ch}}(R_i)\sqrt{\frac{l_{ij}}{R_i}}, \end{equation} (B1)where   \begin{equation} l_{ij}=D_{ij}-d/\cos A_i, D_{ij}=d\frac{\sin \Theta _j}{\sin (\Theta _j-A_i)} \end{equation} (B2)for a characteristic sound speed $$c_{{\rm ch}}(R_i)=\sqrt{\frac{kT_{{\rm ch}}(R_i)}{\mu m_p}}$$ defined from the characteristic temperature   \begin{equation} T_{{\rm ch}}(R_i)= \left( \frac{L}{L_{{\rm cr}}} \right)^{2/3} (R_i/R_{{\rm IC}})^{-2/3}, \end{equation} (B3)where the critical luminosity, Lcr is   \begin{equation} L_{{\rm cr}}= \frac{1}{8} \left(\frac{m_e}{\mu m_p}\right)^{1/2} \left( \frac{m_e c^2}{kT_{{\rm IC}}} \right)^{1/2} L_{{\rm Edd}} \end{equation} (B4)(see Done et al. 2018) and fv is free parameter which is determined by comparing with the results of W96. We calculate the density of each shell nij assuming mass conservation along each streamline.   \begin{equation} n_{ij}= \frac{\Delta \dot{M}_{wind}(R_i)}{m_I v_l(R_i, l_{ij}) 4 \pi D_{ij}^2 (\cos \alpha _i-\cos \alpha _{i+1})}, \end{equation} (B5) where   \begin{eqnarray} \Delta \dot{M}_{{\rm wind}}(R_i)&=& 2\pi \dot{m}(R_i) R_i \Delta R_i \nonumber\\ &&\times\, 2 =4\pi \dot{m}(R_i) R_i d(\tan \alpha _{i+1}-\tan \alpha _i) \end{eqnarray} (B6)The total mass-loss rate at a given luminosity L is $$\dot{M}_{{\rm wind}} = \sum _i \Delta \dot{M}_{{\rm wind}}(R_i) = 2.0\times 10^{19}\, \mathrm{g\,s^{-1}} (L/L_{{\rm Edd}}=0.3), 8.0\times 10^{18} \,\mathrm{g\,s^{-1}} (L/L_{{\rm Edd}}=0.08), \, 2.1 \times 10^{18} \,\mathrm{g\,s^{-1}} (L/L_{{\rm Edd}}=0.01)$$. Finally we, calculate the column density.   \begin{equation} N_H(\Theta _j)=\sum _i n_{ij}\Delta h_{ij}, \end{equation} (B7)where   \begin{equation} \Delta h_{ij}=d\left(\frac{\sin \alpha _{i+1}}{\sin (\Theta _j-\alpha _{i+1})}-\frac{\sin \alpha _i}{\sin (\Theta _j-\alpha _i)}\right). \end{equation} (B8) We assume Keplerian velocity on the disc plane (θ = π/2 which is at j = 40) so that   \begin{equation} v_{\phi i,40}= \sqrt{\frac{GM}{D_{i,40}\sin A_i}} \end{equation} (B9)and assume the angular momentum conversation along stream line so that   \begin{equation} v_{\phi ij}=\frac{v_{\phi i,40}D_{i,40}\sin A_i}{D_{ij}\sin A_i} =\frac{v_{\phi i,40} D_{i,40}}{D_{ij}}. \end{equation} (B10) © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

Monte Carlo simulations of the detailed iron absorption line profiles from thermal winds in X-ray binaries

Loading next page...
 
/lp/ou_press/monte-carlo-simulations-of-the-detailed-iron-absorption-line-profiles-OdQzGHJKT6
Publisher
Oxford University Press
Copyright
© 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty336
Publisher site
See Article on Publisher Site

Abstract

Abstract Blueshifted absorption lines from highly ionized iron are seen in some high inclination X-ray binary systems, indicating the presence of an equatorial disc wind. This launch mechanism is under debate, but thermal driving should be ubiquitous. X-ray irradiation from the central source heats disc surface, forming a wind from the outer disc where the local escape velocity is lower than the sound speed. The mass-loss rate from each part of the disc is determined by the luminosity and spectral shape of the central source. We use these together with an assumed density and velocity structure of the wind to predict the column density and ionization state, then combine this with a Monte Carlo radiation transfer to predict the detailed shape of the absorption (and emission) line profiles. We test this on the persistent wind seen in the bright neutron star binary GX 13+1, with luminosity L/LEdd ∼ 0.5. We approximately include the effect of radiation pressure because of high luminosity, and compute line features. We compare these to the highest resolution data, the Chandra third-order grating spectra, which we show here for the first time. This is the first physical model for the wind in this system, and it succeeds in reproducing many of the features seen in the data, showing that the wind in GX13+1 is most likely a thermal-radiation driven wind. This approach, combined with better streamline structures derived from full radiation hydrodynamic simulations, will allow future calorimeter data to explore the detail wind structure. accretion, accretion discs, black hole physics, X-rays: binaries, X-rays: individual: (GX 13+1) 1 INTRODUCTION Low-mass X-ray binaries (LMXBs) are systems where a companion star overflows its Roche lobe, so that material spirals down towards a compact object which can be either a black hole or neutron star. The observed X-ray emission is powered by the enormous gravitational potential energy released by this material as it falls inwards, lighting up the regions of intense space–time curvature and giving observational tests of strong gravity. This inflow also powers outflows. Accretion disc winds are seen via blueshifted absorption lines from highly ionized material in high inclination LMXBs (see e.g. the reviews by Ponti et al. 2012; Díaz Trigo & Boirin 2016) but the driving mechanism of these winds is not well understood. The potential candidates are acceleration of gas by the Lorentz force from magnetic fields threading discs (magnetic driving: Blandford & Payne 1982; Fukumura et al. 2014), radiation pressure on the electrons (continuum driving: Proga & Kallman 2002; Hashizume et al. 2015) and thermal expansion of the hot disc atmosphere heated by the central X-ray source, which makes a wind at radii which are large enough for the sound speed to exceed the local escape velocity (thermal driving: Begelman, McKee & Shields 1983; Woods et al. 1996, hereafter W96). Recent work has focused on magnetic driving, primarily because of a single observation of dramatic wind absorption seen from the black hole GRO J1655−40 at low luminosity, far below the Eddington limit where continuum driving becomes important, and with a derived launch radius which is far too small for thermal driving (Miller et al. 2006; Luketic et al. 2010; Higginbottom & Proga 2015). Magnetic wind models can fit this spectrum (Fukumura et al. 2017), but then it is difficult to understand why such high column and (relatively) low ionization winds are not seen from other observations of this object or any other high inclination systems with similar luminosities and spectra. Instead, the unique properties of this wind could potentially be explained if the outflow has become optically thick along the line of sight. This would suppress the observed flux from an intrinsically super-Eddington source (Shidatsu, Done & Ueda 2016; Neilsen et al. 2016). Alternatively, this singular wind may be a transient phenomena, not representative of the somewhat lower column/higher ionization winds which are normally seen. What then powers the more typical winds? Thermal driving qualitatively fits the observed properties, as these winds are preferentially seen in systems with larger discs (Díaz Trigo & Boirin 2016). X-rays from near the compact objects (the inner disc emission and any corona and/or boundary layer) irradiate the outer disc, and the balance between Compton heating and cooling heats the surface to the Compton temperature, which is a luminosity weighted mean energy given by $$T_{\mathrm{IC}} = \frac{1}{4}\int E L(E){\rm d}E/\int L(E) {\rm d}E$$ (Begelman, McKee & Shields 1983; Done, Tomaru & Takahashi 2018, hereafter D18). This temperature is constant with radius, as it depends only on the spectrum of the radiation from the central region, but gravity decreases with distance from the source. For a large enough disc, the isothermal sound speed from this Compton temperature is bigger than the escape velocity, so the heated material makes a transition from forming a bound atmosphere, to an outflowing wind. This defines the Compton radius (RIC = GMmp/kTIC = (6.4 × 104/TIC, 8)Rg, where Rg = GM/c2 and TIC, 8 = TIC/108 K) as the typical launch radius for a thermal wind (Begelman et al. 1983; D18). These thermal wind models give an analytic solution for the mass-loss rate from the disc (Begelman et al. 1983; W96). However, to calculate the observables such as column density and ionization structure requires an understanding of the velocity as a function of two-dimensional position (or equivalently, velocity as a function of length along a streamline, together with the shape of the streamlines). D18 assume a very simple two-dimensional density and velocity structure along radial streamlines, and show that this can match the column seen in the full hydrodynamic simulation of W96. This composite model was applied to multiple spectra of H1743−322, matching the wind seen in its disc dominated state, and predicting the disappearance of these wind absorption lines in a bright low/hard state, as observed (Shidatsu & Done, in preparation). This wind disappearance does not occur from simply the change in photoionization state from the changing illumination (Miller et al. 2012). The key difference is that thermal winds respond to changing illumination by changing their launch radius, density and velocity as well as responding to the changing photoionizing spectrum (Shidatsu & Done, in preparation). Here we explore the thermal wind predictions in more detail, using a three-dimensional Monte Carlo simulation code monaco (Odaka et al. 2011) to calculate the radiation transport through the material so as to predict the resulting emission and absorption line profiles. We calculate these for the very simple constant velocity structure assumed in D18, and then extend this to consider an accelerating wind along biconical streamlines, as a more appropriate disc wind geometry (Waters & Proga 2012). We show results from both these geometries for the specific case of L/LEdd = 0.3, TIC, 8 = 0.13, Rout = 5RIC to compare with W96. We apply the biconical wind model to the bright neutron star binary system GX13+1. This system is unique amongst all the black hole and neutron star binaries in showing persistently strong blueshifted absorption features in its spectrum (Ueda et al. 2004; D'Aì et al. 2014). The X-ray continuum is rather stable, with TIC, 8 ∼ 0.13, making it a good match to the simulations, but it has slightly higher luminosity at L/LEdd = 0.5. This means that radiation pressure should become important, decreasing the radius from which the wind can be launched and increasing its column density. We compare results from a hybrid thermal–radiative wind (calculated using the approximate radiation pressure correction from D18) to the detailed absorption line profiles seen in the third-order Chandra grating data from this source. This is the first quantitative, physical model for the wind, and the first exploration of the highest spectral resolution data from this source. While magnetic driving models are always possible, our results match the majority of the observed features, showing that the wind properties are broadly consistent with hybrid thermal–radiative driving. 2 RADIATIVE TRANSFER CODE We use the Monte Carlo simulation code monaco (Odaka et al. 2011) to calculate radiative transfer through the wind. monaco uses its original physics implementation of photon interactions (Watanabe et al. 2006; Odaka et al. 2011) while it employs the Geant4 toolkit library (Agostinelli et al. 2003) for photon tracking in an arbitrary three-dimensional geometry. We consider an azimuthally symmetric density and velocity field, then use the xstar photoionization code (Kallman & Bautista 2001) to calculate the equilibrium population of ions from each element assuming one-dimensional radiation transfer from the central source. We grid in radial distance, r, and θ and assume the illuminating flux is the transmitted part of the central spectrum along the line of sight to the element. This gives an ionization parameter   \begin{equation} \xi (r, \theta ) = \frac{L \exp (-\tau _{{\rm abs}})}{n (r, \theta ) r^2} \end{equation} (1)where τabs is the optical depth to absorption (but does not include electron scattering), and n(r, θ) is the density as a function of position. We then use monaco to track photons through this ionization structure, including their interaction with this material. Photons interacting with ions can be absorbed in photoionization or photoexcitation, and photons generated via recombination and atomic de-excitation are tracked. Doppler shifts of the absorption cross-sections from the velocity structure of the material are included, as is the Compton energy change on interaction with free electrons. Ideally, this calculated radiation field should then be used as input to xstar, the ion populations recalculated, and the process should be iterated until convergence. However, for our simulations here the wind is mostly optically thin, so we do not include this self-consistent iteration. This method of radiation transfer calculation is based on the Hagino et al. (2015), but the geometry and the velocity/density structure we use is reflected on the thermally driven winds in LMXBs whereas Hagino et al. (2015) is focused on the UV-line driven winds in active galactic nucleus. Also, as thermal winds are typically highly ionized, we consider only H-like and He-like ions of Fe and calculate the spectrum only over a restricted energy band of 6.5–7.2 keV. Table 1 details the transitions used. Table 1. Detailed parameters for each line included in these monaco simulations. Note that we list only lines which have larger oscillator strengths than 10−3. Line ID  Energy (keV)  Oscillator strength  Fe XXV Heαy  6.668  6.57 × 10−2  Fe XXV Heαw  6.700  7.26 × 10−1  Fe XXVI Lyα2  6.952  1.36 × 10−1  Fe XXVI Lyα1  6.973  2.73 × 10−1  Line ID  Energy (keV)  Oscillator strength  Fe XXV Heαy  6.668  6.57 × 10−2  Fe XXV Heαw  6.700  7.26 × 10−1  Fe XXVI Lyα2  6.952  1.36 × 10−1  Fe XXVI Lyα1  6.973  2.73 × 10−1  View Large 3 RADIAL STREAMLINES: D18 3.1 Geometry and parameters We first consider the radial streamline wind model of D18. This calculates the analytic mass-loss rate per unit area, $$\dot{m}(R)$$ where R denotes distance along the disc plane. Integrating over the whole disc gives the total mass-loss rate in the wind, $$\dot{M}$$. This is assumed to flow along radial (centred at the origin) streamlines from a launch radius which is 0.2RIC for high L/LEdd, with constant velocity set at the mass-loss weighted average escape velocity. The mass-loss rate along each radial streamline is weighted with angle such that $$\dot{M}(\theta )\propto \dot{M}(1-\cos \theta )$$, and then mass conservation gives n(r, θ) ∝ (1 − cos θ)/r2. D18 show that these assumptions lead to a total column density through the structure which matches to within a factor 2 of that in the hydrodynamic simulations of W96 (see also Section 4). We put this structure into monaco for L = 0.3LEdd with TIC = 1.3 × 107 K and Rout = 5RIC (mass-loss rate $$\dot{M}_w = 2.0 \times 10^{19}\,\mathrm{g\, s^{-1}}$$ for a 10 M⊙ black hole which means the ratio of mass-loss rate to mass accretion rate $$\dot{M}_w/\dot{M}_a = 3.9$$, where the $$\dot{M}_a = L/(0.1c^2)$$, launch radius of 0.2RIC ≈ 105Rg and weighted average vout = 420 km s−1). We include turbulence, assuming vturb = vout, and calculate the rotation velocity along each stream line assuming angular momentum conservation (see Appendix A). We make a grid which follows the symmetry of the assumed structure, i.e. centred on the origin, with 20 linearly spaced spherical shells from 0.2 to 5RIC, and 20 angles, linearly spaced in θ from 7° to 83° (see below). This density structure is shown in the left-hand panel of Fig. 1, while the right-hand panel shows the mean Fe ion state obtained from the xstar calculation. This is constant along each streamline because the ionization parameter ξ = L/(nr2), and the constant velocity radial streamlines mean that density decreases as 1/r2. Fe is almost completely ionized over the whole grid, with a small fraction of hydrogen-like iron remaining only for high inclination streamlines. Figure 1. View largeDownload slide Distribution of density (left) and mean Fe ionization state (right) for the radial streamline model. Figure 1. View largeDownload slide Distribution of density (left) and mean Fe ionization state (right) for the radial streamline model. 3.2 monaco output Fig. 2 shows the resulting spectra at three different inclination angles. These show that the emission lines are always similarly weak, and that the electron scattered continuum flux makes only a ∼0.5 per cent contribution to the total flux, but that the absorption lines strongly increase at higher inclination angles. We calculate the equivalent width (EW) of each emission and absorption line by fitting the continuum outside the emission and absorption regions with an arbitrary function [F(E) = aEbE + c; Odaka et al. 2016]. The EW of each emission and absorption line is then measured by numerical integration of the difference between the model and the simulation data. The left-hand panel of Fig. 3 shows the EW of the He-like (red) and H-like Lyα2 (green) and Lyα1 (blue) absorption lines. The corresponding emission lines always have EW lower than 0.1 eV so are not seen on this plot. The strong increase of the absorption line EW with inclination clearly shows that the wind is equatorial (by construction from the 1 − cos θ density dependence and constant velocity assumptions). At inclinations above 70°, the Doppler wings of the Kα1 and Kα2 absorption lines merge together due to the turbulent velocities, so Fig. 3 shows only a single EW for this blend. Figure 2. View largeDownload slide Spectra computed for the radial streamline model with 1 eV resolution for different lines of sight. Each panel shows the spectrum in a different inclination angle bin (the angular bin sizes are indicated the top of each panel). The total spectrum is shown in black (top), the spectrum direct photons in red and scattered/reprocessed spectrum is blue (bottom). Note that the vertical axis is plotted linearly in the top panels but logarithmically in the bottom panels. Lines are Fe XXV (6.668 keV for Heαy and 6.700 keV for Heαw) and Fe XXVI (6.952 keV for Lyα2 and 6.973 keV for Lyα1). The equatorial density structure of the wind means that the absorption is much stronger at high inclination angles. The emission is more isotropic, so it can clearly be seen at low inclination angles, but is absorbed by the wind at high inclinations. Figure 2. View largeDownload slide Spectra computed for the radial streamline model with 1 eV resolution for different lines of sight. Each panel shows the spectrum in a different inclination angle bin (the angular bin sizes are indicated the top of each panel). The total spectrum is shown in black (top), the spectrum direct photons in red and scattered/reprocessed spectrum is blue (bottom). Note that the vertical axis is plotted linearly in the top panels but logarithmically in the bottom panels. Lines are Fe XXV (6.668 keV for Heαy and 6.700 keV for Heαw) and Fe XXVI (6.952 keV for Lyα2 and 6.973 keV for Lyα1). The equatorial density structure of the wind means that the absorption is much stronger at high inclination angles. The emission is more isotropic, so it can clearly be seen at low inclination angles, but is absorbed by the wind at high inclinations. Figure 3. View largeDownload slide Left-hand panel: EW of the absorption lines as a function of inclination angle, Fe XXV (Heαw, red) and Fe XXVI (Lyα2, green and Lyα1, blue). The EW of all absorption lines increases strongly at higher inclination, showing the assumed equatorial disc wind geometry. The Doppler wings (with width set by turbulent velocity) of the two H-like absorption lines start to merge for inclinations above 70° so above this we show the total EW of the two lines. Right-hand panel: the blueshifted absorption line velocity for each ion species. This clearly shows the assumed constant velocity structure of the radial streamline. Figure 3. View largeDownload slide Left-hand panel: EW of the absorption lines as a function of inclination angle, Fe XXV (Heαw, red) and Fe XXVI (Lyα2, green and Lyα1, blue). The EW of all absorption lines increases strongly at higher inclination, showing the assumed equatorial disc wind geometry. The Doppler wings (with width set by turbulent velocity) of the two H-like absorption lines start to merge for inclinations above 70° so above this we show the total EW of the two lines. Right-hand panel: the blueshifted absorption line velocity for each ion species. This clearly shows the assumed constant velocity structure of the radial streamline. Fig. 3 (right) shows the outflow velocity, as measured from the energy of the deepest absorption lines (with error set by the resolution of the simulation to ±0.5 eV). These velocities are constant within 25  per cent as a function of inclination, again by construction due to the assumption of constant radial velocity along radial streamlines. 4 DIVERGING WIND In Section 3, we considered a wind model with constant velocity along radial streamlines. However, the expected thermal wind geometry is instead much more like an accelerating, diverging biconical wind Waters & Proga (2012). Full streamline structures which give the density and velocity of the wind at all points can only be found by hydrodynamic calculations (but see Clarke & Alexander 2016 for some analytic approximations). Since modern calculations only exist for the singular case of GRO J1655−40, we follow D18 and use the W96 simulation results. W96 does not give full density/velocity structures, but do give total column density through the wind at three different luminosities. We use these to match to our assumed streamline and velocity structure, which is the standard biconical diverging disc wind used in a variety of systems including cataclysmic variables (Knigge, Woods & Drew 1995; Long & Knigge 2002) and Active Galaxies (Sim et al. 2010; Hagino et al. 2015). 4.1 Geometry and parameters The geometry can be defined by three parameters (Fig. 4). Rin = 0.1RIC, the distance from the source to the inner edge of the wind Rout, the distance from the source to outer edge of the wind αmin, the angle from z-axis to the inner edge of the wind Figure 4. View largeDownload slide The geometry of the diverging biconical wind model. Figure 4. View largeDownload slide The geometry of the diverging biconical wind model. The disc wind is fan-shaped, with a focal point offset down from the centre by a distance d = 0.1RIC/tan αmin so that the wind fills the angles from αmax − αmin down to the disc surface. We use R to denote distance along the disc surface, and r, θ denote radial distance and polar angle from the origin, as before. αmin (or equivalently d) is a free parameter, which sets the wind geometry. Streamlines are assumed to be along lines of constant angle α (where αmin < α < αmax) from the focal point. Distance along a streamline which starts on the disc at radius R is l(R) (see Appendix. A). Velocity along the streamline is assumed to be of the form $$v (r, \theta )=f_v c_{{\rm ch}}(r)\sqrt{\frac{l(r, \theta )}{R(r)}}$$, i.e. this wind accelerates with distance along the streamline, with a terminal velocity which is related via a free parameter fv to the characteristic sound speed cch, given by the balance between heating and cooling in the time it takes the wind to reach a height H ∼ R (D18). The density structure is solved by the mass conservation continuity equation along streamlines (see Appendix B). We calculate the wind properties out to a distance which is twice that of the focal point of the wind to Rout. We set the free parameter values, fv and αmin, and calculate the total column along each line of sight, NH(θ), to the central source for parameters matching to the three W96 simulations. These are L/LEdd = 0.3, 0.08 with Rout = 5RIC and L/LEdd = 0.01 with Rout = 12RIC. We adjust fv and αmin to minimize the difference between our model and W96. We find αmin = 7° and fv = 0.25 matches within a factor 2 of the results from W96. Fig. 5 shows results with these parameters (filled circles), compared to the radial wind model of Section 3 (open circles) as well as the W96 results (solid line). This more physically realistic geometry and velocity gives a similarly good match to the simulations as the D18 radial wind. Figure 5. View largeDownload slide The solid lines show column density as a function of the cosine of the inclination angle through the wind resulting from the hydrodynamic simulations of W96 for L/LEdd = 0.01(green), 0.08 (red), 0.3 (black). The filled circles show that resulting from the diverging biconical wind (Section 4) while the open circles show the radial streamline model of D18 (Section 3). Figure 5. View largeDownload slide The solid lines show column density as a function of the cosine of the inclination angle through the wind resulting from the hydrodynamic simulations of W96 for L/LEdd = 0.01(green), 0.08 (red), 0.3 (black). The filled circles show that resulting from the diverging biconical wind (Section 4) while the open circles show the radial streamline model of D18 (Section 3). The resulting density structure from this different geometry and velocity are shown in the left-hand panel of Fig. 6 for L/LEdd = 0.3. Comparing this with the radial wind shows that the density is higher closer to the disc, and lower further away due to the material accelerating away from the disc rather than being at constant velocity. We run xstar as before, and the right-hand panel of Fig. 6 shows that this leads to a lower mean ionization state of Fe than before, and this is no longer constant along the radial sightline due to the different wind geometry. Figure 6. View largeDownload slide Distribution of density (left) and Fe ionization state (right) for the diverging wind geometry. The accelerating flow gives higher density material close to the disc compared to the constant velocity outflow model in Fig. 1, giving lower mean ionization state. Figure 6. View largeDownload slide Distribution of density (left) and Fe ionization state (right) for the diverging wind geometry. The accelerating flow gives higher density material close to the disc compared to the constant velocity outflow model in Fig. 1, giving lower mean ionization state. 4.2 monaco output We calculate the emission and absorption lines resulting from the different wind structure (Fig. 7). The diverging bipolar wind has higher density material closer to the source compared to the radial wind geometry, so it subtends a larger solid angle to scattering. This means that there is more emission line contribution, as well as a higher fraction of electron scattered continuum (around 2 per cent, see the lower panel of Fig. 7). The left-hand panel of Fig. 8 shows the emission line EW (dotted lines) for each ion species (red: Fe XXV w, green: Fe XXVI Lyα2, blue: Lyα1). These can now be of order 1 eV for face on inclinations, decreasing at higher inclination as they are significantly suppressed by line absorption. Figure 7. View largeDownload slide As in Fig 2, but for the diverging biconical wind geometry. Figure 7. View largeDownload slide As in Fig 2, but for the diverging biconical wind geometry. Figure 8. View largeDownload slide As in Fig. 3, but for the diverging wind model. The lower ionization state means that there is also a contribution from the intercombination line of Fe XXV Heαy (black) at the highest inclinations. Figure 8. View largeDownload slide As in Fig. 3, but for the diverging wind model. The lower ionization state means that there is also a contribution from the intercombination line of Fe XXV Heαy (black) at the highest inclinations. The corresponding absorption line EWs are shown as the solid lines (compare to Fig. 3). The lower mean ionization state leads to more He-like Fe, so there is more of this ion seen in absorption than in the radial streamline model. These absorption lines increase as a function of inclination angle as before, but now the Lyα1 and Lyα2 do not merge together at the highest inclination angles due to the different velocity structure (see right-hand panel of Fig. 8). The lines are formed preferentially in the higher density material close to the disc. The assumed acceleration law means that the typical velocities here are lower than in the constant velocity model, as the material has only just begun to accelerate. Thus, the turbulence is also lower, so the Doppler width of the absorption lines is smaller. This also means that the absorption line saturates to a constant EW at lower column density, so the EW of the absorption lines does not increase so strongly as before at the highest inclination angles. 5 COMPARISON WITH GX13+1 We now use the more physically motivated diverging biconical wind geometry to compare with observational data. An ideal source would be one which is not too different from the parameters simulated in the previous sections, as here we know the total column from W96 and know that our assumed velocity/density matches to this. Of the sources listed in Díaz Trigo & Boirin (2016), the neutron star LMXB GX13+1 is the source which has most similar L/LEdd and TIC to that assumed here, and it also has the advantage that it is a persistent source, with relatively constant luminosity and spectral shape, and it shows similarly strong absorption lines in multiple data sets. 5.1 Observational data GX13+1 was observed by the Chandra HighEnergy Transmission Grating (HETG) four times in 2 weeks in 2010 (Table 2). The first-order data are shown in D'Aì et al. (2014) and reveal multiple absorption features from highly ionized elements (see also Ueda et al. 2004 for similar features in an earlier observation). Higher order grating spectra give higher resolution, as demonstrated for the black hole binaries by (Miller et al. 2015). Here, we show for the first time the third-order Chandra data for GX13+1. We extract first- and third-order HEG spectra from these observations, using ciao version 4.9 and corresponding calibration files. We reprocess the event files with $$`\mathrm{chandra\_repro}$$’, and make response files using ‘mktgres’ to make the redistribution and ancillary response files. We run ‘tgsplit’ to get the HEG ±3 spectra, and run ‘combine_grating_spectra’ to combine HEG plus and minus orders for each observation to derive a single first-order spectrum (black), and a single third-order spectrum (red) as shown in Fig. 9. The first-order spectra can resolve the components of the He-like Fe triplet, with a clear dip to the low energy side at the resonance line energy of 6.7 keV, but the H-like Lyα1 and α2 are blended together. The higher resolution of the third-order spectra is able to clearly separate the He-like intercombination and resonance lines, and even the H-like Lyα1 and α2 (Miller et al. 2015). Figure 9. View largeDownload slide HEG spectra of GX 13+1 from first order (black) and third order (red). Figure 9. View largeDownload slide HEG spectra of GX 13+1 from first order (black) and third order (red). Table 2. List of the Chandra HETG observations. OBSID  MODE  Date  Exposure (ks)  11815  TE-F  24/07/2010  28  11816  TE-F  30/07/2010  28  11814  TE-F  01/08/2010  28  11817  TE-F  03/08/2010  28  OBSID  MODE  Date  Exposure (ks)  11815  TE-F  24/07/2010  28  11816  TE-F  30/07/2010  28  11814  TE-F  01/08/2010  28  11817  TE-F  03/08/2010  28  View Large 5.2 Model of GX 13+1 We fit the contemporaneous RXTE spectrum (ObsID 95338-01-01-05) with a model consisting of a disc, Comptonized boundary layer and its reflection. The resulting inverse Compton temperature of the continuum (disc plus Comptonization) is TIC ∼ 1.2 × 107 K, almost identical to the simulation (see also D'Aì et al. 2014). The luminosity is L = 0.5LEdd (Díaz Trigo, Migliari & Guainazzi 2014; D'Aì et al. 2014), similar to the maximum simulation value of L = 0.3LEdd in W96. The simulation also requires Rout, which can be calculated from the orbital period and mass of binary stars. GX 13+1 has 24 d orbital period, and the neutron star and companion have masses of 1.4 and 5 M⊙, respectively (Bandyopadhyay et al. 1999; Corbet et al. 2010). This gives a binary separation a = 4.6 × 1012 cm, for a Roche lobe radius RR/a = 0.27. The disc size is then Rout = 10RIC assuming that Rout = 0.8RR (Shahbaz, Charles & King 1998), double the value assumed in the simulations. D18 shows that this increase in disc size makes the predicted column slightly larger, but the effect is fairly small (Fig. 3: D18). Fig. 10 (blue line) shows the predicted column density through the wind as a function of inclination angle. This is very similar to the column predicted for the fiducial simulations (Fig. 5) Figure 10. View largeDownload slide The column density as a function of the cosine of the inclination angle for the diverging biconical wind calculated for the system parameters of GX13+1. The blue line shows the predictions for a purely thermal wind, while the red includes a very simple treatment of radiation pressure. The source has L/LEdd ∼ 0.5, so the thermal wind can be launched from closer in due to the lower effective gravity. This effect has a large impact on the predicted column, so the details of how this radiation pressure correction affects the velocity and density structure will be important in determining the line profiles. Figure 10. View largeDownload slide The column density as a function of the cosine of the inclination angle for the diverging biconical wind calculated for the system parameters of GX13+1. The blue line shows the predictions for a purely thermal wind, while the red includes a very simple treatment of radiation pressure. The source has L/LEdd ∼ 0.5, so the thermal wind can be launched from closer in due to the lower effective gravity. This effect has a large impact on the predicted column, so the details of how this radiation pressure correction affects the velocity and density structure will be important in determining the line profiles. However, D18 show that radiation pressure should make a rapidly increasing contribution to the wind as L/LEdd increases from 0.3 to 0.7. The GX13+1 luminosity is midway between these two, so radiation pressure should significantly lower the effective gravity, meaning that the wind can be launched from smaller radii. We follow D18 and estimate a radiation pressure correction to the launch radius of $$\bar{R}_{\mathrm{IC}} = (1.0-0.5L_{\mathrm{Edd}}/0.71L_{\mathrm{Edd}})R_{\mathrm{IC}}=0.30R_{\mathrm{IC}}$$, hence $$R_{\mathrm{out}}=33 \bar{R}_{\mathrm{IC}}$$, dramatically larger than assumed in the fiducial simulations. This correction predicts a density which is 11 times larger and column along any sightline which is 3.3 times larger assuming (as in D18) that the velocity structure is unchanged (red line, Fig. 10). This increase in Rout in terms of RIC means that more wind is produced (as in D18), so the wind efficiency increases to 4.0 (from 2.3). The column density goes close to 1024 cm−2 at high inclinations, so electron scattering becomes important. This effect reduces the illuminating ionizing flux by $$e^{-\tau _T}$$ from the central source along the line of sight to each wind element, and also increases the contribution of diffuse and scattered emission from the wind to the ionizing continuum. We include scattering, reducing the xstar illumination by $$e^{-\tau _{T}}$$ along each line of sight, but do not include the diffuse emission as the time-scale to integrate over the entire wind at each point is prohibitive. We run monaco on this wind structure to predict the detailed absorption line profiles for comparison to the third-order HEG data. Fig. 11 shows the result assuming an inclination angle of 80° (Díaz Trigo et al. 2012) which gives the best fit to the data. This gives a fairly good match to the overall absorption, except for the highest velocity material seen in the data. Lower inclination angles give higher blueshift, but lower absorption line EW, while higher inclination gives larger absorption line but lower blueshift (see Fig. 12). Thus, it is not possible to completely reproduce the observed lines in GX13+1 with our simple radiation pressure corrected thermal wind model. This is not surprising, as radiation pressure will almost certainly change the velocity law by radiative acceleration as well as changing the launch radius. Full radiation hydrodynamic simulations are required to predict the resulting velocity and density structure. None the less, our result demonstrates for the first time that hybrid thermal–radiative wind models can give a good overall match to the column and ionization state of the wind in GX13+1, and that current data can already give constraints on the velocity and density structure of this material. Figure 11. View largeDownload slide The model (red) and HEG third-order spectrum (black). The best-fitting inclination angle is i = 80°. This gives roughly the correct column of Fe XXV and XXVI at low velocity, but fails to match the observed higher velocity blue wing to the absorption features. Figure 11. View largeDownload slide The model (red) and HEG third-order spectrum (black). The best-fitting inclination angle is i = 80°. This gives roughly the correct column of Fe XXV and XXVI at low velocity, but fails to match the observed higher velocity blue wing to the absorption features. Figure 12. View largeDownload slide As in Fig. 8, but with the system parameters of GX 13+1 and the simplest radiation pressure correction to make a hybrid thermal–radiative wind. Figure 12. View largeDownload slide As in Fig. 8, but with the system parameters of GX 13+1 and the simplest radiation pressure correction to make a hybrid thermal–radiative wind. 6 DISCUSSION AND SUMMARY We construct a Monte Carlo code to calculate detailed spectra from any given density and velocity distribution of highly ionized material. We use this to explore the absorption and emission lines of H and He-like Fe for the mass-loss rates predicted from thermal wind models. We first use the radial streamline, constant velocity model of D18 which is able to reproduce the column derived from the hydrodynamic calculations of W96, but then extend this to a more realistic disc-wind geometry with gas accelerating along diverging streamlines, again reproducing the column from W96. The different assumed velocity and density structures for the thermal wind mass-loss rates give different predictions for the overall ionization state of the material, the resulting EW of emission and absorption lines, and their velocity shift. These show the potential of observations to test the detailed structure of the wind. We apply the biconical disc wind model to some of the best data on winds from an LMXB. The neutron star GX13+1 shows strong and persistent absorption features in Chandra first-order HETG spectra (Ueda et al. 2004; D'Aì et al. 2014), but here we show for the first time the higher resolution third-order data. We find that while the source is fairly well matched to the parameters of the brightest fiducial simulation in terms of TIC, the higher luminosity (L/LEdd = 0.5 compared to 0.3 for the simulation) makes a significant impact on the predicted wind properties as it puts the source firmly into the regime where radiation pressure driving should become important. We use the simple radiation pressure correction suggested by D18 and calculate the line profiles from a hybrid thermal–radiative wind. The additional radiation pressure driving means that the wind can be launched from much closer to the central source, and has higher mass-loss rate. This is the first detailed test of the absorption line profiles predicted by physical wind models on any source other than the singular wind seen in GRO J1655−40 (Luketic et al. 2010). Our simulations quantitatively match many of the observed features except for the highest velocity material. This is not surprising, given the simplistic assumptions about the effect of radiation pressure. In future work, we will use the wind velocity and density structure determined from full radiation hydrodynamics simulations in order to properly test the thermal–radiative wind models in GX13+1. None the less, our current simulations already show that the thermal–radiative winds can potentially explain all of the wind absorption features seen in GX13+1, so that there is very little room for any additional magnetically driven winds in this source. Acknowledgements We thank K. Hagino for help in setting up the monaco simulations. CD acknowledges STFC funding under grant ST/L00075X/1 and a JSPS long-term fellowship L16581. REFERENCES Agostinelli S. et al.  , 2003, Nucl. Inst. Meth. Phys. Res. A , 506, 250 https://doi.org/10.1016/S0168-9002(03)01368-8 CrossRef Search ADS   Bandyopadhyay R. M., Shahbaz T., Charles P. A., Naylor T., 1999, MNRAS , 306, 417 https://doi.org/10.1046/j.1365-8711.1999.02547.x CrossRef Search ADS   Begelman M. C., McKee C. F., Shields G. A., 1983, ApJ , 271, 70 https://doi.org/10.1086/161178 CrossRef Search ADS   Blandford R. D., Payne D. G., 1982, MNRAS , 199, 883 https://doi.org/10.1093/mnras/199.4.883 CrossRef Search ADS   Clarke C. J., Alexander R. D., 2016, MNRAS , 460, 3044 CrossRef Search ADS   Corbet R. H. D., Pearlman A. B., Buxton M., Levine A. M., 2010, ApJ , 719, 979 https://doi.org/10.1088/0004-637X/719/1/979 CrossRef Search ADS   D'Aì A., Iaria R., Di Salvo T., Riggio A., Burderi L., Robba N. R., 2014, A&A , 564, A62 CrossRef Search ADS   Díaz Trigo M., Boirin L., 2016, Astron. Nachr. , 337, 368 https://doi.org/10.1002/asna.201612315 CrossRef Search ADS   Díaz Trigo M., Sidoli L., Boirin L., Parmar A., 2012, A&A , 543, A50 CrossRef Search ADS   Díaz Trigo M., Migliari S., Guainazzi M., 2014, A&A , 76, 1 Done C., Tomaru R., Takahashi T., 2018, MNRAS , 473, 838 (D18) https://doi.org/10.1093/mnras/stx2400 (D18) CrossRef Search ADS   Fukumura K., Tombesi F., Kazanas D., Shrader C., Behar E., Contopoulos I., 2014, ApJ , 780, 120 https://doi.org/10.1088/0004-637X/780/2/120 CrossRef Search ADS   Fukumura K., Kazanas D., Shrader C., Behar E., Tombesi F., Contopoulos I., 2017, Nat. Astron. , 1, 0062 https://doi.org/10.1038/s41550-017-0062 CrossRef Search ADS   Hagino K., Odaka H., Done C., Gandhi P., Watanabe S., Sako M., Takahashi T., 2015, MNRAS , 446, 663 https://doi.org/10.1093/mnras/stu2095 CrossRef Search ADS   Hashizume K., Ohsuga K., Kawashima T., Tanaka M., 2015, PASJ , 67, 1 https://doi.org/10.1093/pasj/psu132 CrossRef Search ADS   Higginbottom N., Proga D., 2015, ApJ , 807, 107 https://doi.org/10.1088/0004-637X/807/1/107 CrossRef Search ADS   Kallman T., Bautista M., 2001, ApJS , 133, 221 https://doi.org/10.1086/319184 CrossRef Search ADS   Knigge C., Woods J. A., Drew J. E., 1995, MNRAS , 273, 225 https://doi.org/10.1093/mnras/273.2.225 CrossRef Search ADS   Long K. S., Knigge C., 2002, ApJ , 579, 725 https://doi.org/10.1086/342879 CrossRef Search ADS   Luketic S., Proga D., Kallman T. R., Raymond J. C., Miller J. M., 2010, ApJ , 719, 515 https://doi.org/10.1088/0004-637X/719/1/515 CrossRef Search ADS   Miller J. M., Raymond J., Fabian A., Steeghs D., Homan J., Reynolds C. S., van der Klis M., Wijnands R., 2006, Nature , 441, 953 https://doi.org/10.1038/nature04912 CrossRef Search ADS PubMed  Miller J. M. et al.  , 2012, ApJ , 759, L6 https://doi.org/10.1088/2041-8205/759/1/L6 CrossRef Search ADS   Miller J. M., Fabian A. C., Kaastra J., Kallman T., King A. L., Proga D., Raymond J., Reynolds C. S., 2015, ApJ , 814, 87 https://doi.org/10.1088/0004-637X/814/2/87 CrossRef Search ADS   Neilsen J., Rahoui F., Homan J., Buxton M., 2016, ApJ , 822, 20 https://doi.org/10.3847/0004-637X/822/1/20 CrossRef Search ADS   Odaka H., Aharonian F., Watanabe S., Tanaka Y., Khangulyan D., Takahashi T., 2011, ApJ , 740, 103 https://doi.org/10.1088/0004-637X/740/2/103 CrossRef Search ADS   Odaka H., Yoneda H., Takahashi T., Fabian A., 2016, MNRAS , 462, 2366 https://doi.org/10.1093/mnras/stw1764 CrossRef Search ADS   Ponti G., Fender R. P., Begelman M. C., Dunn R. J. H., Neilsen J., Coriat M., 2012, MNRAS , 422, L11 https://doi.org/10.1111/j.1745-3933.2012.01224.x CrossRef Search ADS   Proga D., Kallman T. R., 2002, ApJ , 20, 455 https://doi.org/10.1086/324534 CrossRef Search ADS   Shahbaz T., Charles P. A., King A. R., 1998, MNRAS , 301, 382 https://doi.org/10.1046/j.1365-8711.1998.01991.x CrossRef Search ADS   Shidatsu M., Done C., Ueda Y., 2016, ApJ , 823, 159 https://doi.org/10.3847/0004-637X/823/2/159 CrossRef Search ADS   Sim S. A., Proga D., Miller L., Long K. S., Turner T. J., Reeves J. N., 2010, MNRAS , 408, 1369 https://doi.org/10.1111/j.1365-2966.2010.17215.x CrossRef Search ADS   Ueda Y., Murakami H., Yamaoka K., Dotani T., Ebisawa K., 2004, ApJ , 609, 325 https://doi.org/10.1086/420973 CrossRef Search ADS   Watanabe S. et al.  , 2006, ApJ , 651, 421 https://doi.org/10.1086/507458 CrossRef Search ADS   Waters T. R., Proga D., 2012, MNRAS , 426, 2239 https://doi.org/10.1111/j.1365-2966.2012.21823.x CrossRef Search ADS   Woods D. T., Klein R. I., Castor J. I., McKee C. F., Bell J. B., 1996, ApJ , 461, 767 (W96) https://doi.org/10.1086/177101 (W96) CrossRef Search ADS   APPENDIX A: THE ROTATION VELOCITY FOR RADIAL STREAMLINES Here, we give details of how we calculate the rotation velocity of each element of the wind for radial streamlines (Section 3). We have a linear radial grid, with 20 points from Rin to Rout, so spaced by dR = (Rout − Rin)/20. The inner shell has midpoint R0 = Rin + dR/2. We inject all the mass-loss rate into this radial shell, distributed as (1 − cos θ), on a linear grid of 20 points in θ. Each point on this inner shell is at a horizontal distance of R0sin θ from the black hole. We assume the material has the Keplerian velocity at this horizontal distance i.e. $$v_\phi (R_{0}, \theta ) = \sqrt{GM /(R_{0}\,sin \theta )}$$. Angular momentum conservation along each stream line (of constant θ for these radial streamlines) then gives R0 sin θvϕ(R0, θ) = R sin θvϕ(R, θ) so vϕ(R, θ) = (R/R0)vϕ(R0, θ). APPENDIX B: DENSITY AND VELOCITY STRUCTURE FOR THE DIVERGING STREAMLINES The diverging wind streamlines originate from the focal point which is a distance d below the black hole (see Fig. 4). The innermost edge of the streamlines for the wind is at $$\alpha _{{\rm min}} = \arctan (R_{{\rm in}}/d)$$, and the outer edge is at $$\alpha _{{\rm max}} = \arctan (R_{{\rm out}}/d)$$. We make a linear grid so there are 40 angle elements in the wind, separated by dα = (αmax − αmin)/40 so that αi = α0 + idα for i = 0…40. We have α0 as a free parameter, set by comparison to the results of W96 (see Section 4). The maximum ‘streamline’ length below the disc is from αmax, where $$D=\sqrt{d^2+R_{{\rm out}}^2}$$. We follow this for the same length above the disc. This defines the outer radius of the simulation box which is $$R_{{\rm max}}=2\sqrt{d^2+R_{{\rm out}}^2}$$. We take the inner edge at Rin = 0.1RIC. We superpose a standard θ grid on this (measuring down from the z-axis to radial lines from the centre: Fig. B1). We set θ0 to the point where the innermost streamline edge (at angle α0) reaches Rmax from the origin, and take 41 angles from this to π/2, giving θj(j = 0, 1…40). We make shells using the crossing points of these angles θj with the initial angles αi (Fig. B1). We also define the midpoint angles Ai = (αi + αi + 1)/2 and $$\Theta _j =\frac{1}{2}(\theta _j+\theta _{j+1})$$. Figure B1. View largeDownload slide Details of the model geometry for the diverging wind streamlines. Figure B1. View largeDownload slide Details of the model geometry for the diverging wind streamlines. The velocity along each stream line at a distance lij from its launch point on the disc at radius Ri = d tan Ai is   \begin{equation} v_l (R_i,l_{ij})=f_v c_{{\rm ch}}(R_i)\sqrt{\frac{l_{ij}}{R_i}}, \end{equation} (B1)where   \begin{equation} l_{ij}=D_{ij}-d/\cos A_i, D_{ij}=d\frac{\sin \Theta _j}{\sin (\Theta _j-A_i)} \end{equation} (B2)for a characteristic sound speed $$c_{{\rm ch}}(R_i)=\sqrt{\frac{kT_{{\rm ch}}(R_i)}{\mu m_p}}$$ defined from the characteristic temperature   \begin{equation} T_{{\rm ch}}(R_i)= \left( \frac{L}{L_{{\rm cr}}} \right)^{2/3} (R_i/R_{{\rm IC}})^{-2/3}, \end{equation} (B3)where the critical luminosity, Lcr is   \begin{equation} L_{{\rm cr}}= \frac{1}{8} \left(\frac{m_e}{\mu m_p}\right)^{1/2} \left( \frac{m_e c^2}{kT_{{\rm IC}}} \right)^{1/2} L_{{\rm Edd}} \end{equation} (B4)(see Done et al. 2018) and fv is free parameter which is determined by comparing with the results of W96. We calculate the density of each shell nij assuming mass conservation along each streamline.   \begin{equation} n_{ij}= \frac{\Delta \dot{M}_{wind}(R_i)}{m_I v_l(R_i, l_{ij}) 4 \pi D_{ij}^2 (\cos \alpha _i-\cos \alpha _{i+1})}, \end{equation} (B5) where   \begin{eqnarray} \Delta \dot{M}_{{\rm wind}}(R_i)&=& 2\pi \dot{m}(R_i) R_i \Delta R_i \nonumber\\ &&\times\, 2 =4\pi \dot{m}(R_i) R_i d(\tan \alpha _{i+1}-\tan \alpha _i) \end{eqnarray} (B6)The total mass-loss rate at a given luminosity L is $$\dot{M}_{{\rm wind}} = \sum _i \Delta \dot{M}_{{\rm wind}}(R_i) = 2.0\times 10^{19}\, \mathrm{g\,s^{-1}} (L/L_{{\rm Edd}}=0.3), 8.0\times 10^{18} \,\mathrm{g\,s^{-1}} (L/L_{{\rm Edd}}=0.08), \, 2.1 \times 10^{18} \,\mathrm{g\,s^{-1}} (L/L_{{\rm Edd}}=0.01)$$. Finally we, calculate the column density.   \begin{equation} N_H(\Theta _j)=\sum _i n_{ij}\Delta h_{ij}, \end{equation} (B7)where   \begin{equation} \Delta h_{ij}=d\left(\frac{\sin \alpha _{i+1}}{\sin (\Theta _j-\alpha _{i+1})}-\frac{\sin \alpha _i}{\sin (\Theta _j-\alpha _i)}\right). \end{equation} (B8) We assume Keplerian velocity on the disc plane (θ = π/2 which is at j = 40) so that   \begin{equation} v_{\phi i,40}= \sqrt{\frac{GM}{D_{i,40}\sin A_i}} \end{equation} (B9)and assume the angular momentum conversation along stream line so that   \begin{equation} v_{\phi ij}=\frac{v_{\phi i,40}D_{i,40}\sin A_i}{D_{ij}\sin A_i} =\frac{v_{\phi i,40} D_{i,40}}{D_{ij}}. \end{equation} (B10) © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off