CALIBRATION OF THE NuSTAR HIGH-ENERGY FOCUSING X-RAY TELESCOPE

The Crab is a center filled pulsar wind nebula (PWN) powered by a pulsar with a double peaked profile of period P ∼ 33 ms. It has served as the primary celestial calibration source for many hard X-ray instruments because of its brightness, relative stability, and simple power-law spectrum over the band from 1 to 100 keV (Kirsch et al. 2005). The Crab, however, is too bright for reliable pileup corrections for most CCD based focusing X-ray instruments, and has been replaced with fainter sources, such as the Crab-like PWN G21.5+0.9 (Tsujimoto et al. 2011), the blazar PKS 2155-304 (Ishida et al. 2011) or the quasar 3C 273. The latter two are variable and only serve for cross-calibration among observatories, while G21.5+0.9 would require unfeasible long integration times for NuSTAR to obtain the required statistics to calibrate the response at 80 keV. The Crab therefore remains the best choice for the internal calibration in the X-ray band covered by NuSTAR.

The spatially phase-averaged integrated spectrum of the Crab Nebula+pulsar in the 1–100 keV X-ray band has been well-described by a power-law with photon index of ${\rm{\Gamma }}\sim 2.1$ (RXTE, BeppoSAX, EXOSAT, INTEGRAL/JEM_X; Kirsch et al. 2005). Above 100 keV the hard X-ray instruments (INTEGRAL/SPI/ISGRI, CGRO) measure a softer index of ${\rm{\Gamma }}\;\sim 2.20-2.25$, and below 10 keV the instruments with CCD detectors a harder spectrum (${\rm{\Gamma }}\lt 2.1$). The softening above 100 keV primarily comes from the curvature of the pulsed spectrum (Kuiper et al. 2001). The hardening in CCD X-ray instruments may be due to photon pile-up and although models exist for these instruments to deal with the pile-up the Crab usually still challenges the pile-up corrections and requires special non-standard reductions. Weisskopf et al. (2010) did a study with RXTE/PCA, XMM, and ASCA looking for deviations from a power-law in the 0.2–50 keV energy range and concluded that within the precision of the available instrumentation there is no detectable bend in the phase-averaged integrated spectrum of the Crab. Over the 16 years RXTE was operational and regularly monitoring the Crab, the spectral index was seen to vary by a peak-to-peak variation of ${\rm{\Delta }}{\rm{\Gamma }}\sim 0.025$ (Shaposhnikov et al. 2012). This variation is consistent with the observed spread between instruments, but it is slow and on average over the 16 years the Crab has remained at ${\rm{\Gamma }}=2.1$.

Because the non-piled-up instruments covering the 1–100 keV band agree on a photon index of ${\rm{\Gamma }}=2.1\pm 0.02$ and do not measure any curvature in the spectrum across this band, we have calibrated the effective area against a Crab index of ${\rm{\Gamma }}=2.1$. We model the differential photon spectral model of the phase-averaged integrated spectrum of pulsar+nebula as a simple absorbed power law,

$$\begin{eqnarray}&&\displaystyle \frac{{dN}(E)}{{dE}}=\mathrm{tbabs}(E){{NE}}^{-{\rm{\Gamma }}}\mathrm{ph}\;{{\rm{s}}}^{-1}\;{\mathrm{cm}}^{-2}\;{\mathrm{keV}}^{-1}.\end{eqnarray} \tag{ 1 }$$

Here E is the photon energy, tbabs is the interstellar absorption using Wilms abundances (Wilms et al. 2000) and Verner cross-sections (Verner et al. 1996), Γ is the power-law photon index, and N is the normalization factor. When ${\rm{\Gamma }}=2.1$ and N = 10 keV⁻¹ cm⁻² s⁻¹ at 1 keV, we refer to this as the canonical Crab model.

In terms of flux stability the Crab has been tracked in great detail over the last decade. It remains largely stable, albeit with flux changes as much as 7% (Wilson-Hodge et al. 2011) across the 10–100 keV bandpass. This corresponds to a decay in the flux of 3.5% per year over the period it has been observed. The NuSTAR observations span approximately a year, but due to the optical axis (OA) knowledge (to be discussed in Section 2.4) we expect flux differences between repeated observations to be on the same order (5%), and therefore the intrinsic variability of the Crab is not an issue. In addition, after adjustment of the effective areas of the two mirror modules assuming the canonical Crab model, we adjust the overall normalization by 15% and use the two targets PKS 2155-304 and 3C 273 observed simultaneously with Chandra, NuSTAR, Swift, Suzaku, and XMM-Newton to confirm the final global normalization of our responses with respect to the other observatories.

The NuSTAR Crab calibration campaign measured the Crab spectrum at different locations on the detector modules at various off-axis angles. This resulted in the 39 individual observations listed in Table 1.

Despite the angular extent of the Crab of ∼120'' $\times \;100$'', we treat it as a point source since the center of the nebula where the pulsar resides dominates the emission. The difference in the effective area response between using a point source effective area centered on the pulsar and a count-rate adjusted extended effective area is only ∼1%–2%. This difference is largely constant as a function of energy and is therefore an acceptable error. Treating the Crab as a point source allows us to use the built-in pipeline corrections of the Aperture Stop (AS), Ghost Ray (GR), and PSF corrections (the details of these components will be discussed in Section 2.3.2). These corrections are not applicable to extended source responses, and since they are much larger than the error introduced by assuming the Crab is a point source, they are extremely important to include.

We extracted counts from the largest possible source region, typically 200'', that did not extend beyond the edge of the detector. At this radius about 95% of all photons are contained in the source extraction region. We used the NuSTARDAS pipeline version v1.3.0 to reduce the data, generate the detector response, RMF, and to apply GR, AS and PSF corrections to the base effective area.

Prior to launch, we used ground calibration and modeling to produce a physics-based optics ray-trace model (Westergaard et al. 2012) to generate a set of base effective area and PSF files. Figure 1 shows the ratio of the canonical Crab model folded through the base responses to one of the observations. This comparison shows several features that indicate inaccuracies in the ground calibration.

• 1.  
  Below 10 keV, the response turns over exponentially due to uncertainties in the thickness of the Platinum (Pt) electrode on the detector surface.
• 2.  
  At around 10 keV, the Tungsten (W) L-edge causes residuals (See Table 2). The exact shape of the W edge depends on the optical constants of the multilayer coatings, which are very sensitive to the crystalline structure and density of the W. Between 10 and 14 keV, the Pt and W L-edge complex, along with the top few coatings of the multilayer stack, create residuals.
• 3.  
  Below 20 keV, X-rays reflect off the optics by total external reflection. Above this energy the multilayer coatings become effective and the resonant features, enhancing the reflectivity, are sensitive to the structure of the multilayer that we can only approximate because the coatings are not precisely tabulated for each individual mirror segment of which there are more than 6000.
• 4.  
  Between 30 and 40 keV, the effective area deviates gradually from the model areas. The exact reason for this discrepancy is not known.
• 5.  
  At 69.4 keV, the W K-edge creates a residual for the same reasons as the L-edge. At 78.4 keV the Pt K-edge is just evident.

Zoom In Zoom Out Reset image size

In the following sections we model the detected counts in a given instrumental pulse height bin, $C({PI},\theta )$, according to the equation

$$\begin{eqnarray}&&C({PI},\theta )=\int \displaystyle \frac{{dN}(E)}{{dE}}R({PI},E,\theta ){dE}.\end{eqnarray} \tag{ 2 }$$

Here ${dN}(E)/{dE}$ is the model differential photon spectrum of the observed target as a function of incident photon energy, E (flux in units of photons s⁻¹ cm⁻²), and $R({PI},E,\theta )$ is the response matrix that captures the photon indicent in a given pulse height bin, PI (units of cm²), at an off-axis angle, θ. In practice, this integral is approximated as a finite sum by sampling $R({PI},E,\theta )$ on a fine grid. The off-axis angle, θ, is the angle of incoming X-rays with respect to the OA of the mirrors. As the OA moves with respect to the detector position during an observation, the modeled response is sampled on a finite time grid and then summed for a given exposure.

Zoom In Zoom Out Reset image size

As is typical for X-ray astronomy missions, the response matrix is divided into two components,

$$\begin{eqnarray}&&R({PI},E,\theta )=\mathrm{RMF}({PI},E)A(E,\theta ).\end{eqnarray} \tag{ 3 }$$

RMF(PI, E) is known as the redistribution matrix, which contains detector quantum efficiency and resolution effects, and is unitless (a fraction between 0 and 1, where a value of 1 indicates 100% quantum efficiency). The quantity $A(E,\theta )$ is the effective area, also known as the ancillary response function (ARF), which captures the effective area of the mirror optics, as well as several other factors unique to NuSTAR,

$$\begin{eqnarray}A(E,\theta ) & = & {A}_{0}(E,\theta )V(E,\theta ){detabs}(E)\mathrm{GR}(E,\theta )\\ & & \mathrm{AS}(E,\theta )C(E,\theta ).\end{eqnarray} \tag{ 4 }$$

Here A₀ is the modeled effective area of the mirror segments estimated from the physics-based ray-tracing simulations, and V is the geometric vignetting function, also based on the ray-tracing simulations. The quantities detabs, GR, and AS are the detector dead layer absorption, GR correction, and AS correction, respectively, described in the next sections. C is an empirically derived correction factor, which we will discuss at length in Section 2.3.3.

2.3.1. Detector Absorption Correction

2.3.2. Ghost Rays and Aperture Correction

The NuSTAR observatory has two benches, the Focal Plane Bench (FB) and the Optics Bench (OB) separated by 10.14 m, which move relative to each other due to motions of the mast that connects them and causes the OA to travel across the FPMs. Because the optics act as thin lenses, this track is not equivalent to the path made by the source; the rotational component of a bend in the mast does not project into a movement of the source on the detectors, only the translational part does. In the system there are therefore two paths that must be kept track of as a function of time: the OA path, which is typically on the order of 30''–60'' wide, and the source path, which can be considerably larger and is mainly driven by the pointing stability of the spacecraft bus. The calculation of off-axis angle is carried out as a function of time, and since the resolution of the response files is 10'', the off-axis aspect histogram is binned at this resolution.

Three ASs are attached to each FPM with the top acting as the limiting aperture. The diameter of the aperture opening is 58 mm and the total thickness of the solid part of the aperture is 2 mm, layered into 0.75 mm Al, 0.15 mm Cu and 1.1 mm Sn. The purpose is to block the diffuse background by limiting the FOV of the detectors to the open sky. However, this comes at the cost of a reduction in effective area at off-axis angles due to the clipping of the edges of the optical path between detector and optic. Since the AS is fixed with respect to the detector, and the OA is moving independently, the correction is strongly observation dependent. This clipping is a purely geometrical effect and can be calculated by a ray-trace. We have assembled the result as a function of off-axis angle and distance to AS center in a look up table used by the software that creates the ancillary response file for a given observation (numkarf ). An example of the magnitude of this correction is shown in Figure 2. The correction is azimuthally dependent since the AS and the OA are not centered, and the correction has a spectral dependence, preferentially cutting away low energy photons since the majority of these come from the outer shells of the optic and thus are more prone to being blocked than light focused by the inner shells. Figure 2 shows that the correction is only important for off-axis angles > 3'.

Zoom In Zoom Out Reset image size

The NuSTAR optics are based on a Wolter-I conical approximation, which is a grazing incidence, double mirror system. A properly focused photon will reflect twice off the optics before exiting. However, it is possible for a photon to reflect only once by either the upper conical section or the lower conical section. This occurs at either very shallow or very steep angle and we call these single-bounce photons GR. The pattern of GR is very distinct, as shown in Figure 3, and axi-symmetric. Figure 3 shows the pattern for a series of off-axis angles, starting with the source being on-axis and hitting the center of the detector. In this case the GR are uniformly distributed in a circular halo, which we never see in NuSTAR because it is blocked by the AS. As the source moves off-axis and the spot wanders toward the left side of the detector, the GR can be observed entering from the left until around 6' they overlap with the double bounce photons. Since the GR are now inside the typical extraction region, the net effect is an increased effective area. Like the properly reflected photons, these are also subject to being partially blocked by the AS, and so this must also be taken into account. Due to the complexity of the GR pattern, the corrections were only derived for circular regions. Figure 4 shows an example of the correction, which, like the AS correction, is a function of azimuth, but also depends on source extraction region size and energy. The suppression of the low energies is because of the AS correction to the GR.

Zoom In Zoom Out Reset image size

The energy dependent correction is applied to the ancillary file in numkarf, and because of the complexity of the GR correction files and their dependency on extraction region, this correction is only available for circular and annular regions, and only for point sources.

2.3.3. Empirical Correction and Vignetting

After adjusting for all the above terms (detabs, AS, GR), there still remain systematic residuals in the Crab spectrum which are not accounted for, with typical residuals of a few percent for on axis targets and 10%–20% for off axis targets.

Although it is desirable to have a completely physics-based response matrix derived from a priori information, this is not usually achievable with limited project resources. We found that the residual trends as a function of energy, E, and off axis angle, θ, are very reproducible, and therefore we elected to adjust the response by an empirical correction factor $C(E,\theta )$, using the canonical Crab spectral model. C is a fractional quantity, where a value of 1 indicates no correction. We note that C is a function of energy E (not pulse height), so its effect is equivalent to modifying the effective area term in the ARF, and its effects must be folded through the response matrix to be seen in the model pulse height spectrum. As a purely practical matter, $C(E,\theta )$ is stored along with the vignetting function in the NuSTAR CALDB (i.e., a new vignetting function, ${V}^{\prime }(E,\theta )=V(E,\theta )C(E,\theta )$).

We investigated several functions for $C(E,\theta )$ and decided to use a piece-wise linear interpolation function, in both energy and off-axis angle. This function consists of piece-wise linear segments between fixed control points E_ci in energy. We chose the control points of the spline based on the known features. Figure 5 marks the locations of the nodes with 27 control points between 3.0 and 80 keV. The function is completely specified by the value at each of the control points, which are considered the parameters of the correction model. The angular dependence is treated as linear as well, as described further below. We investigated using Chebyshev polynomials or cubic splines as an interpolation function, but encountered difficulty with ringing and other unconstrained behavior in the neighborhood of the control points. The interpolation approximation results shown here worked well for our purposes to within available statistics of the data. An improved model may be found in the future, but the present correction achieves residuals at the few percent level over almost the entire NuSTAR energy range.

Zoom In Zoom Out Reset image size

We derived the correction spectrum in two stages. First, we fit individual Crab spectra to derive the energy dependence for several off-axis angles. Second, we used linear interpolation to model the angular dependence to derive a single global correction model.

For the first stage of analysis, we combined the observations in 1' bins to improve statistics, resulting in seven spectra between 0' (on-axis) and 7' off-axis. We fit the parameters of the correction model while holding the canonical spectral model fixed, using standard ${\chi }^{2}$ fitting, and requiring at least 200 counts per bin. The result is a linear model for each of the seven spectra. We note that for the highest control points at 77 and 80 keV, there were not enough counts to reliably estimate the parameter values at both points independently, so the value of the response at the 80 keV control point was fixed manually to be 90% of the value at 77 keV. In the final result, we claim calibrated flux levels only between 3 and 78 keV.

The corrections show an angular dependence nearly linear with off-axis angle (Figure 6). This indicates that the physics-based ray-tracing captured most of the angular dependence of the reflectivity of the mirrors, except for a residual trend. We therefore treat the angular dependence of the correction as linear between on-axis and 7' off-axis. For off-axis angles above 7' (the outer parts of the focal plane) we lack the Crab data to constrain the model, so we cannot verify if the linear trend continues. To be conservative, we do not extrapolate the linear trend beyond 7', but simply hold the function fixed at the 7' value.

Zoom In Zoom Out Reset image size

In a second stage of analysis, we fitted the individual correction values at each control point independently with a linear function, using the standard errors from the individual fits, to derive one global correction function for each FPM. We stress that this global correction is an ensemble-derived model using all of the Crab observations; no single Crab spectrum was used to fix the final global model parameters. Thus, the final global correction model consists of 28 correction parameters on-axis and 28 parameters 7' off-axis, for a total number of model parameters of 56 for each FPM. These are derived from 7 × 26 (=182) individual correction spectra, which were in turn derived from 7 × 940 (≃6600) individual Crab count bins for each FPM. Figure 5 shows the correction function $C(E,\theta )$ as a function of energy E and off-axis angle θ for both FPMs. The systematic nature of the adjustments and some of the features mentioned in Section 2.3 are evident.

As noted above, the correction factor $C(E,\theta )$ and geometric vignetting function derived from ray tracing $V(E,\theta )$ were combined into a single effective vignetting function ${V}^{\prime }(E,\theta )$, and stored in the NuSTAR CALDB. This effective vignetting function is shown in Figure 7 (bottom), along with the resulting total effective area for each module (top).

Zoom In Zoom Out Reset image size

We reprocessed the Crab data with the adjusted calibration files to verify that the correction was successful. Figure 8 shows the residual spectrum for all Crab data within 3' of on-axis for both modules. The residuals are typically better than ±2% up to ∼40 keV. Between 40 and 80 keV, residuals are dominated by counting statistics of the Crab data, but are typically 5%–10%, so the systematic errors are less than that.

Zoom In Zoom Out Reset image size

The final correction achieves few percent residuals over the entire NuSTAR energy range of 3–78 keV, and is applicable for off-axis angles of 0'–7'. We do not advise to fit spectra above 78 keV.

Having calculated the correction function, as a consistency check we reprocessed all individual Crab data sets to determine if the specified Crab spectral parameters were recovered. Although the Crab spectral shape and normalization were enforced at the ensemble level during the correction process, we did not required that individual spectra conform, nor did we require that FPMA and FPMB produce the same results. Any deviations at the individual observation level would indicate unmodeled systematic errors, as well as errors in the assumption of a constant Crab. To estimate this, we computed sample biases and standard deviations for various subsets of observations. This is a measure of how repeatable NuSTAR observations are.

We fit each of the spectra with the canonical spectral shape, but allowed the parameters to vary. The parameters were the relative power-law flux normalization and the power law photon index, Γ, while the neutral absorption parameter, ${N}_{{\rm{H}}}$, was fixed. We also performed a separate fit where in addition ${N}_{{\rm{H}}}$ was allowed to vary. Results from FPMA and FPMB were recorded separately.

We recorded the individual parameters and computed statistics of various subsets of observations. The results are shown in Table 4. In addition to the sample variation among "All" observations, we also examined subsets where the target was near on-axis and far off-axis, and on individual detectors (DET0-3). These designations are not mutually exclusive. In a few cases the Crab straddled two detectors, and these are shown separately. We also computed the sample variations of the differences between FPMA and FPMB. For relative flux, we computed the sample bias as well as the standard deviation; for photon index and absorption, we only list the sample standard deviation. The results are also shown graphically in Figure 9 as a function of off-axis angle.

Zoom In Zoom Out Reset image size

Table 4.  Measurement Repeatability Statistics

Subset	Mean	Std. Dev.
Relative Flux
All	−0.4%	2.6%
On-axis (θ < 3')	−0.9%	3.0%
Off-axis (θ > 3')	−0.9%	3.0%
DET0	−0.3%	2.5%
DET1	−0.7%	2.5%
DET2	+0.8%	1.4%
DET3	−0.0%	3.4%
DET0+3 Straddle	−1.7%	4.0%
FPMA–FPMB Diff.	−0.0%	2.6%
Power-law Photon Index, Γ
All	⋯	0.009
FPMA–FPMB Diff.	⋯	0.012
Absorption, ${N}_{{\rm{H}}}$(10${}^{22}\;{\mathrm{cm}}^{-2}$)
All	⋯	0.150

Download table as:  ASCIITypeset image

The results demonstrate that average relative flux biases are small, typically <1%. Scatter in the relative fluxes is typically 2%–3% (1σ). There are no large flux offsets between FPMA and FPMB during the same observation; in the worst case the difference was 3.8% in this sample. There are no strong biases as a function of off-axis angle or detector number. The power law index was stable to about ±0.01 (1σ); this can be taken as an approximate estimate of repeatibility errors in spectral slope. For neutral absorption, NuSTAR is not particularly sensitive to the Crab's hydrogen column; the sample standard deviation is comparable to the canonical value itself.

The fluctuations in flux and slope are understood, and are consistent with the estimated error in OA location. The knowledge of the OA location is accurate to ∼30'', and because the OA is constantly moving due to the mast motions, there will be situations where the agreement between effective areas are better than in others. The flux level is also affected by the understanding of the PSF (to be discussed in the next section), and the fidelity of the instrument map: gaps, dead pixels, etc. The gap width between detectors was measured on the ground by scanning an X-ray beam across the gaps, but due to systematic errors getting the exact physical gap width is challenging, making it difficult to correct the flux exactly in the cases where the core of the PSF falls, or crosses the gap.

Observed differences between the two modules of up to 5% are not uncommon and are simply a result of the OA configuration during the time of observation. In extension, observed flux differences less than 5% between different observations of the same target, should not be interpreted as intrinsic to the source.

Calculating the correct effective area requires a good understanding of the PSF. Prior to launch a ray-trace model of the PSF was developed and tested against limited ground calibration scatter data from individual flight mirrors (Westergaard et al. 2012). This original set of calibration files formed the basis of pre-launch PSFs determined every 05 out to an off-axis angle of $\theta =8\buildrel{\,\prime}\over{.} 5$.

To conduct the in-orbit calibration of the PSF, we observed the strong point sources listed in Table 5. The main point of these observations was to establish the corrections required to the ground based estimates. To this end, we extracted events from the point source observations in several energy ranges (3–4.5, 4.5–6, 6–8, 8–12, 12–20, and 20–78 keV) and produced radial profiles centered at the peak of the intensity map. Table 6 shows the half-power diameters (HPD) for the nine on-axis observations.

Table 6.  Measured Half-power Diameter for the On-axis Observations in Table 5

Energy	FPMA				FPMB
	HDP	${\sigma }_{\mathrm{HPD}}$	min.	max.	HPD	${\sigma }_{\mathrm{HPD}}$	min.	max.
(keV)	('')	('')	('')	('')	('')	('')	('')	('')
3–4.5	70.3	2.4	66.7	75.5	65.6	2.4	62.9	69.6
4.5–6	67.1	1.0	64.7	67.7	62.6	1.2	60.9	64.7
6–8	64.7	1.0	62.8	65.7	60.7	1.4	58.8	63.7
8–12	63.5	1.1	61.8	64.7	59.5	1.5	57.9	62.8
12–20	63.4	1.1	61.8	64.7	60.3	1.2	58.8	62.8
20–79	63.4	1.0	61.8	64.7	62.4	1.5	60.8	65.7

Download table as:  ASCIITypeset image

We compared these radial profiles with the ray-trace model. Since the off-axis angle changes during each observation due to the mast motions, we weighted the ray-traced PSF with the aspect solution before comparing. We fitted the observed radial profiles to the aspect-weighted ray-trace PSF and found that: (i) the central cores of observed radial profiles are slightly broader than that of the model PSF, (ii) there is an unaccounted for wing at $R\gt 100^{\prime\prime}$ in the radial profiles. The broadening of the central core does not affect the FWHM and is an effect of incomplete astrometric aspect reconstruction. The origin of the wider wing is unknown.

To take care of the central broadening we convolved the modeled PSF with a Gaussian function. Since the origin of the wing is unknown, we decided to empirically model it by introducing a component on top of the ray-trace PSF. Because background subtraction is not feasible during the fitting, we had to include a background component as well, which is well known and can be modeled accurately (Wik et al. 2014). The backgrounds were in general small compared to the sources, dominating only above $R\gt 500^{\prime\prime}$.

The two-dimensional (2D) model PSF can be expressed by the following formula:

$$\begin{eqnarray*}{\rm{\Psi }}(x,y) & = & {C}_{0}{{\rm{\Psi }}}_{\mathrm{RT}}(x,y)\otimes G(x,y,\sigma )\\ & & +{C}_{1}{e}^{-a\sqrt{{x}^{2}+{y}^{2}}}+{C}_{2}B(x,y),\end{eqnarray*}$$

where ${{\rm{\Psi }}}_{\mathrm{RT}}$ is the ray-trace PSF, $G(x,y,\sigma )$ is the Gaussian function, ${e}^{a\sqrt{{x}^{2}+{y}^{2}}}$ is the wing component, $B(x,y)$ is the background model, and ${C}_{\mathrm{0,1,2}}$ are normalization constants.

To fit the radial PSF profiles, we first convolved the ray-trace PSFs with the Gaussian function and weighted them using the aspect solution of each observation to produce an aspect-weighted PSF. The wing model was added to construct an effective PSF, and we multiplied it with the exposure map for the observation. Since the point source exposure map is almost flat, the overall shape of the PSF does not change much. However, multiplying with the exposure map is important for some localized effects such as detector gaps.

The fitting parameters are the amplitude of the ray-trace PSF (${C}_{0}$), width of the Gaussian convolution core (σ), amplitude and decay constant for the exponential wing (${C}_{1}$ and a), and amplitude of the background (${C}_{2}$). Since fitting the observations simultaneously was difficult given different aspect histories (i.e., different effective PSFs) and backgrounds, we calculated ${\chi }^{2}$ for various values of fitting parameters for each observation, and found a set of parameters that minimize the combined ${\chi }^{2}$. An example the observed radial profile for a Cyg X-1 observation (obs. ID: 10002003001) and the best-fit PSF model is shown in Figure 10.

Zoom In Zoom Out Reset image size

From the analysis, we found that the PSF sharpens—the HPD decreases—with energy. As a function of increasing energy we find a smaller core-convolution width and narrower wing, as well as a relatively larger amplitude for the core component (see Figure 11). The sharpening saturates at ∼10 keV. The differences in the enclosed count fractions at different energies are modest, with a maximum of ∼3% between the lowest and the highest energy bands for a ∼1' extraction radius region as shown in Figure 12. However, since this is significant enough for small extraction regions, we include the energy dependent evolution into the ARF generation, linearly interpolating between energies. The additional broadening of the wings of the PSF at energies below 10 keV is likely explained by scattering from surface contaminants.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.1. Off-axis PSF

The FWHM of the NuSTAR PSF does not change with off-axis angle, but the shape of the PSF gradually distorts azimuthally, and the geometrical shadowing of the shells causes the PSF to appear elongated, as shown in Figure 13. Probing the PSF at high off-axis angles is complicated by the lower count-rate, and necessitates a 2D approach where we match the contours of the model PSF to the actual PSF.

Zoom In Zoom Out Reset image size

We followed the same procedure (e.g., convolution, aspect-weighting, etc.) as we did for the on-axis studies to produce the model PSF for the off-axis observations, $\theta \gt 3^{\prime}$ (listed in Table 5). The ray-trace PSF includes the effect of geometrical shadowing (vignetting of photons), and successfully reproduces the observed PSF at the core. However, the wing component, which dominates at R > 100'', did not match with the observed contours at large radii. We parametrized the distortion of the PSF wing with an ellipticity factor and a shadow function, where the form of the shadow function was derived using the geometry of the optics:

$$\begin{eqnarray*}{F}_{\mathrm{shadow}}(x,y,\theta ) & = & 1-A(\theta ){\mathrm{cos}}^{2}({\mathrm{tan}}^{-1}(B(\theta )y/x)),\\ {F}_{\mathrm{wing}}(x,y,\theta ) & = & {e}^{-a\sqrt{{x}^{2}+B(\theta ){y}^{2}}},\end{eqnarray*}$$

where $A(\theta )$ is the shadow factor, $B(\theta )$ is the ellipticity factor, and θ is the off-axis angle in arcminutes. The shadow function was multiplied onto the wing component, which was added to the ray-trace PSF. The parameters were determined using the off-axis point source observations to match the 2D contours at large radii, and we found that $A(\theta )=0.025\theta$ and $B(\theta )=1-0.025\theta$ made the contours match at large radii (see Figure 13).

With these correction factors we generated new energy-dependent 2D PSF files for the NuSTAR CALDB at 3–4.5, 4.5–6, 6–8, 8–12, 12–20, and 20–78 keV.

3.1.2. PSF Stability Over Time

Changes in PSF over time can occur from out-gassing of the epoxy used to bond glass substrates, and/or by temperature gradients applied to the optics by the Sun. Both have been intensively studied with ground experiments, which showed that the PSF should be stable over ∼10 years. Out-gassing is expected to be strongest in the first years after launch, and we compared observations of Vela X-1 (Fürst et al. 2014, Table 5) taken ∼300 days apart to look for broadening in the PSF that might suggest contamination or out-gassing had occurred.

We constructed radial profiles for the three observations of Vela X-1 (Fürst et al. 2014) for the energy bands used above, and compared the radial profiles directly for each radial bin. We subtracted the late-time profiles from the reference profile (T = 0, obs. ID 10002007001) in order to see if there was a significant change in the radial profiles. The radial profiles and the differences in the 3–4.5 keV band are shown in Figure 14. We find there is no significant change in the radial profiles.

Zoom In Zoom Out Reset image size

The cross-calibration observations are listed in Table 8 and were processed using the respective standard pipelines.

For Chandra we used CIAO 4.6.1 and CALDB 4.6.1.1. For 3C 273 the data were taken with gratings configuration ACIS+HETG and reprocessed using the CIAO chandra_repro reprocessing script. We combined orders 1–3 for the HEG and MEG arm separately, and binned the data at 30 counts. For PKS 2155-304 the data were taken with gratings configuration ACIS+LETG. We fit orders m1 and p1 simultaneously.

For NuSTAR we used HEAsoft 6.15.1 and CALDB 20131223. The data were processed with all standard settings and source counts extracted from a 30'' radius circular region for both 3C 273 and PKS 2155-304. Background was taken from the same detector. We did not combine modules.

For Suzaku we used CALDB: HXD (20110913), XIS (20140203). For 3C 273 the observation was taken in 1/4 window mode, and we used 100 arcsec radius circular regions for the front-illuminates, FI, detectors (XIS0,3) and a 140 arcsec region for the back-illuminated, BI, detector (XIS1), such that the regions were centered on the source, but were restricted to the operational portions of the detectors. For PKS 2155-304 the observation was also taken in 1/4 window mode. We extracted from 115'' radius circular region for the FI and 135'' from the BI.

For Swift we used HEAsoft 6.15.1 and XRT CALDB 2014 February 04. The data for 3C 273 were taken in "PHOTON" mode and PKS 2155-304 in "WINDOW TIMING" mode. Both were reduced using xrtpipeline. Spectra for 3C 273 were extracted from an annulus region, inner radius 5'' and outer radius 30'' to correct for pileup. The two observations for 3C 273 were combined.

For XMM-Newton we used SAS v. 13.5.0 with CALDB 2014 January 31. The data for 3C 273 were taken in "Small Window" mode, and to corrected for pileup in the MOS we extracted counts form an annulus region with inner radius of 15'' and outer radius of 45''. For the PN we extracted from a circular region of radius 45''. For PKS 2155-304 the observation was also taken in "Small Window" mode, and for all instruments we extracted from an annular region of inner radius 10'' and outer radius 36''.

Due to the relative beating of the South Atlantic Anomaly (SAA) passages and occultation periods between the low Earth orbit observatories (NuSTAR, Suzaku, and Swift), we decided to forego strict simultaneity between the observatories. Instead we compare NuSTAR exclusively to each instrument separately and truncate the observations to have matching START and STOP times, ignoring SAA passages and occultation periods.

For the data analysis we use the XSPEC version 12.8.2 analysis software (Arnaud 1996), fitting with cstat statistics (Cash 1979) and present 1σ errors unless otherwise stated.

We fitted the spectra independently in XSPEC with an absorbed pegged power-law (tbabs×pegpwrlw) since it uses as normalization the flux between two energies and breaks the coupling of the regular model between power-law index and normalization at 1 keV. In this manner we can distinguish between slope differences and overall normalization offsets. The hydrogen column was held fixed in all cases at 1.79 × 10²⁰ cm⁻² for 3C 273 and 1.42 × 10²⁰ cm⁻² for PKS 2155-304 (Dickey & Lockman 1990). The flux normalization of the pegpwrlw was calculated between 3 and 7 keV for the soft instruments and 20–40 keV for the hard instruments. We limited the NuSTAR fitting range to match that of the comparison instrument and set the lower energy for the soft instruments to 3 keV, except for Chandra and Swift where it was necessary to go down to 2 keV. While the non-overlapping energy bands may affect the slope measurement, it does not affect the flux calculation. The upper range was 8 keV for Chandra, 8 keV for Swift, 9 keV for XMM-Newton, and 9 keV for Suzaku.

The results are summarized in Table 9 for both sources, and Figure 17 shows the 1, 2, and 3σ contours of the individual fits of the instruments. In all cases there is consistency between the measured slopes when considering the 3σ confidence contours, but obvious slope differences do exist, most notably between NuSTAR and XMM-newton/MOS instruments and Swift/XRT in 3C 273. It is possible that this discrepancy is in part due to the brightness of 3C 273, which caused pileup in both instruments. While excision of the piled-up region is common practice, the errors in the PSF wings of the instruments start to play a role and may have skewed the results.

Zoom In Zoom Out Reset image size

The ratio of fluxes are shown relative to FPMA and NuSTAR was intentionally calibrated to be roughly in the middle of the spread between instruments with Chandra gratings and Suzaku/HXD yielding the highest fluxes (R${}_{\mathrm{FPMA}}\sim 1.10$) and Swift/XRT generally the lowest fluxes between 3 and 7 keV (R${}_{\mathrm{FPMA}}\sim 0.95$). As discussed in Section 2.4 we expect to see differences between FPMA and FPMB due to the restrictions on the absolute OA knowledge and while for PKS 2155-304 the two modules yielded the same flux, during 3C 273 there was an offset of 3%. It is well known that in general XMM-Newton/pn is low compared to MOS (Read et al. 2014), but for 3C 273 the spectrum appears to be unusually low for reasons that are not understood. This and the slope differences will be investigated in K. K. Madsen (2015, in preparation).

To investigate cross-normalization constants, we made simultaneous fits between each instrument, respectively, and FPMA and FPMB assuming constant×tbabs×pegpwrlw. We tied all parameters together and only allowed the constant between them to float to evaluate a potential skewing the slope differences might cause on the cross-calibration constant compared to the ratio of fluxes. The results are summarized in Table 10. In all cases the constants are very similar to the FPMA normalized flux from Table 9. We do, however, strongly encourage users to first evaluate the parameters of individual fits and calculate fluxes before applying a cross-normalization constant and tying fit parameters together. Because of the possible variations in flux between the NuSTAR FPMs, we advise applying a floating cross-normalization constant, but with the caveat that if it takes on extreme values the data reductions might have to be re-examined, and if significant residuals occur after fitting it may be an indication that the source signal to noise is high enough to make inter-instrumental slope differences significant. Large discrepancies between FPMA and FPMB may be due to the interference of a gap, which will not affect the shape of the spectrum, but cause one module to have a lower flux.

Different cross-normalization constants may occur if fluxes for energy ranges other than those presented here are used. This is due to slope differences between observatories, and if simultaneity between observations is not strictly enforced this may also have an effect. In general, however, typical values found for other sources than those discussed above are within 10% as demonstrated in Walton et al. (2014) for Holmberg IX X-1 (XMM-Newton and Suzaku), Brenneman et al. (2014) for IC 4329A (Suzaku), and Baloković et al. (2014) for NGC 424, NGC 1320, and IC 2560 (XMM-Newton and Swift).

We have in the above not included the ∼3% uncertainty on the NuSTAR fluxes, which are due to the constraints on absolute knowledge of the optical axis location as discussed in Section 2.4. Finally we stress that since both targets are variable and observations were of different lengths between instruments, the non-NuSTAR instrument fluxes should not be compared against each other, only with respect to NuSTAR. A forthcoming paper K. K. Madsen (2015, in preparation) of the IACHEC consortium on these same cross-calibration campaigns, will discuss the cross-calibration between all instruments in far greater detail.

A basic description of the NuSTAR timing system was provided in Harrison et al. (2013) and Bachetti et al. (2015). We elaborate upon the system here and describe efforts to characterize the timing performance. The performance of the components is summarized in Table 11, and described in more detail below. Known biases in the system are shown in one column, and are typically removed by the standard analysis. The unmodeled terms are shown in the second column, and are either unknown at the time of writing, or are stochastic. Based on this analysis, we expect the corrected event time stamps to be accurate at the 3 ms level (95%).

The NuSTAR on-board time reference is a temperature-compensated crystal oscillator with a nominal frequency of 24 MHz. The frequency stability of this oscillator is specified to be ±4 ppm over the full operating temperature range. On-board time synchronization of subsystems is done by broadcasting a one pulse-per-second (PPS) electronic signal. The propagation delay of this signal is expected to be less than a few microseconds. The instrument focal plane modules maintain their own internal high resolution timers: seconds are marked using the PPS signal, and sub-seconds are marked internally with 2 μs precision. The delay between when a photon interacts with the detector and when it is time-stamped by the electronics is expected to be a few microseconds, but is certainly less than 0.1 ms. After the spacecraft clock timestamp is applied to an event, it is telemetered to the ground. Science event data is archived with this timestamp.

The spacecraft clock is routinely correlated with ground time references during most ground contacts, which allows measurement of the spacecraft clock offset and trends. The offset is measured essentially by performing one-way ranging of telemetry frames set 1 s intervals. The time of transmission as recorded on the spacecraft and the time of reception on the ground as recorded by the ground station are used to estimate the clock offset. During downlink contacts, the rising edge of the PPS signal initiates the downlink of a telemetry frame. There is a delay of at least 352.75 ms before this packet is transmitted to the ground. During times when the downlink system is under heavy load, this delay can be larger, so clock offset data are filtered to consider only idle downlink times. One-way range delays are modeled using the known spacecraft ephemeris, which is typically accurate to ∼30 km (∼100 lt-μs). The ground stations maintain their own time standards referenced to UTC via GPS, which are typically accurate at the few nanosecond level. However, the USN Hawaii and KSAT Singapore stations only record receive times with 1 ms precision. The RF propagation and Malindi ground station portions of the system have been validated by analysis of Swift data to better than 0.1 ms (Cusumano et al. 2012).

The clock offset trend is examined on the ground. The NuSTAR system is required to maintain on-board time to within 100 ms of UTC, but typically the offset is maintained to less than 30 ms. As the spacecraft clock's enviroment changes, the true frequency also changes. Mission operations staff can adjust the clock divisor setting from the nominal 24 MHz value in order to steer the clock offset toward zero. The clock divisor is the "effective" frequency, or oscillator ticks required to produce a PPS pulse. The mission-average clock divisor offset is ∼14 ppm, with a variation of about ±0.2 ppm (95%). Clock divisor adjustments are visible in both the science event timestamps and the clock offset measurements, and are included in the clock modeling as described below.

We model the clock offset trend as a piece-wise continuous polynomial function. Control points are inserted at points in time when the clock divisor was adjusted, and also where residuals demand extra compliance. The coefficients of the function are recorded in a clock correction file, which is available to the NuSTAR community via the CALDB system. After modeling, the residuals are less than 2 ms (95%), although in cases where the clock measurements are sparse, the retrieval error may be larger than this. This file can be used by the general-purpose FTOOL barycorr in order to correct event times to the solar system barycenter time system, or TDB. As the NuSTAR orbit is known to approximately 30 km, light time errors would be less than 0.1 ms.

We used observations of known X-ray pulsars in order to characterize the actual performance of the system. This includes absolute timing, as well as long term stability of the timing measurement system. We analyzed the extensive observations of the Crab pulsar described above. Observation of a single pulsar allows one to track the mission-long stability of the clock, modulo the ∼33 ms pulse period. We also observed PSR B1509−58, with a longer period of 151 ms to resolve some of this ambiguity. We describe the results of these observations below. We refer observations of known pulsars to their absolute radio ephemerides.

We used all available NuSTAR observations of the Crab pulsar (Table 1), and extracted pulse profiles in the 3–79 keV range. We used the PRESTO¹⁶ (Ransom 2001) suite of programs to obtain an average pulsed profile and to calculate the time-of-arrival (TOA) of the pulses on timescales of ∼5 minutes. We then compared these TOAs to the relevant monthly timing solution for the Crab produced by the Jodrell Bank radio telescope¹⁷ (Lyne et al. 1993) using the Tempo2 software¹⁸ (Edwards et al. 2006). The radio ephemeris is referred to well-known terrestrial time standards. There is an offset between radio and X-ray pulse profiles of ∼300 μs, which we do not consider further in this analysis because it is negligible in light of the current NuSTAR timing performance.

The results are shown in Figure 18. Over the whole history of NuSTAR Crab observations, the residuals of these TOAs show a trend between –2 and +3 ms. This level of residual is indeed mostly inside the $\pm \sim 3$ ms uncertainty range estimated above.

Zoom In Zoom Out Reset image size

In addition to the long term residual trend, shorter term trends are also visible (see inset of Figure 18). These trends repeat on the timescale of the spacecraft orbital period of about 5400 s (0.07 day), and have an amplitude of about 400 μs. This trend is clearly related to the temperature history of the spacecraft clock oscillator. Oscillator instability on these timescales hinders the detection of pulsars with period shorter than 3 ms, because the pulse peak will be smeared over a spacecraft orbital period.

We also measured the TOAs of another pulsar, PSR B1509−58, with a pulse period of ∼151 ms and verified that the TOAs were aligned with those measured by Swift in SW timing mode to 1 ± 2 ms (details can be found in Mori et al. 2014).

In summary, using two pulsars with known absolute ephemerides, the NuSTAR timing system is performing at the accuracy level of 3 ms over the entire mission to date, which is about the expected level of timing uncertainty based on design (Table 11).

Present research indicates that the spacecraft clock frequency variations are driven almost entirely by its thermal environment. As accurate temperature information about the environment is known, it may be possible to make a more accurate clock variation model. This will be the subject of future work.