© ESO, 2011

1. Introduction

Outflows from active galactic nuclei (AGN) play an important role in the evolution of the super-massive black holes (SMBH) at the centres of the AGN, as well as on the evolution of the host galaxies and their surroundings. In order to better understand the role of photo-ionised outflows, the geometry must be determined. In particular, our goal is to determine the distance of the photo-ionised gas to the SMBH, which currently has large uncertainties. For this reason we have started an extended monitoring campaign on one of the brightest AGN with an outflow, the Seyfert 1 galaxy Mrk 509 (Kaastra et al. 2011). The main goal of this campaign is to track the response of the photo-ionised gas to the temporal variations of the ionising X-ray and UV continuum. The response time immediately yields the recombination time scale and hence the density of the gas. Combining this with the ionisation parameter of the gas, the distance of the outflow to the central SMBH can be determined.

The first step in this process is to accurately determine the physical state of the outflow: what is the distribution of gas as a function of ionisation parameter, how many velocity components are present, how large is the turbulent line broadening, etc. Our campaign on Mrk 509 is centred around ten observations with XMM-Newton spanning seven weeks in Oct.–Nov. 2009. The properties of the outflow are derived from the high-resolution X-ray spectra taken with the Reflection Grating Spectrometers (RGS, den Herder et al. 2001) of XMM-Newton. The time-averaged RGS spectrum is one of the best spectra ever taken with this instrument, and the statistical quality of this spectrum can be used to improve the current accuracy of the calibration and analysis tools. This creates new challenges for the analysis of the data. The methods developed here also apply to other time-variable sources. Therefore we describe them in some detail in this paper.

For this work we have derived a list of the strongest absorption lines in the X-ray spectrum of Mrk 509, and we perform some simple line diagnostics on a few of the most prominent features. The full time-averaged spectrum will be presented elsewhere (Detmers et al. 2011).

2. Data analysis

Table 1 gives some details on our observations. We quote here only exposure times for RGS. The total net exposure time is 608 ks. No filtering for enhanced background radiation was needed for these observations, because the background was very low and stable for the full duration of our campaign.

The campaign consisted of ten different observations. Each observation was pointed a bit differently to obtain slightly different positions of the spectra on the detectors (the multi-pointing mode of RGS; steps of 0, ± 15′′ and ± 30′′). In this way the holes in the detector caused by bad CCD columns and pixels and CCD gaps, changed position over the spectra, allowing all spectral bins to be sampled. In addition, eventually this procedure will limit outliers in individual wavelength bins from isolated noisy CCD pixels, because the noise of these pixels will be spread over more wavelength bins.

After all exposures were separately processed using the standard RGS pipeline of the XMM-Newton data analysis system SAS version 9, noise from remaining noisy CCD columns and pixels was decreased even more using the following process. Two-dimensional plots of the spectral image (cross-dispersion position versus dispersion angle) and the CCD-pulse height versus spectral dispersion were plotted in the detector reference frame for each separate CCD node, but with all exposures combined. Because the spectral features will be smeared out in these plots, possible CCD defects like hot or dead columns and hot or dead pixels are clearly visible. In this way a number of additional bad columns and pixels were manually identified for each detector CCD node and were added to the SAS bad pixel calibration file. After this, all exposures were reprocessed with the new bad pixel tables.

3. Combination of spectra

To obtain an accurate relative effective area and absolute wavelength calibration, we need the signal-to-noise of the combined 600 ks spectrum. Therefore we first describe how the individual spectra are accurately combined, and discuss the effective area and wavelength scale later.

3.1. Combination of spectra from the same RGS and spectral order

Detectors like the RGS yield spectra with lacking counts in parts of the spectrum owing to the presence of hot pixels, bad columns, or gaps between CCDs. This leads to problems when spectra taken at different times have to be combined into a single spectrum, in particular when the source varies in intensity and the lacking counts fall on different wavelengths for each individual spectrum. For RGS this is for example the case when the multi-pointing mode is used like in our present Mrk 509 campaign, or when data of different epochs with slightly different pointings are combined, or if bad columns have (dis)appeared between different observations because of the transient nature of some bad pixels.

Under the above conditions the SAS task rgscombine should not be used. We illustrate this with a simplified example (Fig. 1). Consider two spectra of the same X-ray source, labelled 1 and 2, with seven spectral bins each (Fig. 1). The source has a flat spectrum and the flux per bin in spectrum 1 is twice the flux of spectrum 2. In spectrum 1, bin A is missing (zero counts), and in spectrum 2, bin B is missing. The combined spectrum has 600 counts, except for bins A and B with 200 and 400 counts, respectively. However, when the response matrices of both observations are combined, it is assumed that bins A and B were on for 50% of the time, hence the model predicts 300 counts for those bins. Clearly, the data show a deficit or excess at those bins, which an observer could easily misinterpret as an additional astrophysical absorption or emission feature at bins A and B.

How can this problem be solved? One possibility is to fit all spectra simultaneously. However, when there are many different spectra (like the 40 spectra obtained from two RGS detectors in two spectral orders for ten observations that we have for our Mrk 509 data), this can be a cumbersome procedure because the memory and CPU requirements become very demanding. Another option is to rigorously discard all wavelength bins where even during short periods (a single observation or a part of an observation) data are lacking. This solution will cause a significant loss of diagnostic capability, because many more bins are discarded and likely several of those bins will be close to astrophysically interesting features.

Here we follow a different route. It is common wisdom (and this is also suggested by the SAS manual) not to use fluxed spectra for fitting analyses. The problem is the lack of a well-defined redistribution function for that case. We will show below in Sect. 4 that we can obtain an appropriate response matrix for our case, which allows us to make an analysis with fluxed spectra under the conditions given in that section.

We first create individual fluxed spectra using the SAS task rgsfluxer. We use exactly the same wavelength grid for each observation. We then average the fluxed spectra using for each spectral bin the exposure times as weights. This is repeated for all bins having 100% exposure and for all spectra that are to combined.

For bins with lacking data (like bins A and B of Fig. 1) we have to follow a different approach. If the total exposure time of spectrum k is given by t_k, and the effective exposure time of spectral bin j in spectrum k is given by t_kj, we clearly have t_kj < t_k for “problematic” bins j, and we will include in the final spectrum only data bins with t_kj > ft_k, with 0 ≤ f ≤ 1 a tunable parameter discussed below. However, because of the variability of the source, the average flux level of the included spectra k for problematic bin j may differ from the flux level of the neighbouring bins j − 1 and j + 1 which are fully exposed. To correct for this, we assume that the spectral shape (but not the overall normalisation) in the local neighbourhood of bin j is constant in time. From these neighbouring bins (with fluxes F_{j − 1} and F_{j + 1}) we estimate the relative contribution R to the total flux for the spectra k that have sufficient exposure time t_kj for bin j, i.e. (1)and similar for the other neighbour at j + 1. We interpolate the values for R linearly at both sides to obtain a single value at the position of bin j. The flux for bin j is now estimated as (2)Example: for bin A in Fig. 1 we have R = 200 ×   1/(400 ×   1 + 200 × 1) = 1/3, and hence using the flux measurement in spectrum 2, we have F_A = (200 ×   1)/((1/3) ×   2) = 300 counts s^-1, as it should be.

This procedure gives reliable results as long as the spectral shape does not change locally; because there is less exposure at bin j, the error bar on the flux will obviously be larger. However, when there is reason to suspect that the bin is at the location of a spectral line that changes in equivalent width, this procedure cannot be applied!

To give the user more flexibility, we have introduced the minimum exposure fraction f in (1). For f = 0, we obtain the best results for spectra with constant shape (because every bit of information that is available is used). On the other hand, if f = 1, only those data bins will be included that have no bad bins in any of the observations to be combined. The advantage in that case is that there is no bias in the stacked spectrum, but a disadvantage is that a significant part of the spectrum may be lost, for example near important diagnostic lines. This can be a problem in particular for the multi-pointing mode, the purpose of which is to have different wavelengths fall on different parts of the spectrum and thus to have a measurement of the spectrum at all wavelengths. We illustrate this for the region near the chip gap between CCD 5 and 6 of RGS2 (see Fig. 2). For f < 1 some data points have large error bars owing to the low effective exposure of some bins that contain missing columns, and the corresponding low number of counts N yields large relative errors (~N^-0.5). In our analysis we use f = 0 throughout.

Finally, in the combined spectrum it is still possible because of the binning and randomisation procedures within the SAS task rgsfluxer, that despite careful screening for bad pixels, a few bad bins remain in the fluxed spectrum, often adjacent or close to discarded bins. For this reason, we check for any bins that are more than 3σ below their right or left neighbouring bin (taking the statistical errors of both into account). Typically, the algorithm finds a few additional bad bins in an individual observation, which we also discard from our analysis. Only for very strong isolated emission lines with more than 1500 counts in a single observation our method would produce false rejections near the bending points of the instrumental spectral redistribution function, because then the spectral changes for neighbouring bins are stronger than 3σ statistical fluctuations, but our spectra of Mrk 509 do not contain such sharp and strong features.

We have implemented the procedures presented in this section in the programme rgs_fluxcombine that is available within the publicly available SPEX distribution¹.

Finally, we note that there is in principle another route to combine the spectra. This would be to run rgscombine but to modify by a user routine the column AREASCAL so that the count rates are properly corrected for the missing bins. It can be shown that in the combined spectrum AREASCAL must be scaled by 1/R, with R derived as described above, in order to obtain robust results. It requires the development of similar tools as rgs_fluxcombine, however, to implement this modification.

3.2. Combination of spectra from different RGS or spectral order

In the previous section we showed how to combine spectra from the same RGS and spectral order, but for different observations. However, we also intend to combine RGS1 and RGS2 spectra, and spectra from the first and second spectral order. We do this by simply averaging the fluxed spectra, using for each bin the statistical errors on the flux as weight factors. Because of the lower effective area in the second spectral order, second-order spectra are allotted less weight than first order spectra in this way, their weights are typically between 20% (near 18 Å) and 50% (near 8 Å) of the first-order spectral weights.

In Sect. 4 we show how to create a response matrix for the fluxed RGS spectra. The same weights that are used to determine the relative contributions of the different spectra to the combined spectrum are also used to weigh the contributions from the corresponding redistribution functions. To avoid discontinuities near the end points of second-order spectra or near missing CCDs, we pay attention to small effective area corrections, as outlined in Sect. 5.

4. Response matrix

We have explained above that for time-variable sources standard SAS tasks like rgscombine cannot be used to combine the ten individual spectra of Mrk 509 or any other source into a single response matrix. For that reason, we have combined fluxed spectra as described in the sections above. The only other feasible alternative would be to fit the ten individual spectra simultaneously. But with two RGS detectors, two spectral orders and ten observations, this adds up to 40 individual spectra. The memory occupied by the corresponding response matrices is 2.0 Gbyte. Although our fitting program SPEX is able to cope with this, fitting becomes cumbersome and error searches extremely slow, owing to the huge number of matrix multiplications that have to be performed.

For our combined fluxed spectra, fitting is much faster because we have only a single spectrum and hence only need one response matrix. Because we use fluxed spectra, we can use a simple response matrix with unity effective area, and the spectral redistribution function given by the redistribution part of the RGS matrix. Unfortunately, the RGS response matrix produced by the SAS rgsrmfgen task combines the effective area and redistribution part into a single data file. Furthermore, because of the multi-pointing mode that we used and because of transient bad columns, the matrix for each of the ten observations will be slightly different.

Therefore, we have adopted the following approach. For a number of wavelengths (7 to 37 Å, step size 2 Å) we have fitted the RGS response to the sum of three or four Gaussians. This gives a quasi-diagonal matrix resulting in another speed-up of spectral analysis. For these fits to the redistribution function we have omitted the data channels with incomplete exposure (near chip boundaries, and at bad pixels). The parameters of the Gaussians (normalisation, centre, and width) were then modelled with smooth analytical functions of wavelength. We paid most attention to the peak of the redistribution function, where most counts are found.

The following parameterisation describes the RGS redistribution functions well (all units are in Å):

The relevant parameters are given in Table 2, omitting parameters that are not used (zero). We cut off the redistribution function beyond ± 1 Å from the line centre. We have compared these approximations to the true redistribution function (see Figs. 3, 4), and found that our model accurately describes the core down to a level of about 1% of the peak of the redistribution function. Also, the flux outside the ± 1 Å band is less than ~1% of the total flux for all wavelengths. It should be noted that the above approximation is more accurate than the calibration of the redistribution itself. Currently, the width of the redistribution is known only to about 10% accuracy².

We finally remark that Fig. 3 seems to suggest that the response also contains the second-order spectra, for instance for 9 Å photons, apart from the main peak near 9 Å, there are secondary peaks at 18, 27, and 36 Å. These peaks are, however, not the genuine higher order spectra, but they correspond to higher order photons that by chance end up in the low-energy tail of the CCD redistribution function. For instance, a photon with intrinsic energy 1.38 keV would be seen in first order at 9 Å, but in second order at 18 Å. A small fraction of these second order events end up in the CCD redistribution tail and are not detected near 1.38 keV but near 0.69 keV. Because their (second-order) measured wavelength is 18 Å, they are identified as first-order events (energy 0.69 keV, wavelength 18 Å). As can be seen from Fig. 3, this is only a small contribution, and we can safely ignore them.

5. Effective area corrections

The statistics of our stacked 600 ks spectrum of Mrk 509 is excellent, with statistical uncertainties down to 2% per 0.1 Å per RGS in first order. In larger bins of 1 Å width, the quality is even better by a factor of ~3. These statistical uncertainties are smaller than the accuracy to which the RGS effective area has been calibrated. A close inspection of the first- and second-order spectra of our source for both RGS detectors showed that there are some systematic flux differences of up to a few percent at wavelength scales larger than a few tenths of an Å (Figs. 5, 6), consistent with the known accuracy of the effective area calibration. These differences sometimes correlate with chip boundaries (like at 34.5 Å), but not always. A preliminary analysis shows the same trends for stacked RGS spectra of the blazar Mrk 421.

While a full analysis is underway, we use here a simple approach to correct these differences and obtain an accurate relative effective area calibration. We model the differences purely empirically with the sum of a constant function and a few Gaussians with different centroids, widths, and amplitudes, as indicated in Figs. 5, 6. We then attribute half of the deviations to RGS1 and the other half to RGS2. In this way there is better agreement between the different RGS spectra, in particular across spectral regions where because of the failure of one of the chips the flux would otherwise show a discontinuity at the boundary between a region with two chips and a region with only one chip. Because the effective area of the second-order spectra has been calibrated with less accuracy than the first-order spectra, we adjust the second-order spectrum to match the first-order spectrum.

We obtain the following corrections factors: Here r_ij is the flux of RGS i order j relative to RGS2 order 1 as shown in Figs. 5, 6, and G(λ,N,μ,σ) ≡ Ne ^{− (λ − μ)2/2σ2}. With these definitions, the fluxes need to be multiplied by correction factors a_ij to obtain the corrected values as follows: The above approach may fail whenever the small discrepancy between one spectrum and the other is caused by a single RGS, and not by both. Then there will be a systematic error in the derived flux with the same Gaussian shape as used for the correction factors, but with half its amplitude. This is compensated largely, however, by our approach where we fit the true underlying continuum of Mrk 509 by a spline, and in our analysis we make sure that we do not attribute real broad lines to potential effective area corrections. The adopted effective area corrections (the inverse of the a_ij factors) are shown in Fig. 7. Note that for wavelength ranges with lacking CCDs the effective area corrections are based on the formal continuation of the above equations, and because of that may differ from unity. However, this is still consistent with the absolute accuracy of the RGS effective area calibration.

6. Wavelength scale

6.1. Aligning both RGS detectors

The zero-point of the RGS wavelength scale has a systematic uncertainty of about 8 mÅ (den Herder et al. 2001). According to the latest insights, this is essentially caused by a slight tilting of the grating box because of Solar irradiance on one side of it, which causes small temperature differences and hence inhomogeneous thermal expansion. The thermal expansion changes the incidence angle of the radiation on the gratings and hence the apparent wavelength. The details of this effect are currently under investigation. Owing to the varying angle of the satellite with respect to the Sun during an orbit and the finite thermal conductivity time-scale of the relevant parts, the precise correction will also depend on the history of the satellite orbit and pointing.

Because corrections for this effect are still under investigation, we take here an empirical approach. We start with aligning both RGS detectors. We have taken the nine strongest spectral lines in our Mrk 509 spectrum that can be detected on both RGS detectors, and have determined the difference between the average measured wavelength for both RGS detectors in the combined spectrum. We found that the wavelengths measured with RGS2 are 6.9 ± 0.7 mÅ longer than the wavelengths measured with RGS1 (Fig. 8).

Therefore, we have re-created the RGS2 spectra by sampling them on a wavelength grid similar to the RGS1 wavelength grid, but shifted by +6.9 mÅ. Using this procedure, spectral lines always end up in the same spectral bins. This allows us to stack RGS1 and RGS2 spectra by simply adding up the counts for each spectral bin, without a further need to re-bin the counts when combining both spectra.

After applying this correction, we stacked the spectra of both RGS detectors and compared the first- and second-order spectra. We found that the second-order centroids agree well with the first-order centroids (to within an accuracy of 3 mÅ) and therefore did not apply a correction.

6.2. Strongest absorption lines

To facilitate the derivation of the proper wavelength scale, we have compiled a list of the strongest spectral lines that are present in the spectrum of Mrk 509 (Table 3). We have measured the centroids λ₀ and equivalent widths EW of these lines by fitting the fluxed spectra F(λ) locally in a 1 Å wide band around the line centre using Gaussian absorption lines, superimposed on a smooth continuum for which we adopted a quadratic function of wavelength. We took care to re-formulate this model in such a way that the equivalent width is a free parameter of the model and not a derived quantity based on the measured continuum or line peak: (17)Because in this paper we focus entirely on absorption lines, we report here EWs for absorption lines as positive numbers. Where needed we added additional Gaussians to model the other strong but narrow absorption or emission lines that are present in the same 1 Å wide intervals around the line of interest. The Gaussians include the instrumental broadening. For the strongest 14 lines we could leave the width of the Gaussians a free parameter; for the others we constrained the width to a smooth interpolation of the width from these stronger lines.

For illustration purposes, we show the full fluxed spectrum (after all corrections that we derived in this paper have been applied) in Fig. 9. A more detailed view of the spectrum including realistic spectral fits is given in Paper III of this series (Detmers et al. 2011).

6.3. Doppler shifts and velocity scale

The spectra as delivered by the SAS are not corrected for any Doppler shifts. We have estimated that the time-averaged velocity of the Earth with respect to Mrk 509 is −29.3 km s^-1 for the present Mrk 509 data. Because the angle between the lines of sight towards the Sun and Mrk 509 is close to 90° during all our observations, the differences in Doppler shift between the individual observations are always less than 1.0 km s^-1, with an rms variation of 0.5 km s^-1. Therefore these differences can be ignored (less than 0.1 mÅ for all lines). The orbital velocity of XMM-Newton with respect to Earth is less than 1.2 km s^-1 for all our observations and can also be neglected.

We need a reference frame to determine the outflow velocities of the AGN. For the reference cosmological redshift we choose here 0.034397 ± 0.000040 (Huchra et al. 1993), because this value is also recommended by the NASA/IPAC Extragalactic database NED, although the precise origin of this number is unclear; Huchra et al. (1993) refer to a private communication to Huchra et al. in 1988.

Because we do not correct the RGS spectra for the orbital motion of the Earth around the Sun, we will instead use an effective redshift of 0.03450 ± 0.00004 to calculate predicted wavelengths. This effective redshift is simply the combination of the 0.034397 cosmological value with the average −29.3 km s^-1 of the motion of the Earth away from Mrk 509.

For Galactic foreground lines it is common practice to use velocities in the Local Standard of Rest (LSR) frame, and we will adhere to that convention for Galactic lines. The conversion from LSR to Heliocentric velocities for Mrk 509 is given by v_LSR = v_Helio + 8.9   km   s^-1, cf. Sembach et al. (1999).

6.4. Absolute wavelength scale

The correction that we derived in Sect. 6.1 aligns both RGS detectors and spectral orders, but there is still an uncertainty (wavelength offset) in the absolute wavelength scale of RGS1. We estimate this offset here by comparing measured line centroids in our spectrum to predicted wavelengths based on other X-ray, UV, or optical spectroscopy of Mrk 509. We consider five different methods:

• 1.

  lines from the outflow of Mrk 509;

• 2.

  lines from foreground neutral gas absorption;

• 3.

  lines from foreground hot gas absorption;

• 4.

  comparison with a sample of RGS spectra of stars;

• 5.

  comparison with Chandra LETGS spectra.

6.4.1. Lines from the outflow of Mrk 509

The outflow of Mrk 509 has been studied in the past using high-resolution UV spectroscopy with FUSE (Kriss et al. 2000) and HST/STIS (Kraemer et al. 2003). A preliminary study of the COS spectra taken simultaneously with the Chandra LETGS spectra in December 2009, within a month from the present XMM-Newton observations, shows that there are no major changes in the structure of the outflow as seen through the UV lines. Hence, we take the archival FUSE and STIS data as templates for the X-ray lines. Kriss et al. (2000) showed the presence of seven velocity components labelled 1–7 from high to low outflow velocity. Component 4 may have sub-structure (Kraemer et al. 2003) but we ignore these small differences here. The outflow components form two distinct troughs, a high-velocity component from the overlapping components 1–3 and a low-velocity component from the overlapping components 4–7.

Note that our analysis of ISM lines in the HST/COS spectra for our campaign (Kriss et al. 2011) shows that the wavelength scale for the FUSE data (Kriss et al. 2000) needs a slight adjustment. As a consequence, we subtract 26 km s^-1 from the velocities reported by Kriss et al. (2000), in addition to the correction to a different host galaxy redshift scale.

In Table 4 we list the column-density weighted or equivalent width weighted average velocities. We also give the typical ionisation parameter ξ for which the ion has its peak concentration, based upon the analysis of archival XMM-Newton data of Mrk 509 by Detmers et al. (2010).

In principle, equivalent-width weighted line centroids are preferred to predict the centroids of the unresolved X-ray lines. Because for most of the X-ray lines the optical depths are not very high, using just the column-density weighted average will not make much difference, however.

The strong variation of average outflow velocity as a function of ξ that is apparent in Table 4, in particular between O vi and N v, shows that caution must be taken. Therefore, we consider here only X-ray lines from ions that overlap in ionisation parameter with the ions in Table 4. We list those lines in Table 5. For C v and O v we have interpolated the velocity between the O vi and N v velocity according to their log ξ value. The weighted average shift Δλ ≡ λ_obs − λ_pred is −5.3 ± 3.6 mÅ.

6.4.2. Lines from foreground neutral gas absorption

Our spectrum (see Table 3) contains lines from the neutral as well as the hot phase of the ISM. Here we discuss the neutral phase. The spectrum shows two strong lines from this phase, the O i and N i is–2p lines. Other lines are too weak to be used for wavelength calibration purposes.

The sight line to Mrk 509 is dominated by neutral material at LSR velocities of about +10 km s^-1 and +60 km s^-1, as seen for instance in the Na i D₁ and D₂ lines, and the Ca ii H and K lines (York et al. 1982). Detailed 21 cm observations (McGee & Newton 1986) show four components: two almost equal components at −6 and +7 km s^-1, a weaker component (6% of the total H i column density) at +59 km s^-1 and the weakest component at +93 km s^-1 containing 0.7% of the total column density.

The strongest X-ray absorption line from the neutral phase is the O i 1s–2p transition. The rest-frame wavelength of this line has been determined as 23.5113 ± 0.0018 Å (Kaastra et al. 2010). Based on the velocity components of McGee & Newton (1986), we have estimated the average velocity for this line as well as the UV line at 1302 Å, assuming an oxygen abundance of 5.75 × 10^-4 (proto-solar value: Lodders 2003), adopting that ~50% of all neutral oxygen is in its atomic form. We predict for the λ23.5 Å line an LSR velocity of +24 km s^-1, and for the λ1302 Å line a velocity of +40 km s^-1. For a two times higher column density, these lines shift by no more than +6 and +1 km s^-1, respectively. For comparison, the measured centroid of the λ1302 Å line (Collins et al. 2004, Fig. 6) is +40 km s^-1, in close agreement with the prediction. A systematic uncertainty of ± 10 km s^-1 on the predicted λ23.5 Å line seems therefore to be justified. Taking this into account, we expect the O i line at 23.515 ± 0.002 Å.

Unfortunately, the Galactic O i λ23.5 Å line is blended by absorption from lines of two other ions in the outflow of Mrk 509, O iv and Ca xv. We discuss those lines below.

Contamination by O iv:

The O iv 1s–2p ²P transition, which has a laboratory wavelength of 22.741 ± 0.004 Å (Gu et al. 2005), contaminates the Galactic O i line. In fact, the situation is more complex, because the main O iv 1s–2p multiplet has four transitions. In Table 6 we summarise the theoretical calculations and experimental measurements of the wavelengths of these lines. The adopted values are based on the single experimental measurement and the HULLAC calculations, while the adopted uncertainties are based on the scatter in the wavelength differences.

In the present case, both the O iv 1s–2p ²D_3/2 and 1s–2p ²P_1/2 and ²P_3/2 transitions of the outflow blend with the Galactic O i line, with relative contributions of 37, 42, and 21%, respectively. The precise contamination is hard to determine, because these are the strongest O iv transitions and the other O iv transitions are too weak to determine the column density. Based on the absorption measure distribution of other ions with similar ionisation parameter, we estimate that the total equivalent width of these three O iv lines is 3.6 ± 2.0 mÅ (Detmers et al. 2011). The same modelling predicts that the O iv concentration peaks at log ξ = −0.6, at almost the same ionisation parameter as for C iv. The optical depths of all X-ray lines of O iv is less than 0.1, hence we use the column-density weighted average line centroid of C iv (−197 km s^-1) for O iv. With this, we predict a wavelength for the combined ²D and ²P transitions of O iv of 23.523 ± 0.004 Å.

Contamination by Ca xv:

The other contaminant to the Galactic O i 1s–2p line is the Ca xv 2p–3d (³P₀–³D₁) λ22.730 Å transition of the outflow. This is the strongest X-ray absorption line for this ion. Based on the measured column density of Ca xiv in the outflow, and assuming that the Ca xv column density is similar, we expect an equivalent width of 4.9 ± 2.3 mÅ for this line, about 10% of the strength of the O i blend.

The rest-wavelength of this transition is somewhat uncertain, however. Kelly (1987) gives a value of 22.725 Å, without errors, and refers to Bromage & Fawcett (1977), who state that their theoretical wavelength is within 5 mÅ from the observed value of 22.730 Å, and that the line is blended. Although not explicitly mentioned, their experimental value comes from Kelly & Palumbo (1973) (cited in the Bromage & Fawcett paper). Kelly & Palumbo (1973) give a value of 22.73 Å, but state it is for the ³P₂–³D₃ transition, which is a misclassification according to Fawcett & Hayes (1975). These latter authors list as their source Fawcett (1971), who gives the value of 22.73, but no accuracy. In Fawcett (1970) the typical accuracy of the plasma measurements is limited by Doppler motions of about 10^-4c. If we assume that this accuracy also holds for the data in Fawcett (1971), the estimated accuracy should be about 2 mÅ. Given the rounding of Fawcett (1970), an rms uncertainty of 3 mÅ is more appropriate. However, there is some level of blending with other multiplets. Blending is probably caused by other 2p–3d transitions in this ion, most likely transitions between the ³P–³D, ³P–³P and ¹D–¹F multiplets (Aggarwal & Keenan 2003). A comparison of the theoretical wavelengths of Aggarwal & Keenan (2003) with the experimental values of Kelly (1987) for all lines of the ³P–³D multiplet shows a scatter of 6 mÅ, apart from an offset of 145 mÅ.

Given all this, we adopt here the original measurement of Fawcett (1971), but with a systematic uncertainty of 5 mÅ. As Ca xv has a relatively large ionisation parameter (log ξ = 2.5), we expect that like other highly ionised species in Mrk 509 it is dominated by the high-velocity outflow components, which have a velocity of about −340 km s^-1 (see Table 4). With a conservative 200 km s^-1 uncertainty (allowing a centroid up to halfway between the two main absorption troughs), we may expect the line at 23.488 ± 0.017 Å.

Centroid of the O i blend:

The total equivalent width of the O i blend is 29.1 ± 2.8 mÅ (see Table 3). The velocity dispersion of the contaminating O iv and Ca xv lines is an order of magnitude higher than the low velocity dispersion of the cold Galactic gas, and both contaminants have optical depth lower than unity. The total blend is thus the combination of broad and shallow contaminants with deep and narrow O i components. Thus the contribution from O i itself is obtained by simply subtracting the contributions from O iv and Ca xv from the total measured equivalent width of the blend. This yields an O i equivalent width of 20.6 ± 4.2 mÅ. Taking into account the contribution of the contaminants, we expect the centroid of the blend to be at 23.5115 ± 0.0032 Å. This corresponds to an offset Δλ = 9.9 ± 4.5 mÅ.

N i:

The other strong X-ray absorption line from the neutral phase of the ISM is the N i 1s–2p line. We adopt the same LSR velocity as for the unblended O i line, +24 km s^-1. This is consistent with the measured line profile of the N i λ1199.55 Å line, which shows two troughs at about 0 and +50 km s^-1 (Collins et al. 2004). The shift derived from the N i line is therefore +12.2 ± 8.9 mÅ.

We note, however, that for N i we can only use RGS2, because RGS1 shows some instrumental residuals near this line. The nominal wavelength for RGS1 is 26 mÅ higher than for RGS2. It is therefore worthwhile to do a sanity check. For this we have determined the same numbers from archival data of Sco X-1. That source was observed several times with significant offsets in the dispersion direction, so that the spectral lines fall on different parts of the CCDs or even on different CCD chips. We measure wavelengths for O i and N i in Sco X-1 of 23.5123 ± 0.0015 and 31.2897 ± 0.0030 Å, respectively. Sco X-1 has a heliocentric velocity of −114 km s^-1 (Steeghs & Casares 2002). Given its Galactic longitude of 07, we expect the local foreground gas to have negligible velocity. This is confirmed by an archival HST/GHRS spectrum, taken from the HST/MAST archive, and showing a velocity consistent with zero, to within 5 km s^-1 for the N i λ1199.55 Å line from the ISM. This then leads to consistent observed offsets for the 1s–2p transitions of O i and N i in Sco X-1 of 1.0 ± 2.3 mÅ and 4.0 ± 3.0 mÅ, respectively.

Combining both the O i and the N i line for Mrk 509, we obtain a shift of +10.3 ± 4.0 mÅ.

6.4.3. Lines from foreground hot gas absorption

Our spectrum also contains significant lines from hot, ionised gas along the line of sight, in particular from O vi, O vii, O viii, and C vi. Sembach et al. (2003) give a plot of the O vi line profile. This line shows clear high-velocity components at −247 and −143 km s^-1. Unfortunately, no column density of the main Galactic component is given by Sembach et al. (2003). Therefore we have fitted their profile to the sum of 5 Gaussians (Table 7). The column-density weighted centroid for the full line is at −56 km s^-1 in the LSR frame. We estimate an uncertainty of about 8 km s^-1 on this number.

Note that a part of the high-velocity trough in the blue is contaminated by a molecular hydrogen line at 1031.19 Å, which falls at v_LSR = −206 km s^-1. This biases component 1 of our fit for the hot foreground galactic component by about 10 km s^-1. This component in C iv, Si iv, and N v in the COS spectra is at a velocity of v_LSR = −245 km s^-1 (Kriss et al. 2011), and therefore we adopt that value here.

While this velocity can be readily used for predicting the observed wavelength of the O vi 1s–2p line, the uncertainty for other X-ray lines is larger, because the temperatures or ionisation states of the dominant components 1 and 4 may be different. Moreover, O vi is the most highly ionised ion available in the UV band, while most of our X-ray lines from the ISM have a higher degree of ionisation. Therefore, without detailed modelling of all these components, we conservatively assume here for all ions other than O vi an error of 200 km s^-1 on the average outflow, corresponding to the extreme cases that the bulk of the column density is in either component 1 or 5.

The Galactic O vii 1s–2p line at 21.6 Å suffers from some blending by the N vii 1s–3p line from the outflow. That last line, with a predicted equivalent width of 4.5 mÅ, has a relatively high ionisation parameter, log ξ = 1.4. Therefore we do not know exactly what its predicted velocity should be, and we conservatively assume a value of −100 ± 200 km s^-1, spanning the full range of outflow velocities. With that, the predicted wavelength of the N vii outflow line is 21.624 ± 0.014 Å, while the predicted wavelength of the O vii line is 21.5993 ± 0.0006 Å. The predicted wavelength of the blend then is 21.606 ± 0.004 Å.

In Table 8 we list the lines that we use and the derived line shifts Δλ. The weighted average is +1.7 ± 2.6 mÅ when we only include statistical errors, and +2.8 ± 4.9 mÅ if we include the systematic uncertainties.

We make here a final remark. Combining all wavelength scale indicators including the one above (Sect. 6.4.6), we obtain an average shift of −1.6 mÅ. Using this number, the most significant lines of Table 8 (O vii and O viii 1s–2p) clearly show a velocity much closer to the dominant velocity component 4 (+26 km s^-1, Table 7) than to the column-density weighted velocity of the O vi line at −56 km s^-1. This last number is significantly affected by the high-velocity clouds, which apparently are less important for the more highly ionised ions. Therefore, in Table 8 we also list in Col. 3 the wavelength shifts under the assumption that all lines except O vi are dominated by component 4. With that assumption, and taking now only the statistical errors into account, the weighted average shift is −1.7 ± 2.8 mÅ.

6.4.4. Comparison with a sample of RGS spectra of stars

Preliminary analyses (González-Riestra, priv. comm.) of a large sample of RGS spectra of stellar coronae indicate that there is a correlation between the wavelength offsets and the solar aspect angle (SAA) of the satellite (correlation coefficient about 0.6). At a mean SAA of 90°, the average Δλ for RGS1 is 2.3 mÅ, while for RGS2 it is 7.3 mÅ. The scatter for an individual observation is 6 and 7 mÅ, respectively. The difference between RGS1 and RGS2 is fairly consistent with the 6.9 ± 0.7 mÅ that we adopted for the present work. Because our observation campaign on Mrk 509 spans the full visibility window, the average SAA angle is 90°; it varied between 1085 for the first observation to 733 for the last observation. The wavelengths in our spectrum (averages of the original RGS1 wavelengths and the RGS2 wavelengths minus 6.9 mÅ) should then be off by +1.4 ± 6.5 mÅ.

The above estimate can be refined further. It appears that there is also a correlation of the wavelength offset with the temperature parameters T2007 and T2031 of the grating boxes on RGS1 and RGS2, respectively. However, these temperatures do not correlate well with the SAA. Hence, following González-Riestra et al. (priv. comm.) and using their data on 38 spectra of stellar coronae, we obtain the following correlations: where the SAA is expressed in degrees, and the temperatures in °C. The scatter around this relation is 4.1 and 4.8 mÅ for RGS1 and RGS2, respectively. This scatter is significantly smaller than the scatter using only the correlation with SAA mentioned above.

We apply this now to the 10 individual observations of Mrk 509. We have determined the individual line centroids of the six strongest lines (C vi 1s–2p, O viii 1s–2p and 1s–3p, O vii 1s–2p, Ne ix 1s–2p all from the outflow, and O i 1s–2p from the foreground). We then take the difference with the wavelength in the combined observation, and averaged these residuals. We show the results in Table 9 and Fig. 10.

The model gives an improvement: the rms residuals (after correction for the mean statistical errors on the data points of 2.3 mÅ) are reduced from 5.4 to 4.6 mÅ. This is fully consistent with the remaining residuals in the set of coronal spectra discussed above. It is also clear that the remaining residuals are non-statistical, and although on average the model gives an improvement, in some individual cases the opposite is true, for instance for observations 3, 7 and 10.

Keeping this in mind, we can use the model to predict the absolute wavelength scale. We predict absolute offsets Δλ of −2.76 and +0.10 mÅ for RGS1 and RGS2, respectively. In our extraction procedure, we have subtracted 6.9 mÅ from the RGS2 wavelengths; doing the same for the predicted RGS offset, and averaging the offsets for RGS1 and RGS2, the model predicts a total offset of −4.8 mÅ. We only need to apply a small correction for the velocity of the Earth. It can be assumed that the sample of stars has, on average, zero redshift because of the motion of Earth, because usually each source can be observed twice a year with opposite signs of the Doppler shift. Our data of Mrk 509 were taken with a rather extreme Earth velocity of −29.3 km s^-1. At a typical wavelength of 20 Å, this corresponds to a correction factor of 2.0 mÅ. We finally obtain a predicted offset of −2.8 ± 4.5 mÅ  where the error is the average systematic residual found for the sample of coronae mentioned before.

6.4.5. Comparison with Chandra LETGS spectra

We have also measured the wavelengths of the 1s–2p and 1s–3p lines of O vii and O viii in the Chandra LETGS spectrum of Mrk 509 taken 20 days after our last XMM-Newton observation. LETGS has the advantage above RGS that it measures both positive and negative spectral orders, hence allow us to determine the zero-point of the wavelength scale more accurately. Table 10 lists the wavelengths measured by both instruments. The average wavelength offset of RGS relative to LETGS is −6.9 ± 2.4 mÅ. The above number only includes the statistical uncertainties. We have extracted an LETGS spectrum of Capella, taken around the same time as the Mrk 509 observation. The ten strongest lines in the spectrum of Capella between 8–34 Å show a scatter with rms amplitude of 6 mÅ compared to the laboratory wavelengths, but no systematic offset larger than 2 mÅ. We therefore add a systematic uncertainty of 6/ mÅ in quadrature to the statistical uncertainty of the average Δλ for the four lines in Mrk 509 that we use here.

We only need to make a small correction of −0.4 mÅ corresponding to the higher speed of Earth away from Mrk 509 (−29.3 km s^-1) during the XMM-Newton observation as compared to the Chandra observation (−23.3 km s^-1). This all leads to an offset of −7.3 ± 3.8 mÅ.

6.4.6. Summary of results on the wavelength scale

We now summarise our findings from the previous sections (Table 11). All estimates agree within the error bars, except for the lines from the neutral foreground gas, which is off by about 3σ. This number is dominated by the fairly complex O i blend. Possibly there are larger systematic uncertainties on this number than we assumed. The major factor here is the blending with Ca xv. If we assume that the equivalent width of the contaminating Ca xv line is zero, the predicted blend centroid decreases by 4.4 mÅ, towards better agreement with the other indicators. The O iv blend has a much smaller impact, because its centroid is already close to the centroid of the Galactic O i line.

We note that a lower value of the predicted equivalent width of the Ca xv line may be caused by a somewhat lower Ca abundance; alternatively, if the density is higher than about 10¹⁶ m^-3, the population of Ca xv in the ground state decreases significantly (Bhatia & Doschek 1993), leading also to a reduced equivalent width of this absorption line.

All five methods used here have their advantages and disadvantages, and we have done the best we can to estimate the uncertainties for each of the methods. The five methods are for the major part statistically independent, and this allows us to use a weighted average for the final wavelength correction. Ignoring the blending by Ca xv to the Galactic O i line, and using +26 km s^-1 for the outflow velocity of the hot Galactic gas, this wavelength shift is −1.8 ± 1.8 mÅ. Within their statistical uncertainties, all individual indicators are consistent with this value, except for the lines from foreground neutral gas, which are off by 8.3 ± 3.8 mÅ. We do not have an explanation for that deviation.

In our analysis, we apply the wavelength corrections by shifting the wavelength grids in our program rgs_fmat.

7. AGN outflow dynamics

To show how the optimised wavelength calibration improves the scientific results, we give a brief analysis of the velocity structure of the AGN outflow as measured in X-rays through two Rydberg series of oxygen ions. A full analysis of the spectrum is given by Detmers et al. (2011).

Using the measured outflow velocities for individual lines, we can immediately derive some interesting astrophysical conclusions (Fig. 11).

Most of the lines are in between the high-velocity components 1–3 and the lower velocity components 4–7 as found in the O vi lines (Kriss et al. 2000). This shows that there are multiple velocity components present in the X-ray absorption lines. At longer wavelengths these lines are closer to the low velocity components, showing that there is a higher column density in those components. However, for the shortest wavelength lines, in particular the 1s–2p transitions of Ne ix and Ne x, the centroid is in the range of the high-velocity components 1–3. Because the shorter wavelength lines originate from ions with a higher ionisation parameter, there is a correlation between the ionisation parameter and the outflow velocity.

The relative deviation of the 17.8 Å line (O vii 1s–4p) with respect to the other lines, including the 1s–2p and 1s–3p lines of the same ion, can be easily explained: according to preliminary modelling, the optical depth at line centre for the 1s–2p, 1s–3p and 1s–4p lines of O vii is 8, 1.4 and 0.5, respectively. Thus, if the O vii column density of the low velocity component is higher than that of the high-velocity component, the strongest low-velocity line components (like 1s–2p, 1s–3p) will saturate while the high-velocity components of these transitions will not saturate, hence effectively these strong lines will be blueshifted relative to the weaker 1s–4p line, consistent with what is observed.

We have elaborated on the line diagnostics for two of the most important ions in the outflow of Mrk 509, O viii and O vii, to see if we can detect both velocity components.

For O viii we have measured velocities and equivalent widths for the 1s–2p to 1s–5p lines (Table 3). We have tried different models. The first model consists of a single Gaussian absorption line, with free centroid, width and column density. Model 2 has two Gaussian components, with the centroids frozen to the values given in Table 4 as derived from the O vi troughs, but with the widths and column densities free parameters. The results are shown in Table 12 and Fig. 12. Both models give a good fit and cannot be distinguished on statistical grounds. The total column densities for both models are the same within their error bars, and also the sum of the width of the Gaussians for Model 2 is consistent with the width of the single Gaussian for Model 1. Apparently, as far as the determination of total line width and column density is concerned, the detailed column density distribution as a function of velocity is not very important. This can be understood because the strongest line with the best statistics (1s–2p) has a high optical depth (12.7 and 4.9 for the two components in model 2). The line core is thus black, the line shape rather “rectangular” and to lowest order approximation the width of the line is given by the sum parts of the spectrum where the core is black.

For Model 2 the total width for each of the components agrees well with the total width of components 1–3 and 4–7 of the O vi line. Thus, O viii originates from multiple components, but the fastest outflow components have a higher column density than the slower outflow components.

We have performed a similar analysis for the 1s–2p to 1s–5p lines of O vii (Table 13). The velocity broadening is consistent with what we obtained for O viii, although there is somewhat larger uncertainty. For O vii, however, the low velocity components 4–7 have clearly the highest column density, similar to what was found for O vi (Kriss et al. 2000).

8. Statistics

Local continuum fits to RGS spectra in spectral bands that have few lines yield in general χ² values that are too low. This is an artefact of the SAS procedures related to the rebinning of the data. Data have to be binned from the detector pixel grid to the fixed wavelength grid that we use in our analysis. However, the bin boundaries of both grids do not match. As a consequence of this process, the natural Poissonian fluctuations on the spectrum as a function of detector pixel are distributed over the wavelength bin coinciding most with the detector pixel and its neighbours. In addition, there is a small smoothing effect caused by pointing fluctuations of the satellite. Owing to this tempering of the Poissonian fluctuations, χ² values will be lower for a “perfect” spectral fit.

We have quantified the effect by fitting a linear model F_λ = a + bλ to fluxed spectra in relatively line-poor spectral regions, in 1 Å wide bins in the 7–10, 17–18 and 35–38 Å ranges. The median reduced χ² is 0.72, with a 67% confidence range between 0.65 and 0.79. Over the full 7–38 Å range, we have at least 400 counts (and up to 3000 counts) in each bin of ~0.02 Å width. Therefore, the Poissonian distribution is well approximated by the Gaussian distribution, and the usage of chi-squared statistics is well justified.

The rebinning process conserves the number of counts, hence the nominal error bars (square root of the number of counts) are properly determined. The lower reduced χ² is caused by the smoothing effect on the data. For correct inferences about the spectrum, such a bias in χ² is not appropriate. Because we cannot change the observed flux values, we opt to multiply the nominal errors on the fluxes by to obtain acceptable values for ideal data and models.

As a final check we have made a more robust estimate of the effect. For RGS1, the wavelength channel number c_λ as a function of dispersion angle channel c_β can be approximated well with (20)Locally, if the spectrum on the β-grid is given by a set of random variables X_j and the spectrum on the wavelength grid is given by Y_i, we have (21)for a small number of bins j around the corresponding β-channel given by (20). It is important to note here that while the X_j variables are statistically independent, the Y_i variables in general are not. Typically, the f_ij factors are proportional to the overlap of wavelength bin j with β-bin i. The clue is now that from elementary statistical considerations we have for the expected value of Y_i and the variance For instance, if the single β-bins j were exactly redistributed over two adjacent λ-bins, there would be two f-values, each = 0.5, and we have Var [Y_i] /E [Y_i] = 0.5Var [X_j] /E [X_j] , where the factor 0.5 comes from 0.5² + 0.5². We exactly made this calculation exactly using the grids related through (20), and show our result in Fig. 13.

9. Conclusion

9.1. Wavelength scale and spectral resolution

By using five different methods, we have obtained a consistent wavelength scale with an accuracy of about 1.8 mÅ, four times better than the nominal wavelength accuracy of RGS. The one slightly deviant diagnostic that remains is the neutral gas in our galaxy (O i and N i), which is off from the mean by 8.8 ± 3.8 mÅ (2.3σ). The systematics from blending (for O i) and the lack of sufficiently reliable RGS data for N i may result into slightly larger uncertainties than we have estimated.

Our work shows that in cases like ours a dedicated effort to improve the wavelength scale is rewarding. As we showed, at this level of accuracy, Doppler shifts caused by the orbit of Earth around the Sun become measurable. The improved accuracy also allows us to use X-ray lines as tools to study the dynamics of the outflow and the Galactic foreground at a level beyond the nominal resolution of the instrument (see Sect. 7).

9.2. Effective area corrections

The corrections that we have derived for the relative effective area of both RGS detectors as well as for the relative effective area of both spectral orders are at the few percent level but should be taken into account for the high statistical quality spectrum that we consider here. This is because there are data gaps for each individual RGS, but each time at a different wavelength range. Without this correction, there would be artificial jumps at the boundaries of these regions of the order of a few percent. This leads any spectral fitting package to attempt to compensate for this by enhancing or decreasing astrophysical absorption lines or edges, in the hunt for the lowest χ² value.

A disadvantage of our method is of course that the “true” effective area and hence absolute flux is known only to within a few percent, and the deviations are wavelength-dependent. In our subsequent spectral fitting (Detmers et al. 2011) we compensate for this by using a spline for the AGN continuum. This is clearly justified, as one might worry about the intrinsic applicability of e.g. pure power-law or blackbody components when the flux can be measured down to the percent level accuracy. The “true” continua may contain all kinds of features including e.g. relativistically broadened lines (e.g. Branduardi-Raymont et al. 2001; Sako et al. 2003) that can be discerned only in the highest-accuracy data, such as the present ones.

9.3. Combination of spectra

The reduction of the combined response matrices of our 40 spectra with 2 Gb of memory into a single spectrum and response file of only 8 Mb, a reduction by more than two orders of magnitude, allows us to make complex spectral fits. In the most sophisticated models used by Detmers et al. (2011) the number of free parameters approaches 100, and without this efficient response matrix the analysis would have been impossible.

9.4. Final remarks

Through the high statistical quality of our data, the binning problems mentioned in Sect. 8 became apparent. In many RGS spectra that have been obtained during the lifetime of XMM-Newton, the statistical quality of the spectra is not good enough to detect this effect. Still, this could imply that when users obtain a fit with reduced chi-squared of order unity, the fit is formally not acceptable because an ideal fit would yield a smaller reduced chi-squared. Clearly, more work needs to be done here to alleviate the problem, but the discretisation of wavelength in the detector corresponding to the natural size of the CCD pixels limits the possibilities. Even in the best observations there is always some small scatter in the satellite pointing, making some rebinning necessary.

Overall, however, the calibration of RGS is at a fairly advanced stage. Only for exceptionally good spectra like the Mrk 509 spectrum that we discussed here the sophisticated tools developed in this paper are needed, although our improvements can be beneficial for the analysis of weaker sources as well. Our tools are publicly available as auxiliary programmes (rgsfluxcombine, rgs_fmat) in the SPEX distribution.

---

Acknowledgments

We thank Charo González-Riestra (ESA) for providing the data on the RGS spectra of stars in the context of the wavelength scale calibration. This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). SRON is supported financially by NWO, The Netherlands Organization for Scientific Research. K. Steenbrugge acknowledges the support of Comité Mixto ESO – Gobierno de Chile. Ehud Behar was supported by a grant from the ISF. G. Kriss gratefully acknowledges support from NASA/XMM-Newton Guest Investigator grant NNX09AR01G. M. Mehdipour acknowledges the support of a Ph.D. studentship awarded by the UK Science and Technology Facilities Council (STFC). P.-O. Petrucci acknowledges financial support from CNES and the French GDR PCHE. Ponti acknowledges support via an EU Marie Curie Intra-European Fellowship under contract No. FP7-PEOPLE-2009-IEF-254279. G. Ponti and S. Bianchi acknowledge financial support from contract ASI-INAF No. I/088/06/0.

References

All Tables

All Figures