# Concordance: In-flight Calibration of X-ray Telescopes without Absolute
References

Herman L. Marshall  Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, 77 Massachusetts Ave., Cambridge, MA
02139, USA  Yang Chen  U. Michigan, Ann Arbor, MI 48109, USA  Jeremy J. Drake
Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA
02138, USA  Matteo Guainazzi  ESTEC, Keplerlaan 1 2201AZ Noordwijk, The
Netherlands  Vinay L. Kashyap  Harvard-Smithsonian Center for Astrophysics, 60
Garden St., Cambridge, MA 02138, USA  Xiao-Li Meng  Harvard University,
Cambridge, MA, 02138, USA  Paul P. Plucinsky  Harvard-Smithsonian Center for
Astrophysics, 60 Garden St., Cambridge, MA 02138, USA  Peter Ratzlaff
Harvard-Smithsonian Center for Astrophysics, 60 Garden St., Cambridge, MA
02138, USA  David A. van Dyk  Imperial College, London, UK SW7 2AZ  Xufei Wang
Harvard University, Cambridge, MA, 02138, USA

(Accepted August 30, 2021)

###### Abstract

We describe a process for cross-calibrating the effective areas of X-ray
telescopes that observe common targets. The targets are not assumed to be
“standard candles” in the classic sense, in that we assume that the source
fluxes have well-defined, but a priori unknown values. Using a technique
developed by Chen et al. (2019) that involves a statistical method called
shrinkage estimation, we determine effective area correction factors for each
instrument that brings estimated fluxes into the best agreement, consistent
with prior knowledge of their effective areas. We expand the technique to
allow unique priors on systematic uncertainties in effective areas for each
X-ray astronomy instrument and to allow correlations between effective areas
in different energy bands. We demonstrate the method with several data sets
from various X-ray telescopes.

††journal: AJ††facilities: CXO(ACIS), CXO(HRC), CXO(HETGS), CXO(LETGS), XMM-
Newton(pn), XMM-Newton(MOS)††software: Concordance
https://github.com/astrostat/Concordance, ciao (Fruscione et al., 2006), SAS
(Gabriel et al., 2004), MCCal (Drake et al., 2006), PINTofALE (Kashyap &
Drake, 2000)

##  1 Introduction

We address a perennial issue in instrument performance, when estimated fluxes
using two or more instruments disagree. While many instruments can safely rely
on calibration traceable to established standards, the performance of a space-
based telescope usually cannot be recalibrated, as the instruments are not
returned to the lab. Furthermore, space-based instruments may be affected by
the physical rigors of launch into space and instrument performance can change
with time due to gas leakage, filter deterioration, component failure,
contamination buildup and other reasons. Without absolute standards that may
be observed while in space, astronomers generally resort to the use of
secondary, astronomical standards that can be observed while operating.
Nonvariable sources are generally chosen as secondary standards so that they
can be reused by the telescope team and observed by others. For example, in
the 2-8 keV energy range, the Crab Nebula was frequently used as a standard,
especially the observation by Toor & Seward (1974); for a recent use of the
Crab Nebula for cross-calibrating X-ray telescopes, see Kirsch et al. (2005)
for a comparison between many missions and Madsen et al. (2017b) for a case of
assessing the NuSTAR telescope effective model. However, most astronomical
sources vary, even the Crab Nebula, at levels detectable by current
instruments, necessitating another approach: joint in-flight observations for
cross-calibration of instruments.

Our work is set upon the foundation established by the International
Astronomical Consortium for High Energy Calibration (IACHEC). The IACHEC was
formed primarily to assist X-ray telescope teams who cross-calibrate
instruments and to understand the sources used for this purpose. See recent
IACHEC reports for summaries of recent activity (Madsen et al., 2019, 2020).
All of the IACHEC working groups address issues of cross calibration of X-ray
telescopes and often find discrepancies between results for the same source.
Examples of IACHEC work include observations of the supernova remnants (SNRs)
G21.5 (Tsujimoto et al., 2011) and 1E 0102$-$7219 (Plucinsky et al., 2017),
spectra of galaxy clusters (Nevalainen et al., 2010; Kettula et al., 2013;
Schellenberger et al., 2015), spectra of white dwarfs and isolated neutron
stars (Beuermann et al., 2006), and simultaneous observations of active
galaxies such as 3C 273 and PKS 2155$-$304 (Ishida et al., 2011; Madsen et
al., 2017a). The studies had a common problem: assessing how much any
particular instrument’s effective area should be adjusted so that the
measurements might agree. Ordinary weighting of measurements based on photon
counting statistics would give the observations with large effective areas and
exposure times the greatest influence on the result but without consideration
of possible systematic errors. The overarching goal of IACHEC is thus to bring
the competing adjustments to instrumental effective areas into concordance,
while simultaneously including prior knowledge of possible systematic errors.

Here, we further develop the “shrinkage” method pioneered by Chen et al.
(2019), hereafter referred to as Paper I, to compute objective corrections to
effective areas of several high-energy instruments. See Section 2.1 for the
model setup and notation. The method is extended to account for systematic
uncertainties specific to each instrument and incorporate systematic
correlations across passbands. We present the method here in some generality,
without the details of the mathematical machinery presented in Paper I
(Section 2.2), along with extensions added (Section 2.3), and apply it to a
variety of datasets (Section 3). We present the results in Section 4 along
with simulation studies that aim to validate the method, and discuss the next
steps in Section 5. We note that while we examine the case of X-ray telescopes
in this paper, the method is extendable to most types of telescopes.

##  2 Method

###  2.1 Background

We start with an idealized calibration data set where $M$ objects are observed
by each of $N$ instruments, obtaining photon counts, $c_{ij}$, where
$j=1,\ldots,M$ indexes the sources and $i=1,\ldots,N$ indexes the
instruments.111We consider each unique combination of telescope and detector
with any filters or gratings to be an “instrument”. Each observation is
described by a set of observational parameters (e.g., exposure time) that we
encapsulate in a matrix $\mathbf{T}=\\{T_{ij}\\}$.

Denoting the true effective area of instrument $i$ by $A_{i}$ and the true
flux of source $j$ by $F_{j}$, the expected counts $C_{ij}$ for each
object/detector combination is given by

$C_{ij}=T_{ij}A_{i}F_{j},\quad 1\leq i\leq N,\quad 1\leq j\leq M,$ (1)

where $T_{ij}$ has units of seconds $\times$ counts per photon, $F_{j}$ has
units of photons per unit area per second, and $A_{i}$ has units of area. We
assume for now that the true exposure factors have negligible error and equal
the observed values, i.e., $T_{ij}=t_{ij}$. We follow the notation of Paper I
by using lower case to indicate measured quantities and upper case to indicate
the “true” values to be estimated. Note that in fact, the multiplicative
constant $T_{ij}$ contains not only the exposure time but also other factors
that can be calculated precisely for any given observation.

We aim to estimate $F_{j}$ using the calibration data, i.e., the observed
counts $c_{ij}$, and an external (prior) estimated effective area, $a_{i}$ for
each instrument. A naive procedure simply substitutes each quantity in Eq. (1)
with its measured counterpart, to obtain the “estimating equation”,

$c_{ij}=t_{ij}a_{i}f_{ij}$ (2)

and solving Eq. (2) for $f_{ij}$ for each instrument-source combination, thus
yielding $N$ different estimates of each flux. The resulting estimator
$f_{ij}=c_{ij}/(a_{i}t_{ij})$ is a ratio estimator, which is known in
statistical literature to be both seriously biased and highly variable. More
precisely, the variability in the denominator can cause large uncertainty
(considering dividing by a value close to zero). Furthermore, the average,
$\frac{1}{N}\sum_{i=1}^{N}f_{ij}$, is a biased estimator of $F_{j}$. We return
to the issue associated with ratio estimators in §4.4.1.

We can do much better by analyzing all data together using more principled and
sophisticated statistical methods, such as the one given in Paper I. We can
then achieve our goal to obtain better estimates of the instrumental effective
areas, $A_{i}$, that bring our flux estimates, $f_{ij}$, closest to the
$F_{j}$ in some strict statistical sense. Fig. 1 gives a schematic
representation of our goal.

Figure 1:  The goal of Concordance illustrated schematically. This schematic
supposed $3$ sources and $3$ instruments, and plots expected count rates,
$C_{ij}/T_{ij}$, on the vertical axis and plots true and estimated fluxes on
the horizontal axis. The three instruments are represented by different
symbols, with solid symbols representing the naive estimates of $F_{j}$, i.e.,
$f_{ij}=c_{ij}/(T_{ij}a_{i})$, and open symbols representing the same, but
with $a_{i}$ replaced by the true effective area, $A_{i}$. The estimates are
systematically biased because the $a_{i}$ are inaccurate estimates of $A_{i}$;
compare the solid and dashed lines. (In the plot EA is used as an abbreviation
for effective area.) Our aim is to estimate the $A_{i}$ values that yield best
agreement between instruments for each source.

In practice, estimating a flux from an observation is not as simple an
operation as merely counting events. A strictly correct approach takes into
consideration the response of the instrument, especially when attempting to
measure the flux in a specific bandpass. It is often necessary to use a
forward-folding method, such as implemented in xspec, that takes a count
spectrum as input and effective area and response functions as externally (and
accurately) provided. Here, we approximate this process using Eq. (1) because
the flux (as defined here) is often a robustly computed quantity for any given
observation. In fact, we do not actually compute $c_{ij}$ but a proxy for
$c_{ij}$ as described below.

We now carefully examine the flux measurement process to justify using the
simplistic model represented by Eq. (1). Source $j$ is assumed to have a
photon spectrum $n_{E}(\Theta_{j})$ in units of photons per unit area per unit
energy at energy $E$. Each $\Theta_{j}$ represents a vector of spectral
parameters for a source, which we assume to include an overall normalization
$n_{j}$ with the same units as $n_{E}$, so $q(E;\theta_{j})\equiv
n_{E}(\Theta_{j})/n_{j}$ defines the spectral shape as a function of the
remaining spectral parameters, $\theta_{j}$. In the case of a spectrometer
measuring an emission line, the bandpass may be so narrow that $q$ is well
approximated by a delta function, in which case $F_{j}=n_{j}$. Otherwise, the
flux in a band from $E_{1}$ to $E_{2}$ is

$F_{j}=\int_{E_{1}}^{E_{2}}n_{E}(\Theta_{j})dE=n_{j}\int_{E_{1}}^{E_{2}}q(E;\theta_{j})dE$
(3)

giving

$n_{j}=\frac{F_{j}}{\int_{E_{1}}^{E_{2}}q(E;\theta_{j})dE}\equiv\frac{F_{j}}{\tilde{q}_{j}\Delta
E}$ (4)

where $\Delta E=E_{2}-E_{1}$.

For most X-ray telescopes, the instrumental effective area can be separated
into two parts. There is the geometric area of the optics, i.e., $A_{i}^{g}$,
with losses due to mechanical obscuration, reflections, and transmissions of
the optics (including filters) given by $r_{i}(E)$. The other part is due to
the quantum efficiency of the detector, $Q(E)$, which gives probability of a
detecting a photon. We characterize the (true) effective area of instrument
$i$ as a function of energy in the bandpass of interest by
$\tilde{A}_{i}(E)=A_{i}^{g}r_{i}(E)Q_{i}(E)\equiv A_{i}\alpha_{i}(E)$, where
$\alpha_{i}(E)$ now gives the shape of the effective area in the band and is
defined such that its integral over the band is unity, i.e.,

$\alpha_{i}(E)=\frac{r_{i}(E)Q_{i}(E)}{\int_{E_{1}}^{E_{2}}r_{i}(E)Q_{i}(E)dE}.$
(5)

Consequently,

$A_{i}=A_{i}^{g}\cdot\int_{E_{1}}^{E_{2}}r_{i}(E)Q_{i}(E)dE$ (6)

is the scale of the effective area and does not depend on $E$. We take this
approach because the model of the effective area through the band is generally
better known than its absolute value.

The detector provides counts in $K$ channels that are related to the true
photon energy, $E$, via a response function $\Phi_{k}(E)$, where

$\sum_{k=1}^{K}\Phi_{k}(E)=1.$ (7)

A given observation has an accurately determined exposure time $t_{ij}$ that
sets the expected count in channel $k$:

$\displaystyle C_{ijk}$ $\displaystyle=$ $\displaystyle
t_{ij}\int\tilde{A_{i}}(E)n_{E}(\Theta_{j})\Phi_{k}(E)dE$ (8) $\displaystyle=$
$\displaystyle t_{ij}A_{i}n_{j}\int\alpha_{i}(E)q(E;\theta_{j})\Phi_{k}(E)dE$
(9) $\displaystyle=$ $\displaystyle
t_{ij}A_{i}F_{j}\frac{\int\alpha_{i}(E)q(E;\theta_{j})\Phi_{k}(E)dE}{\tilde{q}_{j}\Delta
E},$ (10) $\displaystyle\equiv$ $\displaystyle\tilde{T}_{ijk}A_{i}F_{j},$ (11)

where Eq. (10) follows from Eq. (4). It is important to point out that Eq. (8)
is well-known in the astronomy literature but Eq. (10) is an innovative way of
simplifying the expressions, which is essential for tackling the current
problem. Finally, the counts that are associated with the energy band
$(E_{1},E_{2})$ are mostly in the a range of channels $(k_{1},k_{2})$, which
is chosen so that we have a reliable estimate of flux in the band of interest,
giving the expected count relevant to a particular flux measurement

$C_{ij}=\sum_{k=k_{1}}^{k_{2}}C_{ijk}=A_{i}F_{j}\sum_{k=k_{1}}^{k_{2}}\tilde{T}_{ijk}\equiv
T_{ij}A_{i}F_{j},$ (12)

which defines $T_{ij}$ in terms of $t_{ij}$, a normalized sum over the
response function, and the two shape functions. This is the formal basis for
Eq. (1).

In actual data analysis, the observed counts, $c_{ijk}$, not the expected
counts, $C_{ijk}$, are input into the iterative routine that estimates the
fluxes, $f_{ij}$, and, again, the estimated effective area, $a_{i}$, must be
used because the true effective area is unknown. Thus, in analogy to Eq. (2),
we replace $C_{ij}$, $A_{i}$, and $F_{j}$ in Eq. (12) with their observed
counterparts to obtain the estimating equation

$c_{ij}=\sum_{k=k_{1}}^{k_{2}}c_{ijk}=a_{i}f_{ij}T_{ij}=a_{i}f_{ij}\sum_{k=k_{1}}^{k_{2}}\tilde{T}_{ijk}.$
(13)

Finally, naive, instrument-specific flux estimates are obtained by solving Eq.
(13),

$f_{ij}=\frac{c_{ij}}{T_{ij}a_{i}},$ (14)

which is essentially a full derivation of the ratio estimator.

It is important to our analysis that the $T_{ij}$ values be independent of
$A_{i}$ and $F_{j}$. There are two circumstances that may cause a problem in
this regard: 1) when there is some nonlinearity of the detector response, such
as pileup, where $Q_{i}(E)$ depends on the source flux, and 2) when
$\alpha_{i}(E)$ is highly nonlocal and the bandpass of integration in Eq. (8)
is large, involving portions of the spectrum where $q(E;\theta_{j})$ or
$\alpha_{i}(E)$ are poorly determined. By choosing instrument data in the
linear regime, avoiding pileup, and restricting the bandpasses of interest, we
mitigate these potential issues.

The goal of our statistical modeling is to determine the best estimates of
$A_{i}$ and $F_{j}$ consistent with the data, $c_{ij}$ with uncertainties
$\sigma_{ij}$. If all the effective areas were precisely correct, i.e. if
$A_{i}=a_{i}$, then we could estimate the $F_{j}$ via a relatively trivial
regression model. However, we know that there are systematic errors in our
estimated effective areas because observations show that the $f_{ij}$ do not
cluster about any given value within their individual uncertainties. We
addressed this problem in Paper I by introducing estimates of the systematic
errors on the estimated effective areas and applying a statistical method
called shrinkage estimation.

###  2.2 Statistical Model and Shrinkage Estimates.

Paper I proposes a linear regression model for the log count rates, denoted
$y_{ij}\equiv\log(c_{ij}/T_{ij})$, such that

$y_{ij}=B_{i}+G_{j}-\frac{1}{2\sigma_{i}^{2}}+e_{ij},$ (15)

where the $B_{i}\equiv\log A_{i}$, and the $G_{j}\equiv\log F_{j}$ are the
quantities to be estimated, and the $e_{ij}$ are noise terms that are assumed
normally distributed with mean $0$ and (known) variance $\sigma_{i}^{2}$. The
$-\frac{1}{2\sigma_{i}^{2}}$ term in Eq. (15) is a half-variance correction
that is included to maintain the multiplicative mean modeling in Eq. (12).
This correction ensures that $\hbox{E}(c_{ij})=C_{ij}=T_{ij}A_{i}F_{j}$
because if $\log x\sim\hbox{N}(\mu,v)$, then $\hbox{E}(x)=e^{0.5v+\mu}$.

Paper I adopts a Bayesian hierarchical modeling approach to estimate the
unknown quantities $B_{i}$, $G_{j}$, and $\sigma_{i}^{2}$. This involves
setting independent Gaussian prior distributions for the $B_{i}$, with prior
mean $b_{i}=\log a_{i}$ and prior variance $\tau_{i}^{2}$ and setting
independent flat priors for the $G_{i}$; and setting independent conjugate
priors for the $\sigma_{i}^{2}$. Paper I goes on to show how a Hamiltonian
Monte Carlo algorithm can be used to obtain a sample from the joint posterior
distribution of all unknown quantities in Model (15) and cross-checks this
computation with a blocked Gibbs sampler. The resulting Monte Carlo sample can
be transformed back to the effective areas and source fluxes on their original
scale to obtain their posterior distributions, estimates, and error bars.

A Bayesian perspective allows us to combine multiple sources of information –
in this case the information from the calibration data, the $y_{ij}$, and from
the prior estimates of the effective areas, the $a_{i}$. By replacing the
estimated effective areas with prior distributions that reflect their
uncertainty, we are able to update the estimated effective areas and their
error bars in light of the calibration data. This approach provides improved
estimates (e.g., in terms of mean squared error) of the effective areas and
the fluxes simultaneously. As a result, we obtain estimates of the sources’
true fluxes that combine the instrument-specific estimates in a statistically
principled manner.

The updated estimates of the effective areas are weighted averages of their
priors and the best values based on the current calibration data. In the
context of Model (15), we work on the log scale. The estimates of ${B}_{i}$
and ${G}_{j}$ are given by

$\widehat{B}_{i}=W_{i}(\bar{y}_{i\cdot}^{\prime}-\bar{G}_{i})+(1-W_{i})b_{i}\quad\text{
and }\quad\widehat{G}_{j}=\bar{y}_{\cdot j}^{\prime}-\bar{B}_{j},$ (16)

where

$\bar{y}_{i\cdot}^{\prime}=\frac{\sum_{j=1}^{M}y^{\prime}_{ij}\sigma^{-2}_{i}}{\sum_{j=1}^{M}\sigma^{-2}_{i}},\quad\bar{y}_{\cdot
j}^{\prime}=\frac{\sum_{i=1}^{N}y^{\prime}_{ij}\sigma^{-2}_{i}}{\sum_{i=1}^{N}\sigma^{-2}_{i}},\quad\bar{G}_{i}=\frac{\sum_{j=1}^{M}\widehat{G}_{j}\sigma^{-2}_{i}}{\sum_{j=1}^{M}\sigma^{-2}_{i}},\quad\bar{B}_{j}=\frac{\sum_{i=1}^{N}\widehat{B}_{i}\sigma^{-2}_{i}}{\sum_{i=1}^{N}\sigma^{-2}_{i}},$
(17)

$y^{\prime}_{ij}={y}_{ij}+0.5\sigma^{2}_{i}$, and the weights are given by
$W_{i}=M\sigma^{-2}_{i}/(\tau^{-2}_{i}+M\sigma^{-2}_{i})$. If the prior
estimate of the effective area of a particular instrument is very precise
relative to its calibration data, i.e., $\tau_{i}^{2}\ll\sigma_{i}^{2}/M$,
then $W_{i}\approx 0$, and the updated estimate of that instrument’s effective
area is dominated by the prior estimate, resulting in $\widehat{B}_{i}\approx
b_{i}$. In contrast, if the calibration data are much more precise, then the
weights are near unity and the updated estimate of the effective area is
dominated by the calibration data, giving
$\widehat{B}_{i}\approx\bar{y}_{i\cdot}^{\prime}-\bar{G}_{i}$.

The estimates $\widehat{B}_{i}$ in Eq. (16) are often called “shrinkage
estimates” due to their historical use for “shrinking” several estimates
together toward a common prior mean (Efron & Morris, 1972, 1973) when, for
example, the $b_{i}$ are all the same. Because the prior means, $b_{i}$, are
different for different instruments, the $\widehat{B}_{i}$ are simply a
sensible combination of prior knowledge captured by $b_{i}$ and data
represented by $\bar{y}_{i\cdot}^{\prime}$, weighted by their respective
precisions, which are the reciprocals of their variances (assuming the
$\sigma_{i}$ are known). This combination allows our model to weigh the prior
estimate of the effective area for a given instrument against deviations
between the observed fluxes of the same sources from different instruments and
ultimately to obtain the joint estimates of the true fluxes and effective
areas that are most consistent with the calibration data and the priors on the
effective areas.

Paper I further describes how to handle the case where all sources are not
observed by all instruments and presents a robust version of Model (15) that
allows for outliers among the measured fluxes (or source counts) by replacing
the $\log$ Normal error model with a $\log t$ error model.

###  2.3 Extensions of the Model

Paper I proposes modeling calibration data using Model (15) and its
extensions, derived computational methods for fitting these models, and
validated their statistical properties. However, application of Model (15) to
real data requires relaxing some of its basic assumptions. Here, we illustrate
how this is accomplished via two extensions to the method.

####  2.3.1 Heterogeneous Uncertainties in Effective Area Priors

IACHEC scientists recognize that the quality of ground-based calibrations
varies significantly from instrument to instrument, resulting in perceived
differences in the reliabilities of the estimated effective areas. The
formalism laid out in Paper I allows for instrument specific prior
distributions for for the $B_{i}$, as explained in Section 2.2, given by
Gaussian distributions with instrument-specific variances $\tau_{i}^{2}$. In
the numerical examples in Paper I, however, the $\tau_{i}^{2}$ were set to
$\tau^{2}$ for each $i$, essentially assuming that all modelled calibrations
are equally uncertain in percentage terms. Here we allow for heterogeneous
$\tau^{2}_{i}$ values, as covered by the theory given in Paper I. At IACHEC
meetings in 2017, 2018, and 2019 (Madsen et al., 2019, 2020), we asked
instrument calibration scientists to specify values of $\tau_{i}$ for their
instruments in each of a specific set of bandpasses. These values are given in
Tables 1 and 2; instruments with significant effective area below 1 keV appear
in Table 1 and other instruments appear in Table 2.

Of course, in practice it is difficult even for experts to quantify the
$\tau_{i}$ precisely. Thus, it is important that we examine the sensitivity of
our results to the specified values. Often, however, there is a body of
experience and expert knowledge on the reliability of ground-based standards
that allows rough estimation of systematic errors.

Table 1: Effective Area Uncertainty Priors ($\tau_{i}$)aaThe $\tau_{i}$ values are given as percentages. The ellipses indicate bandpasses where the instrument has an negligible effective area. | Energy Bands (keV)
---|---
Instrument | 0.15-0.33 | 0.33-0.54 | 0.54-0.8 | 0.8-1.2 | 1.2-1.8 | 1.8-2.2 | 2.2-3.5 | 3.5-5.5 | 5.5-10
Astrosat SXT |  $\cdots$ | 15 | 15 | 10 | 10 | 10 | 10 | 10 | 10
Chandra ACIS | 3 | 3 | 3 | 3 | 2.6 | 3.3 | 3.3 | 4.9 | 5
Chandra HETGS |  $\cdots$ |  $\cdots$ | 10 | 5 | 4 | 4 | 4 | 5 | 7
Chandra LETGS | 5 | 7 | 7 | 7 | 7 | 7 | 7 | 10 | 10
ROSAT PSPC | 10 | 10 | 10 | 10 | 10 | 10 |  $\cdots$ |  $\cdots$ |  $\cdots$
Suzaku XIS1 |  $\cdots$ | 20 | 15 | 10 | 10 | 15 | 5 | 5 | 5
Suzaku XIS0,2,3 |  $\cdots$ |  $\cdots$ | 15 | 10 | 10 | 15 | 5 | 5 | 5
Swift PC/WT |  $\cdots$ | 15 | 10 | 7.5 | 7.5 | 10 | 5 | 5 | 5
XMM MOS1,2 | 20 | 10 | 6 | 6 | 6 | 6 | 6 | 6 | 10
XMM pn | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 2 | 3
XMM RGS |  $\cdots$ | 8 | 5 | 5 | 5 |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$
Table 2: Effective Area Uncertainty Priors ($\tau_{i}$)aaThe $\tau_{i}$ values are given as percentages. | Energy Bands (keV)
---|---
Instrument | 2.2-3.5 | 3.5-5.5 | 5.5-10 | 15-25 | 25-50 | 50-100 | 100-300
Astrosat CZTI |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ 20 | 20 | 20 | 25
Astrosat LAXPC |  $\cdots$ | 15 | 15 | 15 | 15 | 20 |  $\cdots$
INTEGRAL IBIS |  $\cdots$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | 8 | 15 | 20
INTEGRAL SPI |  $\cdots$ |  |  $\cdots$ |  $\cdots$ | 5 | 5 | 5
NuSTAR |  $\cdots$ | 4 | 3 | 3 | 15 | 20 |  $\cdots$
RXTE PCA | 5 | 10 | 3 | 3 | 10 | 50 |  $\cdots$
RXTE HEXTE |  $\cdots$ |  $\cdots$ |  $\cdots$ | 5 | 5 | 5 |  $\cdots$
Suzaku HXD |  $\cdots$ |  $\cdots$ |  $\cdots$ | 20 | 20 | 20 | 20
Swift BAT |  $\cdots$ |  $\cdots$ |  $\cdots$ | 15 | 4 | 4 | 12

####  2.3.2 Prior Correlations among Effective Areas

The second extension allows correlations between the effective areas in
different energy bands for each instrument. In Paper I, we treated different
energy bands as separate (independent) instruments, while in reality their
effective areas can be strongly correlated. Continuing to work on the log
relative scale given in Eq. (23), we denote the effective areas as a function
of the energy band $\mathcal{E}=[E_{1},E_{2}]$ by
${B}(\mathcal{E},\vec{\xi})=\log\int_{E_{1}}^{E_{2}}A(E,\vec{\xi})~{}dE$,
where $\vec{\xi}$ parameterizes the effective area and includes quantities
such as the geometric area, filter thicknesses, and chemical compositions that
are initially estimated during ground calibration. Uncertainties in
$\vec{\xi}$ are quantified via the prior distribution, $p(\vec{\xi})$;
examples generated and used for ACIS analyses can be found in Drake et al.
(2006); Kashyap et al. (2008); Lee et al. (2011) and Xu et al. (2014). We
suppress the subscript $i$ throughout this section for notational simplicity
because we are restricting consideration to an arbitrary instrument.

Following the discussion in Section 2.2, we specify the prior distribution on
(the logarithm of) the effective areas of $U$ energy bands,
$\\{{B}(\mathcal{E}_{u},\vec{\xi})\\}$ for ${u=1,\ldots,U}$, as a multivariate
Gaussian distribution with expected values, $\beta_{u}$, and variances,
$\Lambda_{u}$, for each energy band $u$, and with correlations, $\rho_{uv}$,
between all pairs of distinct energy bands, $u\neq v$. Among the $b_{u}$, the
$\Lambda_{u}$, and the $\rho_{uv}$, only the correlations, $\rho_{uv}$, (or,
more precisely, the Monte Carlo estimates of the $\rho_{uv}$) are used in our
data analyses. To distinguish the Monte Carlo estimates of the prior means and
variances of the $B_{i}$ obtained here from those we actually use, we
introduce new notation for these quantities that differ from those used in
Paper I and in Section 2.2. The current methods for computing $\Lambda_{u}$
are still somewhat experimental, so we rely instead on the $\tau_{i}$ values
from IACHEC scientists.

For a given instrument and energy band $u$, the expectation is

$\beta_{u}=\int{B}(\mathcal{E}_{u},\vec{\xi})~{}p(\vec{\xi})~{}d\vec{\xi}.$
(18)

In practice, calibration scientists set $a_{i}$ to be the prior estimate of
$A_{i}$ (on the original scale) based on their best information and
experience. Transforming to a logarithmic scale, we might set our prior
estimate of $B_{i}$ to $\log a_{i}$, which is the choice of prior mean we
suggest in Section 2.2. Eq. (18) can be used in the absence of such intuition
or if we prefer to use a parameterized model for ${B}$. (This is particularly
relevant for the correlations, since we have less intuition for them.) In this
case, a reasonable strategy is to proceed via Monte Carlo integration of
Equation (18). More precisely, we obtain a calibration sample,
$\\{{B}^{(k)}(\mathcal{E}_{u}),k=1,\ldots,K\\}$, that quantifies prior
uncertainty in $B(\mathcal{E}_{u},\vec{\xi})$, for example by obtaining a
sample of size $K$ from $p(\vec{\xi})$ and computing
$B(\mathcal{E}_{u},\vec{\xi}^{(k)})$ for each sample, $\vec{\xi}^{(k)}$, from
$p(\vec{\xi})$. The $B^{(k)}(\mathcal{E}_{u})$ is expressed as
$\log\left[\int_{E_{1}}^{E_{2}}A^{(k)}(E)~{}dE\right]$ for each $u$ (given all
Monte Carlo samples $A^{(k)}(E)$ for $E\in\mathcal{E}_{u}$). The Monte Carlo
version of Eq. (18) is then

$\hat{\beta}_{u}=\frac{1}{K}\sum_{k=1}^{K}{B}^{(k)}(\mathcal{E}_{u}).$ (19)

The prior variance of ${B}(\mathcal{E}_{u},\vec{\xi})$ and its Monte Carlo
estimate are

$\displaystyle\Lambda_{u}=\int[B(\mathcal{E}_{u},\vec{\xi})-\beta_{u}]^{2}~{}p(\vec{\xi})~{}d\vec{\xi}\quad\hbox{
and
}\quad\hat{\Lambda}_{u}=\frac{1}{K-1}\sum_{k=1}^{K}\left[{B}^{(k)}(\mathcal{E}_{u})-\hat{\beta}_{u}\right]^{2},$
(20)

respectively. Finally, the prior correlation between
${B}(\mathcal{E}_{u},\vec{\xi})$ and ${B}(\mathcal{E}_{v},\vec{\xi})$ is the
covariance normalized by the respective standard deviations,

$\rho_{uv}=\frac{1}{\sqrt{\Lambda_{u}\Lambda_{v}}}\int\left[B(\mathcal{E}_{u},\vec{\xi})-\beta_{u}\right]\left[\tilde{B}(\mathcal{E}_{v},\vec{\xi})-\beta_{v}\right]~{}p(\vec{\xi})~{}d\vec{\xi}$
(21)

with Monte Carlo estimate,

$\hat{\rho}_{uv}=\frac{1}{(K-1)\sqrt{\hat{\Lambda}_{u}\hat{\Lambda}_{v}}}\ \
\sum_{k=1}^{K}\left[{B}^{(k)}(\mathcal{E}_{u})-\hat{\beta}_{u}\right]\left[{B}^{(k)}(\mathcal{E}_{v})-\hat{\beta}_{v}\right].$
(22)

###  2.4 Practical Implementation

This section discusses practical implementation of our methods, specifically,
in terms of normalization of observed counts/fluxes and computation of
correlation matrices.

####  2.4.1 Normalization in Practice

Because X-ray data analysis packages such as xspec return the $f_{ij}$ and
their uncertainties, it is convenient to rewrite Eq. (15) in terms of

$\tilde{y}_{ij}=\log\frac{f_{ij}}{\tilde{f}}=\log\frac{c_{ij}}{T_{ij}a_{i}\tilde{f}},\quad\tilde{B_{i}}=\log\frac{A_{i}}{a_{i}},\quad\hbox{and}\quad\tilde{G_{j}}=\log\frac{F_{j}}{\tilde{f}}$
(23)

to obtain

$\tilde{y}_{ij}=\tilde{B_{i}}+\tilde{G_{j}}-\frac{1}{2\sigma_{i}^{2}}+e_{ij},$
(24)

where $\tilde{f}$ is a fiducial flux (usually the maximum of the $f_{ij}$)
used to normalize the data to the range $[0,1]$. Model (24) is functionally
equivalent to Model (15). This definition of $\tilde{G}_{j}$, normalized by
$\tilde{f}$, which depends on data, is only introduced for computational
convenience and does not affect the model or its interpretation. Technically,
using a data-dependent “parameter”, here $\tilde{G}_{j}$, implies a data-
dependent prior distribution, which is generally not legitimate from a
Bayesian viewpoint. However, because using a flat prior on $G_{j}$ is the same
as using a flat prior on $\tilde{G}_{j}=G_{j}-\log{\tilde{f}}$, there is no
actual effect in our implementation. Thus, Model (24) has the same form as
Model (15), so we can embed it into the Bayesian hierarchical model described
in Paper I to obtain the full posterior distribution, estimates, and error
bars for the $\tilde{B}_{i}$ and $\tilde{G}_{j}$. Paper I suggests setting
$b_{i}=0$ (because $\tilde{B}_{i}=0$ implies that $A_{i}=a_{i}$).

####  2.4.2 Deriving Correlations in Practice

We proceed by computing numerous instances of instrument effective areas that
are varied in controlled ways dictated by current knowledge of uncertainties
in calibration. The basis of the method is to generate a so-called calibration
sample of areas that represents the range of uncertainties of the effective
area, including all the correlations between different energies. The approach
to generating the calibration samples for the different instruments we study
here is common to all instruments, with some additional complexity built into
the Chandra samples. We describe the method in brief below and refer the
reader to Drake et al. (2006) for more complete description.

We devise a “perturbation function” that comprises piecewise cubic segments
that stretch between the natural absorption edges of the different materials
encountered along the optical path for a given instrument. This function
varies about unity by random amounts but is constrained within fixed limits
based on specified calibration uncertainties by the cubic function whose
parameters are randomly drawn from a truncated Gaussian distribution. The
perturbation function is applied as a multiplicative factor to the different
subassembly component contributions to the effective area, which are
considered on a case-by-case basis. For instance, in the case of Chandra/ACIS,
six plausible mirror effective areas are used, uncertainties in the optical
blocking filter and contamination and contamination transmittance are modeled
by altering the optical depth of each chemical component within their known
uncertainties and recomputing ensemble transmittance, and CCD quantum
efficiencies are computed for different realizations of depletion depth and
SiO2 layers. Details of this process are given in (Drake et al., 2006, see
also Drake et al. 2021, in preparation).

Chandra/HETGS and Chandra/LETGS use a similar approach with additional
perturbation functions applied for the transmission grating diffraction
efficiences. In the case of the other instruments considered here, the
perturbation function approach alone is used.

##  3 Observations and Data Processing

Three data sets are considered in Paper I and we add a fourth in this paper.
Here, we detail how the data sets are handled and the required data
processing.

###  3.1 Supernova Remnant 1E0102.2-7219

As in Paper I, the fluxes of the emission line complexes of O and Ne in the
X-ray spectra of SNR 1E0102.2$-$7219 are taken from the detailed comparison of
13 instruments by Plucinsky et al. (2017). Briefly, the spectra of each non-
dispersive instrument are fit with a model with five free parameters: an
overall normalization and four emission line fluxes of O vii, O viii, Ne ix
and Ne x. Because the emission lines of the same element have similar
energies, their effective areas are comparable and highly correlated, so we
combined the O and Ne line fluxes to create two fluxes for each instrument. In
our statistical analysis, these two fluxes are treated as “sources”: one for O
and another for Ne. The data are normalized to the O or Ne fluxes obtained by
the XMM-Newton pn instrument. By requiring that each instrument analysis uses
the same model, except primarily for the strengths of the emission lines, the
$q(E;\theta_{j})$ values do not depend on the instrument, nor on the measured
line fluxes.

Table 3 shows the $\tau_{i}$ values assigned to the effective areas for each
instrument considered. The values were taken from Table 1 using the bandpass
that covers the lines of interest: 0.54-0.80 keV for O and 0.8-1.2 keV for Ne.
For the correlation matrix, there is only one off-diagonal term, which we set
to 0.88 for ACIS instruments and 0.82 for XMM instruments, as derived as in §
2.3.2.

Table 3: Heterogeneous $\tau_{i}$ Values for 1E0102 AnalysisaaValues for $\tau_{i}$ are in percentages for each combination of instrument and line complex, using $\tau_{i}$ values taken from Table 1.  Instrument | Oxygen | Neon
---|---|---
XMM/RGS1 | 5 | 5
XMM/MOS1 | 6 | 6
XMM/MOS2 | 6 | 6
XMM/pn | 2 | 2
ACIS-S3 | 3 | 3
ACIS-I3 | 3 | 3
ACIS/HETG | 3 | 3
Suzaku/XIS0 | 15 | 10
Suzaku/XIS1 | 15 | 10
Suzaku/XIS2 | 15 | 10
Suzaku/XIS3 | 15 | 10
Swift/XRT-WT | 10 | 7.5
Swift/XRT-PC | 10 | 7.5

###  3.2 Sources from the 2XMM Catalog

We select a sample of X-ray sources from the 2nd European Photon Imaging
Camera (EPIC) Serendipitous Source (2XMM) Catalog (Watson et al., 2009). EPIC
consists of three X-ray cameras with CCD sensors mounted on the ESA spacecraft
XMM-Newton (Jansen et al., 2001); the EPIC-pn (Strüder et al., 2001) and two
EPIC-MOS (Turner et al., 2001) cameras observe celestial sources quasi-
simultaneously within their co-aligned fields of view. For this analysis,
v14.0 of the Science Analysis System (SAS; Gabriel et al., 2004) is used, as
well as the calibration files as available in 2016. A description of the data
reduction and spectral extraction procedure appears in Read et al. (2014).
Soft, medium, and hard bands are defined to be the 0.5-1.5, 1.5-2.5, and
2.5-10 keV bands, respectively. Due to variability,different observations of
the same source are treated as separate sources for a total of 41 observations
of 35 distinct sources. The normalizing (maximum) fluxes in the soft, medium,
and hard bands are 0.138, 0.000701, and 0.00223 photons cm-2 s-1, the
brightest sources in their respective lists. The data are provided in tables
A-A in the Appendix. There are two primary features of the analysis procedure
that make the results suitable for Concordance analysis: 1) the count spectra
for the different instruments (pn, MOS1, MOS2) are fit simultaneously to power
laws, so that the spectral slopes, $\theta_{j}=\Gamma_{j}$ (where
$n_{E}[\Theta_{j}]=n_{j}E^{-\Gamma_{j}}$) do not depend on the instrument
combination, and 2) the sources are faint enough that pileup is not an issue.

Table 4 gives the $\tau_{i}$ values assigned to the effective areas for the pn
and MOS instruments and the three bandpasses. The values are taken from Table
1 using the 0.54-0.8 and 0.8-1.2 keV $\tau_{i}$ values for the soft band, the
1.2-2.2 keV $\tau_{i}$ value for the medium band, and an average of the 2.2-10
keV $\tau_{i}$ values for the hard band. Table 5 gives the correlation matrix
values, $\rho_{mn}$, computed for the pn and MOS instruments and the three
bandpasses used in the 2XMM catalog while Table 6 provides these values for
the bandpasses used in the XCAL analysis. For simplicity, we assume that the
$\rho_{mn}$ are the same for each instrument.

Table 4: Heterogeneous $\tau$ Values for 2XMM and XCAL AnalysesaaValues for $\tau$ are percentages for each combination of instrument and line complex, using $\tau$ values from Table 1. Instrument | Soft band | Medium band | Hard band
---|---|---|---
pn | 2 | 2 | 2.3
MOS1 | 6 | 6 | 7.3
MOS2 | 6 | 6 | 7.3
Table 5: Correlation matrix for 2XMM Analyses  Band | Soft band | Medium band | Hard band
---|---|---|---
Soft band | 1 | 0.61 | 0.13
Medium band | 0.61 | 1 | 0.53
Hard band | 0.13 | 0.53 | 1
Table 6: Correlation matrix for XCAL Analyses  Band | Soft band | Medium band | Hard band
---|---|---|---
Soft band | 1 | 0.63 | 0.20
Medium band | 0.63 | 1 | 0.52
Hard band | 0.20 | 0.52 | 1

###  3.3 Active Galaxies from the XCAL Sample

Another set of EPIC spectra used to validate the model is the so-called “XMM-
Newton Cross-Calibration” (XCAL) sample.222Details of the XCAL processing are
available in Section 4 of the XMM calibration memo XMM-SOC-CAL-TN-0052, by
Stuhlinger et al. (2010). The memo is available at
https://xmmweb.esac.esa.int/docs/documents/CAL-TN-0052.ps.gz. This is a sample
of radio-loud Active Galactic Nuclei (AGN), primarily blazars, observed
routinely by XMM-Newton in the framework of its in-flight calibration program
(Guainazzi et al., 2015). As with the 2XMM sample described in §3.2, the
sources are variable. In this case, there are more blazars and high signal
observations, for a total of 108, 103, and 94 observations of 22 distinct
sources in the soft, medium, and hard bands that exceeded a flux limit
criterion without highly discrepant fluxes between the three instruments. The
normalizing (maximum) fluxes in the soft, medium, and hard bands are set to
0.126, 0.0156, and 0.0154 photons cm-2 s-1, the brightest sources in their
respective lists. Data are provided in tables A-A) in the Appendix. As with
the 2XMM sources, the pn, MOS1, and MOS2 data were fit simultaneously to power
law spectra so that $\Gamma_{j}$ is the same for each instrument. However,
compared to the 2XMM sources, the XCAL sources are bright, often exceeding the
count rate threshold beyond which the fraction of events affected by pile-up
is no longer negligible (Jethwa et al., 2015). Spatial regions on the detector
affected by pile-up are removed by excising the core of the telescope Point
Spread Function (PSF) up to an observation-dependent radius. This radius is
determined on the basis of the ratio between non X-ray diagonal and standard
X-ray “patterns” (measure of the event shape in the CCD) in EPIC-MOS (Jethwa
et al., 2015); and by visual inspection of the pattern distribution curves in
EPIC-pn using the SAS task epatplot in EPIC-pn. This “PSF core excising
method”, while unavoidable to retain the highest possible fidelity of event
spectral calibration, may introduce excess variance in the $f_{ij}$ via
systematic uncertainty in the energy-dependent correction for the fraction of
events scattered into or out of the annular spectral extraction region by the
PSF (the so-called “Encircled Energy Correction” fraction). Values for
$\tau_{i}$ and $\rho_{mn}$ are the same as for the 2XMM sample.

###  3.4 Active Binary Capella

Capella ($\alpha$ Aur AB; G1III+G8III; 13pc) is a spectroscopic binary that is
the brightest line-dominated source accessible to non-Solar X-ray missions. It
is remarkably steady for a coronal source, having never exhibited significant
flaring. While it does vary over timescales of months, it does not show any
evidence of flux variability over timescales of weeks or less. Consequently,
it has often been used as a calibration target, in particular with Chandra. It
has been observed several times with different detector and grating
combinations in close proximity (see Table A in the Appendix). These
observations allow us to carry out an assessment of the internal cross-
calibration of the Chandra grating spectrometers.

We estimated the total fluxes in each of several strong lines: the highly-
ionized lines of Fe xvii (at 15Å and 17Å) whose formation temperatures overlap
the peak emission measures of Capella, and the hydrogenic lines of Ne x
(12.13Å) and O viii (18.96Å). For the purposes of this calculation, we treat
each of the four emission lines as different sources. We then form 21 epoch
groups, comprised of observations that are within 0.1 yr of each
other333Aug/Sep99, Mar00, Feb01, Apr02, Oct02, Sep03, Sep04, Mar05, Oct05,
Apr06, Apr07, Apr08, Apr09, Nov09, Nov/Dec10, Dec11, Dec13, Dec14, Jul16,
Sep16, and Dec18; see Table A., giving a total of 84 sources. Similarly, $+1$
and $-1$ grating orders are treated distinctly for each of four
grating/detector combinations, ACIS-S/HEG, ACIS-S/MEG, ACIS-S/LEG, and
HRC-S/LEG, for a total of eight instruments. The fluxes were normalized to the
maximum values for each emission line: 28.02, 60.23, 49.66, and 23.33 $\times$
$10^{-13}$ erg s-1 cm-2 for O viii, Fe xvii 17Å, Fe xvii 15 Å, and Ne x,
respectively. We use CIAOv4.11 to extract the dispersed spectra, and compute
the effective areas using the contamination corrections as in CALDBv4.8.0.1.

The values of $\tau_{i}$ are taken from Table 8 and the correlation matrix is
given in Table 9.

##  4 Results

Here we present results from new measurements and extensions to the results in
Paper I. In each case, we generate 10,000 Monte Carlo replicates from the
respective posterior distributions as the basis for our statistical
inferences.

###  4.1 1E0102

These data provide a illustration of a case where there are many instruments
that obtain data on the same source, shown in Figs. 2 and 3. The effect of
allowing heterogeneous $\tau$ values is apparent in both cases. Generally,
when the prior distribution on an instrument’s effective area is more
uncertain than average, giving a relatively large value of $\tau_{i}$, then
the data for that instrument is given less weight, so the posterior estimate
of its effective area is more likely to deviate from the prior estimate by
comparison to when all instruments have equally uncertain prior estimates. In
addition, the posterior range of the deviation is more likely to be large when
$\tau_{i}$ is larger. The Suzaku results for the O lines all show this effect.
When effective areas are correlated between the O and Ne data sets, the
ACIS-I3 point is particularly affected due to the discrepant results obtained
when Ne and and O data are considered independently.

Figure 2: Results of the Concordance analysis for the data from the SNR 1E0102
for the combination of fluxes of the lines of O vii and O viii. The
$\tau=0.025$ (in black) and $0.05$ (in yellow) results are the same as given
in Paper I and are shown to elucidate the effects of including heterogeneous
$\tau$ values (in blue) and adding effective area correlations (in red). The
error bars represent the 90% (5%-95%) confidence regions on the posterior
estimate of $A_{i}/a_{i}$, as defined in §2.2. When effective areas are
correlated between the O and Ne data sets, the ACIS-I3 point is particularly
affected due to the discrepant results obtained when Ne and and O data are
considered independently. Figure 3: Same as Fig. 2 except for the combination
of fluxes of the lines of Ne ix and Ne x (see text). Unlike the case for the O
lines, the posterior estimates of the effective area tend to be more stable in
this band and relatively independent of the uncertainties in the effective
area priors.

###  4.2 XMM Samples

Figures 4 and 5 show the results of the Concordance analysis for the two XMM-
Newton data. In this example, there are many sources and few instruments, in
contrast with the 1E0102 data set. The 2XMM results show a high degree of
consistency between the instruments, consistently favoring 3-5% increases to
the effective areas of the MOS detectors across all bands and a corresponding,
slight decrease to the pn effective area. With the use of individualized
$\tau$ values, the Concordance analysis drives the pn effective areas toward
the prior, as one might expect due to the significantly smaller $\tau$
assigned to the pn compared to the MOS detectors.

The XCAL sample, shows similar trends to that of the 2XMM sample with with a
more significant indication that the MOS2 detector’s effective area should be
increased 2-3% more than that of the MOS1 detector.

Figure 4: Concordance results for the 2XMM sample. Results are color-coded as
in Fig. 2. When the $\tau$ values are allowed to vary by instrument, the
“heterogeneous” case, the posterior for the pn centers on the prior, due to
the smaller value of $\tau$ than used for the MOS detectors. At the same time,
higher effective areas for the MOS detectors are indicated across all bands.
Figure 5: Same as Fig. 4 except for the XCAL sample. These results are
generally consistent with those from the 2XMM sample but with a somewhat
stronger indication that the MOS2 effective area should be increased relative
to MOS1

###  4.3 Capella Line Fluxes for Chandra Grating Spectrometers

Results from the Concordance analysis as applied to the Capella data are shown
in Fig. 6. There are several features of interest. First, the effective area
corrections for the LETGS (ALEG and HLEG) are generally negative while those
of the HETGS (HEG and MEG) are generally positive. These corrections are
consistent with preliminary results on independent data where the instruments
are cross-calibrated with alternating observations of Mk 421. Second, the $+1$
and $-1$ orders generally agree well for all instruments and wavelengths.
Third, when the effective area correlations are included, the posterior
effective areas for the longer wavelengths (O viii and Fe xvii $\lambda 17$)
more strongly deviate from their priors.

Figure 6: Concordance analysis using measurements of four emission lines using
the the Chandra grating spectrometers. Emission lines used in the various
panels are: a) Ne x, b) Fe xvii $\lambda 15$, c) Fe xvii $\lambda 17$, d) O
viii. Color coding of results are as in Fig. 2. Features to note are 1) that
the effective area corrections for the LETGS (ALEG and HLEG) are generally
negative while those of the HETGS (HEG and MEG) are generally positive and 2)
the $+1$ and $-1$ orders generally agree well.

###  4.4 Method Validation and Assessment

Paper I includes a series of numerical studies that explore the statistical
properties of our method. For example, Figure 2 of Paper I illustrates that
our posterior distributions of the effective areas cover the true effective
areas in a simulation study. Figure 7 of Paper I then goes on to contrast the
estimated 95% intervals for log-fluxes constructed using the standard
instrument-specific estimates with the combined estimate based on our
posterior distribution, illustrating how our Bayesian process achieves a
single consistent estimate for each flux but with smaller errors than the
standard estimates. Here we consider additional ways to evaluate how robust
our method’s results may be, supplementing the simulations performed in Paper
I.

####  4.4.1 Simulation Studies

The method developed and applied in Paper I produces Bayesian posterior
distributions for each estimated quantity. The main quantities of interest
here are the fractional corrections to instrument effective areas, given by
$\log A_{i}-\log a_{i}$, and the fractional corrections to the estimated
fluxes of sources, given by $\log F_{j}-\log f_{ij}$.

We demonstrate how the Concordance method yields accurate and reliable
estimates of $A_{i}$ and $F_{j}$ with a simulation study. The simulation
involves 40 simulated sources observed by each of five instruments. We set the
prior means of the effective areas of the instruments to differ from their
actual effective areas by $\log A_{i}-\log a_{i}=[0,1,-1,2,-2]$, for $i=1..5$,
respectively. For example, for instrument 4, the true effective area is
systematically higher than the prior mean by a factor of $e^{2}$, resulting in
flux estimates that are systematically too high compared to the true values.
The $\tau_{i}$ values are all 1.0 in this simulation, indicating large
uncertainties in prior estimates of the effective areas, and the measurement
uncertainties (i.e., $\sigma_{ij}$) are all set to 0.5 on the log scale,
similarly indicating large uncertainties, except for instrument 5, for which
$\sigma=0.1$. This simulation setup is designed to test the robustness of the
Concordance method to data from an instrument with high signal/noise but a
systematically biased effective area. We replicated this entire set up 200
times and processed each replicate with our Concordance method. A
representative replicate is shown in Fig. 7. The simulations demonstrate that
the Concordance method provides source flux estimates that are substantially
better than would be obtained by simply using the prior means as estimates of
the effective areas. The replicate simulations indicate that the 95% equal-
tailed posterior intervals cover the true values of the effective areas and
source fluxes over 99% of the time. Thus, we not only obtain better estimates
of the source fluxes, but also estimate the effective area corrections well.
Note that the instrument-specific flux estimates can deviate substantially
from the true values for any given source, so that when using a set of such
estimates, the weighted average ratio estimators for the flux (i.e.,
$(\sum_{j}f_{ij}/\sigma_{j}^{2})/\sum_{j}1/\sigma_{j}^{2}$) would be generally
biased toward the instrument with the highest signal/noise, regardless of the
accuracy of the instrument’s effective area.

Figure 7:  Posterior distributions of the logarithms of the effective area
corrections (top row) and on the logarithms of the source fluxes (bottom row,
for five sources) for one of 200 replicate simulations, each involving five
instruments and 40 sources. In the simulation, all sources have a true value
of $\log F_{j}=1$ but the prior means of the instrument effective areas are
off by factors of $\exp([0,1,-1,2,-2])$ for $i=1,\ldots,5$. The top row shows
that the Concordance method generates posterior distributions of the effective
area corrections that are well centered on the true values for each instrument
(shown as vertical solid lines in various colors). The bottom row shows that
the posterior distributions of the source fluxes are well centered on the true
values (vertical solid black lines), even while the instrument-specific
estimates based on the prior means of the effective areas (vertical dotted
lines of colors corresponding to those of the instruments in the top row) can
be individually erroneous by large factors.

With an additional pair of simulations, we also quantified the improvement
that can be obtained with the Concordance method. In the first case, we
simulated $M=3$ instruments with $A_{i}/a_{i}=[1,1,0.9]$ for $i=[1,2,3]$,
$\tau_{i}=0.05$, and $N=20$ sources with true fluxes all equal to 1. There
were 200 independent simulations and analyses for each setup. The sources were
assumed to have good signal/noise as might be expected for calibration
observations: $\sigma_{i}=0.03$ (i.e., the $c_{ij}$ were drawn from a Poisson
distribution with a mean of about 1100). The second case is the same as the
first except there is a higher statistical precision for observations with
instrument 3: $\sigma_{3}=0.003$. Source fluxes were estimated for each
simulation using the Concordance method and also using the above-mentioned
weighted average of ratio estimators:
$(\sum_{j}f_{ij}/\sigma_{j}^{2})/\sum_{j}1/\sigma_{j}^{2}$. The 95%
uncertainty bounds on the flux estimates and the coverage fraction where the
true flux is included within the uncertainty bounds are shown in Table 7. The
Concordance estimator generated uncertainty intervals that were accurate and
covered the ground truth, in contrast to the ratio estimator, which, despite
the widths of the confidence intervals being nominally smaller, generated
biased estimates and did not cover the true fluxes. The situation was worse
for the second case, where the Concordance intervals did not change
appreciably, but the ratio estimators were biased low by $\approx$10%,
reflecting the higher statistical weight given to an instrument with a biased
estimate of its effective area. Indeed, in this latter case, there was only
one source (of 20) in only one simulation (of 200) where the ratio estimator
confidence interval included the true value. This pair of simulation setups
illustrates the robustness of the Concordance method to erroneous effective
area priors, showing that it is appropriate to use in calibration work where
robustness and accuracy are highly valued.

Table 7: Results from Two Concordance Simulations  Simulation SetupaaSetups 1 and 2 are the same except that instrument 3 (of 3) has 3% statistical errors for setup 1 and 0.3% statistical errors for setup 2. The prior for the effective area of instrument 3 is 10% higher than its true value in both setups. See Section 4.4.1 for details. | Flux Estimation | 95% Flux RangebbAverage 95% confidence intervals for source flux estimates; the true fluxes for all sources are set to 1. | $r_{95}$ccFraction of flux estimates covering true fluxes at the 95% confidence level out of 200 simulations of 20 sources each.
---|---|---|---
| Method | $F_{\rm lo}$ | $F_{\rm hi}$ |
1 | Concordance | 0.903 | 1.033 | 0.964
1 | Ratio Estimator | 0.927 | 0.994 | 0.372
2 | Concordance | 0.905 | 1.031 | 0.990
2 | Ratio Estimator | 0.896 | 0.907 | 0.000

####  4.4.2 Posterior Histograms

We have found that the posterior distributions of the effective area
corrections are typically well described by Gaussians, as shown in Fig. 8,
parts A-C. These three examples were randomly chosen among the dozens of such
histograms generated in our analysis of the data from §4.1-4.3. Occasionally,
however, there are histograms that are not obviously Gaussian, so we also show
three “bad” examples. In one case, Fig. 8, part D, there is a distinct “notch”
in a side of the distribution and in two cases (Fig. 8, parts E and F), there
is noticeable skew – tails to large fractional corrections. These three cases
were quite rare but give warning that there may be inconsistencies in the
underlying data. One known source of error that is not accounted for in our
analysis is in the shape of the response function, $\Phi_{k}(E)$. For
instruments like ACIS and the EPIC detectors, the low energy response is
somewhat uncertain and difficult to calibrate.

Figure 8: Posterior histograms for the fractional variation of the effective
area. The data sets are (A) 2XMM data, hard band, XMM/pn, correlated $\tau$
values; (B) 2XMM data, hard band, XMM/MOS2, heterogeneous $\tau$ values; (C)
XCAL data, medium band, XMM/pn, correlated $\tau$ values; (D) XCAL data, soft
band, XMM/MOS1, $\tau=0.05$ for all instruments; (E) 1E0102 data, O lines,
Chandra/ACIS-S3, $\tau=0.025$ for all instruments; and (F) 1E0102 data, O
lines, Swift XRT/PC, $\tau=0.05$ for all instruments. Histograms A-C are
typical, chosen randomly from several dozen; the distributions are well
approximated as Gaussians. Histograms D-F are atypical, showing skew or other
non-Gaussian shapes.

####  4.4.3 Sensitivity to Uncertainties in Priors

The specification of priors is typically under scrutiny for Bayesian analysis
in practice. Typically researchers conduct sensitivity analysis to study the
outcome sensitivity with respect to small perturbations of the priors. In our
setting, the sensitivity of the results with respect to $\tau$ values are
revealed by the comparison between heterogeneous $\tau$ values versus the two
homogeneous $\tau$ value choices. However, for the correlation matrix in the
prior distribution, we adopted Monte Carlo estimates, which is subject to
random variations. Thus, we undertook a simple example of the test of
sensitivity in this paper, but it is feasible for any user of the concordance
tools. Namely, for the Capella data, instead of adopting the full correlation
matrix as given in Tables 9 and 10, we only keep the correlations between the
two Fe bands. Again, we can test out different variations of the correlation
matrix with the same procedure and similar analysis. Thorough sensitivity
analysis requires extensive testing on a carefully designed set of variations
of the prior distributions. See Fig. 9 for the results of applying this
variation to the Capella data. By comparing the results of Fig. 9 with the
original, Fig. 6, we can reveal the sensitivity of the results as opposed to
variations of correlations between Fe bands and others, and the correlations
of others (Ne and O) within their own. We can see that in Fig. 9, the O and Ne
shows nicely aligned results under correlated effective areas versus
uncorrelated effective areas. But this is not true for Fig. 6, where O and Ne
are still correlated with each other and with other bands, especially for
HEG+. Furthermore, while the adjustments for the two Fe chanels are similar
across the two figures, the adjustments for Ne are very different across the
two figures. The resulting adjustments of effective areas not only deviate
more significantly from zero but also have smaller error bars when
correlations are taken into account. This makes intuitive sense because the
benefit of accounting for correlations among effective areas is to obtain
sharper or more informative estimates.

Figure 9: Same as Fig. 6 except that only the correlations between the two Fe
bands are kept while others are set to be zero, for both ACIS and HRC
instruments.

##  5 Conclusion and Directions for Development

These data sets provided an excellent foundation for the Concordance project,
whose goal is to determine quantitative and objective evidence for making
effective area adjustments in order to improve agreement between instrument
measurements. The process applied here is available for use in studies such as
we have undertaken. There are some avenues to explore for expanding this
particular implementation of the Concordance analysis.

###  5.1 Correlations between source bandpass fluxes

Several types of calibration sources have simple spectra, which is why they
are often used in cross-calibration. Examples are isolated neutron stars with
blackbody spectra, blazars with power law spectra, and supernova remnants and
clusters of galaxies with thermal spectra. To the extent that these spectra
can be characterized by only a few parameters, such as a power law slope, then
the flux in one band is closely related to that in an adjacent band.
Furthermore, many types of source have smoothly continuous spectra – their
spectral fluxes are tightly correlated on small scales. Modeling a many
bandpasses of a blazar spectrum with a series of power laws with different
slopes would lead to unphysical discontinuities at bandpass boundaries. Thus,
it would be advantageous to take advantage of this astrophysical knowledge and
include spectral band correlations due to spectral continuity and simplicity.
The XMM-Newton-Chandra blazar XCAL sample is an excellent data set to examine
next, involving three Chandra configurations and all four XMM-Newton X-ray
detectors and covering the energy range from 0.1 to 10 keV using simultaneous
observations of active galaxies obtained over 20 years of operation.
Preliminary results have been reported at various IACHEC meetings.

###  5.2 Secular variations of instrument responses

While many sources may well vary erratically, instrument behavior can often be
subject to gradual degradation. With adequate modeling of many observations,
one avenue to explore is how to link instrument effective areas over time
within the Concordance framework.

###  5.3 Nonlocal instrument responses

There are definite difficulties that are encountered when the detector energy
response $\Phi_{k}(E)$ has a non-Gaussian component, a broad asymmetry, or
bimodality because systematic errors in the response function can appear in an
apparently unrelated bandpass. Response function errors may be responsible for
some of the artifacts in our posterior histograms (see §4.4.1). The response
functions of solid state detectors have escape peaks that can generate events
at a significantly different apparent energy than that of the incoming photon.
Modeling the effects of systematic errors in response functions is possible in
principle, especially with methods such as used to determine the effective
correlation function (see § 2.3.2). One approach for dealing with this issue
would be to expand the Concordance mathematical model to include a term to
account for variance of the $T_{ij}$ values.

Support for this work was provided in part by the National Aeronautics and
Space Administration (NASA) through the Smithsonian Astrophysical Observatory
(SAO) contract SV3-73016 to MIT for support of the Chandra X-Ray Center (CXC),
which is operated by SAO for and on behalf of NASA under contract NAS8-03060.
Support was also provided by NASA under contract NAS8-39073 to SAO. This work
was conducted in collaboration with the CHASC International Astrostatistics
Center. CHASC is supported by NSF DMS-18-11308, DMS-18-11083, and
DMS-18-11661. DvD was supported in part by a Marie-Skodowska-Curie RISE
(H2020-MSCA-RISE-2015-691164, H2020-MSCA-RISE-2019-873089) Grants provided by
the European Commission. .

Table 8: Values of $\tau_{i}$ for Capella AnalysesaaValues for $\tau$ are in percentages. Line | Ne x $\lambda$12 | Fe xvii $\lambda$15 | Fe xvii $\lambda$17 | O viii $\lambda$19
---|---|---|---|---
HEG | 4 | 5 | 5 | 10
MEG | 4 | 5 | 5 | 10
LEG/A | 4 | 5 | 5 | 10
LEG/H | 7 | 7 | 7 | 7
Table 9: ACIS Correlation matrix used for Capella  Line | Ne x $\lambda$12 | Fe xvii $\lambda$15 | Fe xvii $\lambda$17 | O viii $\lambda$19
---|---|---|---|---
Ne x $\lambda$12 | 1 | 0.96 | 0.92 | 0.89
Fe xvii $\lambda$15 | 0.96 | 1 | 0.99 | 0.97
Fe xvii $\lambda$17 | 0.92 | 0.99 | 1 | 0.99
O viii $\lambda$19 | 0.89 | 0.97 | 0.99 | 1
Table 10: HRC Correlation matrix used for Capella  Line | Ne x $\lambda$12 | Fe xvii $\lambda$15 | Fe xvii $\lambda$17 | O viii $\lambda$19
---|---|---|---|---
Ne x $\lambda$12 | 1 | 0.84 | 0.71 | 0.62
Fe xvii $\lambda$15 | 0.84 | 1 | 0.83 | 0.74
Fe xvii $\lambda$17 | 0.71 | 0.83 | 1 | 0.91
O viii $\lambda$19 | 0.62 | 0.74 | 0.91 | 1

##  Appendix A Measured fluxes

Here we list all the fluxes used in the calculations described above: for
sources in the 2XMM catalog (see Section 3.2) in the soft (Table A), medium
(Table A), and hard (Table A) bands; for active galaxies from the XCAL sample
(Section 3.3) in the soft (Table A), medium (Table A), and hard (Table A)
bands; and the line fluxes measured during the various Capella grating
observations with Chandra (Section 3.4; Table A). Note that the data used for
the analysis of SNR 1E0102.2-7219 (Section 3.1) are given in Chen et al.
(2019).

Table A.1: 2XMM Concordance Fluxes – Soft BandaaFluxes are normalized to 0.138
photons cm-2 s-1.

Target | pn | MOS1 | MOS2
---|---|---|---
| $f_{ij}$ | $\sigma_{ij}$ | $f_{ij}$ | $\sigma_{ij}$ | $f_{ij}$ | $\sigma_{ij}$
1127$-$145 | 0.0075 | 0.0009 | 0.0074 | 0.0009 | 0.0077 | 0.0010
1E0919$+$515 | 0.0068 | 0.0013 | 0.0081 | 0.0015 | 0.0083 | 0.0015
4C06.41 | 0.0024 | 0.0001 | 0.0026 | 0.0001 | 0.0026 | 0.0001
APM 08279$+$5255 | 0.0035 | 0.0006 | 0.0035 | 0.0006 | 0.0036 | 0.0006
CAL 87 | 0.0125 | 0.0011 | 0.0136 | 0.0012 | 0.0134 | 0.0012
Cen X-4 | 0.0180 | 0.0027 | 0.0180 | 0.0027 | 0.0190 | 0.0028
CoD-33 7795 | 0.1523 | 0.0222 | 0.1578 | 0.0230 | 0.1618 | 0.0236
ESO323$-$G077 | 0.0011 | 0.0004 | 0.0012 | 0.0005 | 0.0012 | 0.0005
GRB080411 | 0.0113 | 0.0009 | 0.0128 | 0.0010 | 0.0128 | 0.0010
Holmberg IX | 0.0152 | 0.0026 | 0.0154 | 0.0027 | 0.0154 | 0.0027
IRAS 13197-1627 | 0.1018 | 0.0402 | 0.1066 | 0.0419 | 0.1062 | 0.0418
LBQS 1228+1116 | 0.0046 | 0.0001 | 0.0048 | 0.0001 | 0.0048 | 0.0001
M31 NN1 | 0.0032 | 0.0005 | 0.0034 | 0.0005 | 0.0035 | 0.0005
MS0205.7+3509 | 0.0097 | 0.0014 | 0.0100 | 0.0015 | 0.0101 | 0.0015
MS1229.2+6430 | 0.0089 | 0.0007 | 0.0095 | 0.0008 | 0.0096 | 0.0008
NGC 1313 | 0.0123 | 0.0013 | 0.0124 | 0.0013 | 0.0127 | 0.0013
NGC 4278 | 0.0090 | 0.0009 | 0.0094 | 0.0010 | 0.0096 | 0.0010
NGC 5204 X$-$1 | 0.0062 | 0.0010 | 0.0063 | 0.0010 | 0.0068 | 0.0011
NGC 5204 X$-$1 | 0.0099 | 0.0012 | 0.0096 | 0.0012 | 0.0099 | 0.0012
NGC 5252 | 0.0008 | 0.0000 | 0.0008 | 0.0000 | 0.0008 | 0.0000
NGC 6251 | 0.0161 | 0.0020 | 0.0168 | 0.0020 | 0.0172 | 0.0021
NGC 7172 | 0.0008 | 0.0030 | 0.0010 | 0.0034 | 0.0009 | 0.0033
PG1351$+$64 | 0.0025 | 0.0001 | 0.0026 | 0.0001 | 0.0027 | 0.0001
PG1407$+$265 | 0.0038 | 0.0003 | 0.0040 | 0.0004 | 0.0041 | 0.0004
PKS0237$-$23 | 0.0036 | 0.0004 | 0.0036 | 0.0004 | 0.0039 | 0.0004
PKSB1334$-$127 | 0.0078 | 0.0012 | 0.0082 | 0.0013 | 0.0082 | 0.0012
RBS 1055 | 0.0083 | 0.0008 | 0.0087 | 0.0009 | 0.0089 | 0.0009
RBS 1423 | 0.0074 | 0.0007 | 0.0076 | 0.0006 | 0.0076 | 0.0007
RX J0136.9$-$3510 | 0.0103 | 0.0008 | 0.0102 | 0.0008 | 0.0105 | 0.0009
RXJ0228$-$40 | 0.0053 | 0.0001 | 0.0053 | 0.0002 | 0.0057 | 0.0002
RX J0944.5$+$0357 | 0.0056 | 0.0009 | 0.0055 | 0.0009 | 0.0059 | 0.0009
UZ LIB | 0.0708 | 0.0116 | 0.0708 | 0.0116 | 0.0762 | 0.0125
V410 Tau | 0.0726 | 0.0145 | 0.0754 | 0.0150 | 0.0762 | 0.0152
V410 Tau | 0.1199 | 0.0145 | 0.1086 | 0.0131 | 0.1144 | 0.0138
VB 50 | 1.0000 | 0.1343 | 0.9509 | 0.1280 | 0.9958 | 0.1339
VV Sco | 0.1161 | 0.0184 | 0.1163 | 0.0184 | 0.1175 | 0.0186
VV Sco | 0.1157 | 0.0167 | 0.1193 | 0.0172 | 0.1218 | 0.0176
X Comae | 0.0052 | 0.0001 | 0.0055 | 0.0001 | 0.0056 | 0.0001
X Comae | 0.0070 | 0.0001 | 0.0072 | 0.0001 | 0.0075 | 0.0001
X Comae | 0.0051 | 0.0001 | 0.0052 | 0.0001 | 0.0053 | 0.0001
X Comae | 0.0055 | 0.0001 | 0.0056 | 0.0001 | 0.0058 | 0.0001
X Comae | 0.0072 | 0.0001 | 0.0074 | 0.0002 | 0.0075 | 0.0002

Table A.2: 2XMM Concordance Fluxes – Medium BandaaFluxes are normalized to
0.000701 photons cm-2 s-1.

Table A.3: 2XMM Concordance Fluxes – Hard BandaaFluxes are normalized to
0.00223 photons cm-2 s-1.

Table A.4: XCAL Concordance Fluxes – Soft BandaaFluxes are normalized to 0.126
photons cm-2 s-1.

Table A.5: XCAL Concordance Fluxes – Medium BandaaFluxes are normalized to
0.0156 photons cm-2 s-1.

Table A.6: XCAL Concordance Fluxes – Hard BandaaFluxes are normalized to
0.0154 photons cm-2 s-1.

Table A.7: Capella grating observations with Chandra

## References

  * Beuermann et al. (2006) Beuermann, K., Burwitz, V., & Rauch, T. 2006, A&A, 458, 541
  * Chen et al. (2019) Chen, Y., Meng, X.-L., Wang, X., et al. 2019, Journal of the American Statistical Association, 114, 1018.  https://doi.org/10.1080/01621459.2018.1528978
  * Drake et al. (2006) Drake, J. J., Ratzlaff, P., Kashyap, V., et al. 2006, in Proc. SPIE, Vol. 6270, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 62701I
  * Efron & Morris (1972) Efron, B., & Morris, C. 1972, Biometrika, 59, 335
  * Efron & Morris (1973) —. 1973, Journal of the American Statistical Association, 68, 117
  * Fruscione et al. (2006) Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in Proc. SPIE, Vol. 6270, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 62701V
  * Gabriel et al. (2004) Gabriel, C., Denby, M., Fyfe, D. J., et al. 2004, in Astronomical Society of the Pacific Conference Series, Vol. 314, Astronomical Data Analysis Software and Systems (ADASS) XIII, ed. F. Ochsenbein, M. G. Allen, & D. Egret, 759
  * Guainazzi et al. (2015) Guainazzi, M., David, L., Grant, C. E., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 047001
  * Ishida et al. (2011) Ishida, M., Tsujimoto, M., Kohmura, T., et al. 2011, PASJ, 63, S657
  * Jansen et al. (2001) Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1
  * Jethwa et al. (2015) Jethwa, P., Saxton, R., Guainazzi, M., Rodriguez-Pascual, P., & Stuhlinger, M. 2015, A&A, 581, A104
  * Kashyap & Drake (2000) Kashyap, V., & Drake, J. J. 2000, Bulletin of the Astronomical Society of India, 28, 475
  * Kashyap et al. (2008) Kashyap, V. L., Lee, H., Siemiginowska, A., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7016, Observatory Operations: Strategies, Processes, and Systems II, ed. R. J. Brissenden & D. R. Silva, 70160P
  * Kettula et al. (2013) Kettula, K., Nevalainen, J., & Miller, E. D. 2013, A&A, 552, A47
  * Kirsch et al. (2005) Kirsch, M. G., Briel, U. G., Burrows, D., et al. 2005, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5898, Proc. SPIE, ed. O. H. W. Siegmund, 22–33
  * Lee et al. (2011) Lee, H., Kashyap, V. L., van Dyk, D. A., et al. 2011, ApJ, 731, 126
  * Madsen et al. (2017a) Madsen, K. K., Beardmore, A. P., Forster, K., et al. 2017a, AJ, 153, 2
  * Madsen et al. (2017b) Madsen, K. K., Forster, K., Grefenstette, B. W., Harrison, F. A., & Stern, D. 2017b, ApJ, 841, 56
  * Madsen et al. (2019) Madsen, K. K., Natalucci, L., Belanger, G., et al. 2019, arXiv e-prints, arXiv:1901.00934
  * Madsen et al. (2020) Madsen, K. K., Terada, Y., Burwitz, V., et al. 2020, arXiv e-prints, arXiv:2001.11117
  * Nevalainen et al. (2010) Nevalainen, J., David, L., & Guainazzi, M. 2010, A&A, 523, A22
  * Plucinsky et al. (2017) Plucinsky, P. P., Beardmore, A. P., Foster, A., et al. 2017, A&A, 597, A35
  * Read et al. (2014) Read, A. M., Guainazzi, M., & Sembay, S. 2014, A&A, 564, A75
  * Schellenberger et al. (2015) Schellenberger, G., Reiprich, T. H., Lovisari, L., Nevalainen, J., & David, L. 2015, A&A, 575, A30
  * Strüder et al. (2001) Strüder, L., Briel, U., Dennerl, K., et al. 2001, A&A, 365, L18
  * Toor & Seward (1974) Toor, A., & Seward, F. D. 1974, AJ, 79, 995
  * Tsujimoto et al. (2011) Tsujimoto, M., Guainazzi, M., Plucinsky, P. P., et al. 2011, A&A, 525, A25
  * Turner et al. (2001) Turner, M. J. L., Abbey, A., Arnaud, M., et al. 2001, A&A, 365, L27
  * Watson et al. (2009) Watson, M. G., Schröder, A. C., Fyfe, D., et al. 2009, A&A, 493, 339
  * Xu et al. (2014) Xu, J., van Dyk, D. A., Kashyap, V. L., et al. 2014, ApJ, 794, 97