High-resolution X-ray spectra from stellar objects are an important tool to infer their plasma structures^[1], physical conditions^{[2, 3]}, wind speed in X-ray binary systems^[4] (hereafter W06), coronal abundances (e.g., Blondin et al. 1990^[5]) and other important physical information about astrophysical objects.

Conventionally, the fitting of X-ray spectra is based on reduced method (e.g., Sako et al. 1999^[6], Goldenstein et al. 2004^[7], hereafter S99 and G04) or maximum likelihood (e.g., W06). Among those studies, most continuum models (e.g., power-law and black body) of X-ray spectra are nonlinear with a large number of free parameters (). As to the spectral lines, the effects of satellite lines have not been considered during the fitting of He-like triplets. For the fitting of both continua and lines, it includes many free parameters, which makes the fitting nonlinear and parameters degenerate. In order to decrease the number of free parameters, some parameters were kept fixed during the fitting process, e.g., G04, based on assumptions and previous understanding. In such cases, pitfalls may exist in the estimate of number of degrees of freedom as well as the corresponding error estimates^{[8, 9]}. In Figure 1, we show an example of reduced fit for He-triplets. It can be seen that the model can fit and lines (see Section 3.2.2) but hard to find signals of if a three-Gaussian model was used to fit the three lines simultaneously. Besides, the error bars resulting from the fit were rather large. If the satellite lines are taken into account, i.e., a five/six-Gaussian model is taken to find satellite lines, this method fails to find any line signals.

In recent years, modern Bayesian inference, a method of statistical inference in which Bayes’ theorem is used to update our knowledge of a physical parameter, was introduced in astrophysical studies to fit complex models and to interpret observations (e.g., Reichart et al. 1999^[10], Benítez 2000^[11], Buchner et al. 2014^[12], Walker et al. 2015^[13]). Moreover, Buchner et al.^[12] adopted the Bayesian analysis for X-ray spectra of deep field Active Galactic Nuclei (AGN), and proved its efficiency. For a Bayesian parameter estimate, Markov Chain Monte Carlo (MCMC) is commonly used in numerical applications. In particular, the method is efficient in probing model-based parameters, estimating uncertainties and describing the nonlinear dependencies among the parameters in cases of a large number of model parameters and low signal-to-noise observations^[14]. The advantages of Bayesian analysis over previously used methods are the inclusion of prior information of data, providing possible values and confidence intervals for each parameter. In the present paper we apply an MCMC method to estimate the physical parameters of X-ray spectra, taking the observations of the X-ray binary Vela X-1 as an example.

This paper is organized as follows. Section 2 presents how we determine the parameters of continuum and line emission model through the Bayesian approach. In Section 3, we apply this method to parameterized continuum and He-like triplets’ models to the three different phases for Vela X-1, in order to get the posterior probability density functions (PDFs) of each parameter, to probe the degeneracies among parameters, and hence to get continua and He-like triplets with reasonable fitting error bars. The fitting results are shown and discussed in Section 4, with a summary in Section 5.

In the present paper, we aim to marginalize the posterior PDFs to derive the uncertainties for the parameters of a model fitted to a high-resolution X-ray spectrum. In the Bayesian approach, the posterior PDFs of the model parameters for a given observed spectra, , is^[15]

 (1)

is the likelihood function of given the parameters , whose logarithm form can be written as

 (2)

 (3)

is the prior function. Normally, it represents previous knowledge from other experiments and physical limiting conditions, or any other prior beliefs. Here we use a step function over each parameter

 (4)

To approach our propose, we adopt a stable, well-tested MCMC algorithm, emcee proposed by Foreman-Mackey et al.^[14]. The algorithm is based on the Metropolis–Hasting method (Section 15.8 in Ref. [16]) and an affine-invariant ensemble sampling algorithm called the ‘stretch move’^[17]. The key feature of emcee is that it makes an affine transformation to probability density of parameter during the sampling, which makes it insensitive to covariances among parameters. Therefore, emcee requires hand-turning of only one or two parameters, while the tuning of parameters scales as in an -dimensional space for traditional algorithms. It thus reduces computational costs during the fitting procedure. A complete discussion of the MCMC methods and the algorithms can be found in Ref. [14] and references therein.

In this section, our method is applied to the spectral study of Vela X-1. Vela X-1 is a well-studied X-ray binary, with extensive astronomical simulations of observations^{[4, 6, 7]}, several linked laboratory experiments^[18] and related theoretical modeling^{[19, 20]}, which provides a chance to test the performance of this method through a thorough comparison. Vela X-1 is a pulsing, eclipsing high mass X-ray binary (HMXB) system with highly ionized gases, which was first identified by Ulmer et al.^[21]. Previous observations and theoretical studies have given a reasonably complete description of the system properties. The global structure of the stellar wind in the system was first modeled by S99 using the spectra in the Eclipse phase observed by Solid-State Image Spectrometer. Other physical properties of Vela X-1 were provided by the high-resolution X-ray spectra during different phases observed with the Chandra High-Energy Transmission Grating Spectrometer (HETGS)^{[22, 23]}. G04 first presented the variation of spectra in three orbital phases: Eclipse, , and . In addition, W06 probed the stellar wind dynamics and ionization structures by a quantitative analysis of Doppler shift and line intensities with the same X-ray spectra. The previous literature studies are considered to be sufficiently convincing to be used as benchmarks. With the results of these studies, it is possible for us to verify our analysis method.

The datasets of Vela X-1 used in our analysis were observed with the Chandra HETG/ACIS-S – the High-Energy Transmission Grating/Advanced Imaging Spectrometer, during three different orbital phases centered on Eclipse, and in 2001, with Chandra ObsIDs 1926, 1928 and 1927, respectively. The observation and instrument details can be found in Ref. [24], Ref. [22] and G04. Table 1 summarizes the observation dates and exposure times. The full datasets including counts spectra, ancillary response files (ARFs), and redistribution matrix file (RMFs) are all obtained from public archive Transmission Grating Catalog and Archive (TGCat)¹. Only the first-order events of the high-energy grating (HEG) data are used in our analysis. The data of positive and negative spectral orders are combined to decrease the error bars of observation. Flux correction was done with the Interactive Spectral Interpretation System (ISIS, version 1.6.2-30)² for the observed data to have the same dimensions as the models and to improve the agreement with physical predictions. The background for the observations is not measured independently, which is difficult to be resolved from the source of interest. Practically, the background could be subtracted from the observation signals, as was done in W06. But it is treated as a correction of observation error in the present work because subtracting the background could cause loss of useful signals.

The spectra of Vela X-1 are considered to be composed of continuum and line emission (W06). Therefore, in the present work, for the likelihood function Equation (2), a linear sum of parameterized continuum () and line (), that is, , is used to model the spectra of Vela X-1. For , we mainly focus on He-like triplets which have relatively high -ratio and are commonly used in photoionization plasma diagnostics.

As to the continuum emissions of Vela X-1, we adopt a two-component model (S99), having the form

 (5)

Each set of He-like triplets is described by a multi-Gaussian model, that is,

 (6)

In the continuum fittings of S99 for the Eclipse phase, the photon indices of the two components were set equal because they assumed that photon scattering is elastic, so that it does not alter the spectral shape. In G04, the photon index was fixed to 1.7 in both power laws in the Eclipse phase, while in other phases and was varied. However, the stellar wind is not strictly spherical because of the wake structures including accretion and photoionized wake (see Figure 8 of Ref. [27], a sketch of the different structures in the stellar wind of Vela X-1). The works of Blondin et al.^{[5, 28]} showed that line-of-sight column density, of the stellar wind, is a function of orbital phase (from phase 0.1 to 0.9). That means the indices and may change in different phases. Thus, in our analysis we set and as free parameters in the data-model fitting. In Equation (1), we introduce a parameter which could be considered as an adjustment for the observation error to take the background into account. Therefore, the predictions for are based on seven parameters for each phase, where .

He-like triplets can be detected in the phases of Eclipse and . Therefore, we only fit He-like triplets for the two phases. Among these emission lines, Mg xi and Si xiii with high ratio are adopted as the fitting sample in the present work. In our analysis, the main and satellite lines are fitted simultaneously by multi-Gaussian components for each set of He-like emission lines. Here we describe them as follows.

∙

Mg xi From the observed spectra, two groups of satellite lines can be seen around and lines. Two Gaussians with different widths are used to represent them. In principle, four Gaussians with the same line width should be adopted to model the main lines of Mg xi He-like triplets because the lines correspond to transitions between the shell and the ground state shell. However, the theoretical central energies of and are closed, i.e.,  eV, while the observed full widths at half maximum of the two lines are larger than 1 eV. In this situation, the two lines cannot be resolved. Therefore, three Gaussians with two different line widths are used to represent main lines, that is, the line widths of and are equal. And the line shifts of and lines will not be adopted in the following stellar wind velocity calculation. Finally, the predictions for (Mg xi) are based on 15 parameters .

Selecting suitable upper limits for each parameter in prior functions is based on previous studies (e.g., S99, G04, W06) and observations for Vela X-1. Considering the observed flux and previous studies and experiments, the upper limits of each parameter are listed in Table 2 for each phase. In the estimate of W06, the background event contributions were found to be for the observed flux; thus, the upper limit of is set to for the Eclipse phase and for other phases.

The best-fit parameters and their uncertainties of and are listed in Tables 3–5. Figures 2, 3 and 4 show and implied by our fits using the models of Equations (5) and (6), respectively. The black lines and gray shadows are the observed flux-corrected spectra and observation errors, respectively. The red thick lines show and for the most likely parameters, and the band of light red lines shows a sampling of the posterior PDFs returned by emcee for and . For , the contributions from scatter and direct components are individually shown in blue and purple dashed lines, respectively. Error regions of parameters are directly derived from 68% confidence interval of the posterior PDFs of fitted parameters. It can be seen that the residuals between models and observations scatter around zero, which proves the fits are unbiased.

The results of previous studies are listed in Table 6 and the last two columns of Tables 4 and 5. We note that the errors of the Eclipse phase are larger than those of the other two phases, because the count of this phase is low and of the same order of magnitude of the observation error.

Our best-fit continua are systematically lower () than those of G04 probably due to exclusion of emission lines in our fitting and 28% larger than the ones of W06. For the Eclipse phase, our derived flux is 8% higher than that of S02, but 132% lower than the one of S99. Figure 5 puts the with the parameters from literatures into context.

Considering the confidence ranges of our fitting, our results are in good agreements with previous literature studies (i.e., S99, G04, and W06). For Eclipse, both the results of S99 and G04 fall in our fitting error range in which observed uncertainties are concluded, but our result is systematically lower and the shapes of the continua are different especially in the energy range of 2–4 keV where some He-like triplets such as Si and Mg occur. Our lower results are mainly caused by excluding emission lines during the continuum fitting process. The shape of continuum is determined by photon indices and which are both set as free parameters in the present study. In S99, both of and were fixed to 1.7, and G04 set as free for convenient reason. This is a main cause for such a difference. We adopt Bayesian inference to probe all possible values of photon indices, which help us well in understanding the shape of continua and even the structure of the stellar wind of the HMBX system. For and , in the energy range of 3–10 keV the continua are in very good agreements with those of G04 and W06. But in the range below 3 keV the results are much higher than that in W06. It may be because they used a one-power-law model in their studies. Similarly, excluding emission lines is the main reason why the result of G04 is larger than ours in lower energy range for . The fitted continuum of G04 in is lower than lower limit given by our analysis, and it may be caused by fixing the photon-index value during the fitting process. As we discuss in Section 4.2, when the photon index and other parameters are degenerate, fixing one parameter may cause inaccuracy of parameters. For , all modeled continua cannot fit the observed ones in the range of 8–10 keV. From Figure 2(b), one can see that in this energy range the continuum should be dominated by . However, the models predict a larger contribution from . That is why there is a gap between models and observation in the 8–10 keV range for . Fortunately, as there are no emission lines detected there, the gap would not affect our final results.

for the three phases are all equal to in G04 which implies a spherical stellar wind. In our present fitting, varies from phase to phase, e.g., , and for Eclipse, and , respectively. It implies a non-spherical structure of the stellar wind, which agrees with the simulation in Ref. [28].

As shown in Tables 4 and 5, line flux derived in the present work is systematically lower than that in W06. This is mainly due to two reasons. (1) W06 considered the corrections of interstellar gas absorption in the flux calculation by assuming a hydrogen column density of  cm which corresponds to the density of and the distance of 1.9 kpc. The flux values listed in W06 were compensated by absorption values, which shower higher values. (2) The contributions of satellite lines are subtracted during our flux calculation. Although there is a systematic difference in absolute line flux between the calculations of W06 and ours, after calculating the -ratio^[2] which is used to estimate the plasma temperature with

 (7)

Also shown in Tables 4 and 5, the velocity stellar wind derived from the line shift shows better self-consistency, that is, , and lines of each set of He-like triplets show consistent velocity results in error region, compared with the calculations in W06. The main reason is the inclusion of the satellites during our fitting process. Using Si xiii as an example (see Figures 4(a) and 4(b)), satellite lines gather around 1.845 keV and locate in the blue side of line. In this situation, the derived velocity would be reduced in the Eclipse phase and increased in if satellite lines were not separated from main He-like triplets during the fitting.

Next we estimate the background effects. As we stated before, we did not subtract the background because their source is not well understood. Any inaccuracy in their subtraction would seriously affect the line intensity determination. We treat it as a correction for the observation error.

Its effect, thereby, results in relatively larger uncertainties of parameters compared with previous literature results. From the posterior PDFs (Figure 6, which is made by Python module corner^[30]) of , using the case of  as an example, it is possible for us to estimate the average contributions of the backgrounds. The value is for the and for the , which are larger than the estimates of 3% in W06. The value for the Eclipse phase is . It is in agreement with the maximum of 5% in W06, in which they derive this value from the adjacent region to the dispersed event region of the observed spectrum.

From Figure 6, one can see that all the parameters except the background show degeneracy. It would lead to the exclusion of some reasonable values for the free parameters that are degenerate with the fixed parameters. For instance, in the case of Vela X-1, as Figure 6 shows, the parameters and are degenerate. Setting to be same for all the phases results in the same values of for all the phases, which consequently implies a spherical stellar wind. In the present work, we set all parameters free and find the most likely values in parameter space. The fitted varies from phase to phase due to setting as free. For the photoionized plasma in other astronomical objects, non-fixing any parameters would provide a general approach to deal with their X-ray observed spectra.

In this paper, we introduce the Bayesian approach, which is applied to the archive spectra of Vela X-1 with three different phases: Eclipse, and . We adopt a parameterized two-component power-law model of and a multi-Gaussian model of to predict the continua and He-like triplets, respectively, for all three phases, by setting all parameters as free. Then we fit the observed continua and He-like triplets of Mg xi and Si xiii, by using an MCMC algorithm, emcee to recover the posterior PDFs of all parameters of , and the background. Then we derive best-fit parameters and associated uncertainties for which propagation from the observational errors, uncertainty in the background and the errors from fitting process are all considered. In our results, the column density of scatter component varies from phase to phase, which implies a non-spherical structure of stellar wind. Moreover, our measured wind velocities show very good self-consistency, which provides a reliable approach for the diagnostics of photoionized plasma in the future.