FAR INFRARED VARIABILITY OF SAGITTARIUS A*: 25.5 hr OF MONITORING WITH HERSCHEL^*

SPIRE data were collected in two 12.75 hr blocks: the first from 2011 August 31 22:04 UT through 2011 September 01 10:51 UT; and the second from 2011 September 01 20:33 through 2011 September 02 9:20. Each interval includes 668 scans across the Galactic Center. Table 1 shows the observation identifiers (ObsIDs) downloaded from the Herschel Science Archive for this work.

We reprocessed the SPIRE data products using the Herschel HIPE software to include the extended flux density gain calibration data products in version 3.1 of the HIPE calibration tree. This step normalizes the response of each bolometer integrated over the beam area, rather than to the peak flux density, which is more appropriate for fields with extended emission. We also chose to include the scan turnarounds in our reprocessing and map making. This option provides additional points on the sky where bolometers make overlapping measurements, increasing the constraints on the calibration algorithms.

For each interval, we concatenated all the Level 1 scans from each SPIRE array for each ObsID array into a single Level 1 context to feed to the HIPE destriper⁸ (Schulz, B., et al. 2016, in preparation). The destriper iteratively determines offsets for all scans crossing the mapped region on the sky. Each scan consists of many detector readouts. The iterations stop when the variances of the readouts within the boundaries of the map pixels cannot be further improved. By running the destriper in "perScan"-mode, individual offsets—we used a 0° polynomial—are fitted for each scan of a given detector, compensating for any long-term variations in the scans. We ran the destriper twice, using the output diagnostic table and destriped scans as inputs for the destriper on the second iteration. This provided small improvements.

Although our relative detector calibration is optimized for extended sources, we produced maps calibrated in units of Jy beam⁻¹. To do this, we binned the scans into groups of 4, and made a single map for each bin, resulting in a time resolution of 4.6 minute per map. We assigned the same sky coordinates to each pixel in each map, taking care to center the location of Sgr A* in the central pixel.

Preliminary review of the maps revealed motion of the flux density distribution with respect to the pixels. This motion is due to insufficiently reconstructed pointing drifts of the telescope that result in inaccurate sky coordinates associated with each bolometer readout. Uncorrected, these drifts limit the precision with which we can calibrate the bolometers and extract light curves.

We solved for pointing offsets as a function of time by shifting each map to best align with the first map produced for each observing interval. Total drifts over the 12.75 hr observing intervals were $\sim 2^{\prime\prime}$ and $\sim 1^{\prime\prime}$ for the observations starting 2011 August 31 and 2011 September 01, respectively. This is consistent with pointing uncertainties given by Sánchez-Portal et al. (2014).

After solving for the best-fit shifts, we updated the coordinates of the SPIRE scans in HIPE and re-ran the destriper and our mapping routine. We iterated the whole process once, and the results showed that our shifts had converged. A small residual drift, $\sim 0\buildrel{\prime\prime}\over{.} 2$, remains in the data.

Due to the bright extended emission of the Sgr A complex at the Galactic Center, and because of the relatively large beam size delivered by Herschel ($18^{\prime\prime}$, $25^{\prime\prime}$, and $36^{\prime\prime}$ at 0.25 mm, 0.35 mm, and 0.5 mm, respectively) Sgr A* is not separated from its surroundings. However, Sgr A* is expected to be the only intrinsically variable source in our maps (there is a magnetar, SGR J1745-2900, that is about $2\buildrel{\prime\prime}\over{.} 4$ away from Sgr A*, but the magnetar was in its quiescent phase when these X-ray/submillimeter observations were carried out; Kennea et al. 2013). Therefore, we extract variability light curves from difference images, subtracting the mean map of each 12.75 hr observation from the 4 minute sub-maps.

We performed photometry on the difference maps by scaling a 1 Jy beam⁻¹ reference point-spread function to best fit our observations. We downloaded the reference beams for each waveband from the SPIRE public wiki.⁹ Specifically, we combined our observed difference maps (D_ij), a variance map created by calculating the variance of all bolometer readouts contributing to a given map pixel (${\sigma }_{{ij}}^{2}$), and the PSF (P_ij) as follows:

$$\begin{eqnarray}&&\bar{f}=\displaystyle \frac{{\displaystyle \sum }_{{ij}}\left(\tfrac{{P}_{{ij}}^{2}}{{\sigma }_{{ij}}^{2}}\right){f}_{{ij}}}{{\displaystyle \sum }_{{ij}}\left(\tfrac{{P}_{{ij}}^{2}}{{\sigma }_{{ij}}^{2}}\right)},\end{eqnarray} \tag{ 1 }$$

which is the inverse variance weighted mean of the scale factor

$$\begin{eqnarray}&&{f}_{{ij}}=\displaystyle \frac{{D}_{{ij}}}{{P}_{{ij}}}.\end{eqnarray} \tag{ 2 }$$

The variance of our measured scale factor can then be computed using

$$\begin{eqnarray}&&{\bar{\sigma }}^{2}=\displaystyle \frac{{\displaystyle \sum }_{{ij}}\left(\tfrac{{P}_{{ij}}^{2}}{{\sigma }_{{ij}}^{2}}\right){({f}_{{ij}}-\bar{f})}^{2}}{{\displaystyle \sum }_{{ij}}\left(\tfrac{{P}_{{ij}}^{2}}{{\sigma }_{{ij}}^{2}}\right)}.\end{eqnarray} \tag{ 3 }$$

We extracted light curves from the location of Sgr A* and several reference locations. Reference light curves should show no intrinsic variability so serve as indicators of time variable systematic problems. We chose reference locations as follows: First, we generated 100 random locations within 2 arcmin of Sgr A*. From that set we excluded any points whose average flux density was not within a factor of two of the average flux density at the position of Sgr A*. We also excluded points where the local spatial gradient had a magnitude that was not within a factor of two of the gradient at the location of Sgr A*. We then searched the remaining locations for a maximal set with no two references within $40^{\prime\prime}$. This last criterion ensures that the 0.5 mm beam does not overlap at the 50% level for any of our reference locations. Using this approach we identified 12 locations on the map for use as references.

We noticed that many of our reference light curves were affected by a small linear trend across the observing interval. This trend was largest in the 0.25 mm band where the average slope was measured to be $\sim -0.02$ Jy h⁻¹. The steepest slope removed from our reference light curves was $\sim -0.1$ Jy h⁻¹. These values are consistent with the trends expected given the residual pointing drift that remains in the maps and the flux density gradients at the locations of our references. Since these drifts strongly affect the appearance of the inter-band cross-correlations (resulting in relatively high power over a large range of lags) we subtract a best-fit line from each light curve. Unfortunately, this correction precludes a meaningful test of whether our light curves are stationary (a time series is stationary if there are no changes in its mean value or variance and there are no periodic components, Chatfield 1989). However, in Section 4 we quantify changes in the variance of our light curves in subintervals of four-hour length.

Our calculated errorbars (Equation (3)) were over estimated for each location, including Sgr A*. This was obvious given the magnitude of the point to point variations in the light curves and the much larger size of the calculated errorbars. The over-sized errorbars result from the way that our variance maps are produced. In HIPE, variance maps are produced by binning all bolometer readouts that occur within a given pixel without respect to where within a pixel a readout occurs. In regions of complex structure, such as Sgr A, spatial flux density gradients will lead to variations in flux density values within a pixel, inflating the variance. To account for this, we scaled the errorbars for each reference light curve to provide a good fit to a constant zero-flux model (reduced ${\chi }^{2}=1$). Typical scale factors were ∼0.3. We took the mean scale factor and applied it to the errorbars for the light curve of Sgr A*. This approach provides empirically accurate errorbars that maintain appropriate relative size as a function of time and location on the map.

Ground-based observations with the SHARC II camera at the CSO provided 0.85 mm monitoring from 2011 September 1 04:25 UT through 09:04 UT, and from 2011 September 2 03:35 UT through 09:00 UT, overlapping each of our Herschel observing intervals. A $3^{\prime}$ field surrounding Sgr A* was observed with Lissajous scanning of the telescope with an amplitude of $100^{\prime\prime}$ and a period of 20 s (Yusef-Zadeh et al. 2006, 2008). On both evenings, the conditions were suitable for observation for the full periods, with clear skies or light cirrus, low wind, and moderate humidity. The zenith atmospheric opacity at 225 GHz was ∼0.14 on September 1 and ∼0.10 on September 2. The telescope focus was monitored and, as needed, adjusted during separate observations of point sources, accounting for the gaps in the light curves; the larger gap around 6:00 UT on September 2 was due to a brief observation of Sgr A* at 0.35 mm which did not yield useful results.

Data analysis, including absolute calibration, followed the method described by Yusef-Zadeh et al. (2009). Sgr A* is not well resolved from surrounding dust emission with the $19^{\prime\prime}$ resolution of CSO at 0.85 mm. This adds ∼1 Jy uncertainty to the absolute flux level of Sgr A* measured at 0.85 mm, but the measurement of variations is much more precise. Uncertainties for each 0.85 mm measurement were derived from the rms in the image, from which the mean image and a Gaussian at the position of Sgr A* have been subtracted.

XMM-Newton data were collected in two blocks, the first (ObsID 0658600101) from 2011 August 31 at 23:37 UT to September 1 at 12:58 UT, and the second (ObsID 0658600201) from 2011 September 1 at 20:26 to September 2 at 10:42 UT. Sgr A* was placed at the center of the XMM-Newton/EPIC field of view (away from any chip gaps). The medium filter, and full-frame mode, were used for all three EPIC instruments. The more sensitive pn camera (Strüder et al. 2001) had exposures of 41.9 and 45.2 ks in the two observations, respectively. The less sensitive MOS1 and MOS2 cameras (Turner et al. 2001) had exposures of 48.6 and 52.3 ks in the two observations. Below we focus on results from the pn camera; the MOS results were similar.

The data were processed with the XMM Science Analysis Software (version 11.0.0) to select PATTERN $\leqslant$ 12, energies between 2 and 10 keV, and FLAG = 0. We extracted light curves at 300 s binning from a radius of 10 arcsec (as typical for Sgr A*, e.g., Porquet et al. 2003), centered at the location of Sgr A*. This radius only encloses 50% of the emission from Sgr A* (Read et al. 2011), yet includes a significant amount of contamination from unrelated sources, both diffuse and point-like. In fact, the quiescent flux of Sgr A* (2–10 keV ${L}_{X}\sim 2.4\times {10}^{33}$ erg s⁻¹) is only ∼10% of the flux enclosed within $10^{\prime\prime}$ (Baganoff et al. 2003).

No statistically significant (3σ) flares were observed in any of the EPIC light curves, and the highest points in each lightcurve did not correspond with the highest points in other light curves. The most interesting possible peak occurred at 4.85 hr into the first observation, reaching 0.153 ± 0.023 counts s⁻¹, compared to an average rate of 0.10 counts s⁻¹. We can thus set an upper limit on the background subtracted Sgr A* flare luminosity during our observations (for flares of 300 s in length) of 7.6 times the quiescent value, or L_X (2–10 keV) $\lt \quad 1.8\times {10}^{34}$ erg s⁻¹; longer flares have stricter upper limits ($\lt 6\times {10}^{33}$ erg s⁻¹ on average for 1 ks flares). We assume an absorbed power-law spectrum for the X-ray flares, as seen for the quiescent Sgr A* spectrum (Baganoff et al. 2003), with photon index of 2.7 and ${N}_{{\rm{H}}}=9.8\times {10}^{22}$ cm⁻². Fits to the spectra of X-ray flares from Sgr A* span a range of spectral indices, from 1.7 to 3.2; changing the assumed spectral index in this range of photon indices alters the upper limits by 10% up or down.

In an attempt to provide overlapping 1.3 mm data, Sgr A* was also observed with the SMA. Unfortunately, it was afternoon in Hawaii during our Herschel and XMM-Newton observations. SMA observing conditions are typically worst in the afternoon because the unstable atmosphere corrupts the interferometer phases. Given the poor quality of the data we can only put a ∼30% upper limit on the amplitude of variations of Sgr A* during our observations. This corresponds to ∼1 Jy, which is about the size of the largest variations seen at 1.3 mm (Dexter et al. 2014).

We observe strong variations in all three SPIRE bands with similar amplitude in each. We do not know the absolute flux density of Sgr A* at any of these wavelengths due to confusion with the surrounding dust emission (though interferometer measurements at 0.43 mm in Marrone et al. (2006) show a minimum flux density of 2 Jy in 4 epochs). However, the negative deviation around 5 UT in the first interval implies that for all three bands there must be a minimum time-averaged flux density of at least 0.5 Jy, even at 0.25 mm where the SED is not otherwise constrained.

We now attempt to provide an empirical comparison of Sgr A* variability at SPIRE wavelengths to variability at 1.3 mm. Dexter et al. (2014) provide a detailed analysis of Sgr A* variability at 1.3, 0.8, and 0.43 mm. They demonstrated consistent variability amplitude characteristics between the bands, though the 0.8 mm and especially the 0.43 mm characteristics were poorly constrained due to the smaller number of measurements at those wavelengths. To compare our light curves to 1.3 mm observations composed of many shorter intervals of irregularly sampled data, we devised the approach described below.

Using the SPIRE light curves, as well as the ∼70 hr of 1.3 mm light curves compiled by Dexter et al. (2014), we determine most-likely variability amplitudes for overlapping 4 hr subsets of the data at each wavelength and compile them into a distribution function. The segmentation time is chosen to span typical variations in the light curves, though our results are not very sensitive to the choice. We used overlapping segments, whose start times were separated by at least 1.33 hr (66% overlap), to provide additional measurements of σ. Explicitly, for each four-hour block, we binned the data to 20 minute time resolution for better signal-to-noise, subtracted the bin mean, and then constructed the function

$$\begin{eqnarray}&&{P}_{\mathrm{block}}(\sigma | {f}_{i},{\delta }_{i})={{\rm{\Pi }}}_{i}\left(\displaystyle \frac{1}{\sqrt{2\pi }({\sigma }^{2}+{\delta }_{i}^{2})}\right){e}^{\left(\tfrac{-{f}_{i}^{2}}{2({\sigma }^{2}+{\delta }_{i}^{2})}\right)},\end{eqnarray} \tag{ 4 }$$

for the probability that σ is the typical variability amplitude within the block, given the binned mean-subtracted flux density measurements (f_i) and their uncertainties (${\delta }_{i}$). We evaluated each ${P}_{\mathrm{block}}$ at $\sigma =0$ to 2 Jy using 0.001 Jy steps. This range and step size is wide enough to capture the most probable value in each block and to finely sample the function. By taking the most likely variability amplitude from each block, we can create a single cumulative distribution function (CDF) for each band. In order to illustrate the range of possible CDFs that are consistent with our data, we numerically sample the probability function for each block 100 times and create 100 additional CDFs. Figure 5 shows the region containing the central 68 CDFs created in this way for each wavelength.

Zoom In Zoom Out Reset image size

Numerical simulations predict constant or increasing fractional variability with decreasing wavelength in the millimeter/submillimeter portion of the SED (e.g., Goldston et al. 2005; Dexter & Fragile 2013; Chan et al. 2015). While the allowed regions for the SPIRE CDFs of the variability amplitude overlap in Figure 5, the 0.25 and 0.5 mm regions are mostly distinct, and there is a clear trend of decreasing variability amplitude with wavelength, suggestive of a falling SED from 0.5 to 0.25 mm.

The distribution of 1.3 mm variability is notably different from the SPIRE curves in Figure 5. There is a long tail of high-amplitude variations not seen in the submillimeter light curves. This may be the result of catching Sgr A* during a quiet state, as past ground-based measurements of Sgr A* in the 0.45 and 0.35 mm atmospheric windows differed by several Jy (Dent et al. 1993; Serabyn et al. 1997; Pierce-Price et al. 2000; Marrone et al. 2006, 2008; Yusef-Zadeh et al. 2006). A quiet state is also suggested by the XMM-Newton light curves shown in Figures 2 and 3, which show no flare event with a 2–10 keV luminosity greater than $1.8\times {10}^{34}$ erg s⁻¹. The absence of a larger flare in our X-ray light curves is consistent with the flare rate inferred from more than 800 hr of Chandra monitoring of Sgr A* (∼1 day⁻¹ above 10³⁴ erg s⁻¹; Neilsen et al. 2013). The 1.3 mm light curve has a duration nearly three times that of the SPIRE light curve, and also extends over many years and therefore includes more of the long term fluctuations that can be expected from AGN variability.

In Figure 6 we again plot the SPIRE light curves, but now show how the 0.5–0.35 mm flux density difference changes with time. In the first interval, we see a red color following our largest observed feature, from 05:00 to 08:00 UT. Notably, the flare falls off even faster at 0.25 mm than it does at 0.35 or 0.5 mm. The reddened color as flux density decreases is suggestive of a cooling process. The expanding-blob model for features in the submillimeter light curve of Sgr A* should exhibit a blue-first fall off in flux reminiscent of what we observe, yet the simultaneous rise at all wavelengths is not consistent with the most naive blob models. It is somewhat harder to predict how color changes should manifest in more complex occultation models where both absorption and adiabatic expansion take place simultaneously (Yusef-Zadeh et al. 2010).

Zoom In Zoom Out Reset image size

During our second observing interval, the light curves exhibit a clear pattern of relatively red local maxima, and relatively blue local minima. This pattern indicates a larger absolute amplitude of variation for the longer-wavelength band. This pattern is less evident in the first interval, though there is some hint of a similar pattern in the smaller flares and the largest flare is brightest at 0.5 mm for most of its duration. This seeming change in the spectrum of the flaring emission between the first and second intervals may indicate that the quiescent 0.5–0.25 mm spectrum also changed, but without a way to directly measure the absolute flux density of Sgr A* we can only speculate.

We analyzed the power spectra of variations in the SPIRE Sgr A* light curves. To do this we combined the spectral averaging technique of Welch (1967)—to provide a high fidelity estimate of the power spectrum—and the Monte Carlo fitting approach of Uttley et al. (2002)—to properly account for the effects of our sampling pattern, aliasing, and red-noise leak on the shape of our power spectra.

Welch's method involves dividing a time series into overlapping segments, apodizing and Fourier transforming each, and then averaging the power at each frequency. We chose to use segment lengths 1/3 as long as our 12.75 hr observing intervals in order to closely match the timescale we used in our analysis in Section 4.1. We used 50% overlap as suggested by Press et al. (2002). After applying Welch's method to each of our two 12.75 hr intervals, we combined the results to produce one average power spectrum. Figure 7 shows the final power spectrum for our 0.5 mm light curve. For clarity we do not show the 0.35 mm or 0.25 mm power spectra. They are similar though, as expected based on the shape of the light curves, but less well determined.

Zoom In Zoom Out Reset image size

We modeled our light curve as arising from a power law noise process with $P(\nu )\propto {\nu }^{-\beta }$. To find the best-fitting power law slope, β, we followed Uttley et al. (2002) and simulated a large number of light curves with a given slope, sampled them according to the sampling pattern defined by our observations, computed their power spectra in the same way as for our Sgr A* light curves, and then used the distribution of simulated spectra to define a goodness of fit metric.

For a given β, we used the method of Timmer & Koenig (1995) to simulate light curves. To ensure we captured the effects of aliasing and red noise leak in our simulated spectra, we produced synthetic light curves that were 200 times as long, sampled 10 times as frequently as our observed light curves. These were then sampled at the same rate as our observations and divided into 100 pairs of 12.75 hr light curves. From each light curve, a best fit-slope was subtracted, as this was a necessary step in the reduction of our observed light curves (Section 2.1). Finally, a single power spectrum for each pair was computed using the same approach used for our observed light curves. This process is then repeated 100 times to yield 10⁴ simulated spectra.

For each of the 10⁴ simulated spectra, we computed the quantity

$$\begin{eqnarray}&&{\chi }_{\mathrm{dist}}^{2}=\displaystyle \sum _{\nu }\;\displaystyle \frac{{[{P}_{\mathrm{sim},i}(\nu )-\bar{{P}_{\mathrm{sim}}}(\nu )]}^{2}}{{\sigma }_{\mathrm{sim}}{(\nu )}^{2}},\end{eqnarray} \tag{ 5 }$$

where ${P}_{\mathrm{sim},i}$ is a single simulated spectrum, $\bar{{P}_{\mathrm{sim}}}(\nu )$ is the mean of all the simulated spectra, and ${\sigma }_{\mathrm{sim}}(\nu )$ is the standard deviation of the spectra at each frequency (Uttley et al. 2002). We then computed a similar quantity for our observed power spectrum, ${P}_{\mathrm{obs}}(\nu )$,

$$\begin{eqnarray}&&{\chi }_{\mathrm{dist},\mathrm{obs}}^{2}=\displaystyle \sum _{\nu }\;\displaystyle \frac{{[{P}_{\mathrm{obs}}(\nu )-\bar{{P}_{\mathrm{sim}}}(\nu )]}^{2}}{{\sigma }_{\mathrm{sim}}{(\nu )}^{2}},\end{eqnarray} \tag{ 6 }$$

scaling $\bar{{P}_{\mathrm{sim}}}(\nu )$ and ${\sigma }_{\mathrm{sim}}(\nu )$ by a common factor in order to minimize ${\chi }_{\mathrm{dist},\mathrm{obs}}^{2}$. Then we compared ${\chi }_{\mathrm{dist},\mathrm{obs}}^{2}$ to the distribution of ${\chi }_{\mathrm{dist}}^{2}$ defined by the 10⁴ simulated spectra. The rejection probability for a given β is taken as the percent of the simulated light curves with ${\chi }_{\mathrm{dist}}^{2}$ smaller than ${\chi }_{\mathrm{dist},\mathrm{obs}}^{2}$. We plot one minus the rejection probability versus β in Figure 8. Our best-fit power law slope is $\beta =2.40$ with a 95% confidence interval that spans from $\beta =2.16$ to $\beta =2.73$. We show our best fit model spectrum and associated errorbars in Figure 7. For comparison with our Monte Carlo-based approach, we also performed a basic fit of a line to our observed spectrum in log–log space. In this case we recover $\beta =2.25$, which is less steep than we find following Uttley et al. (2002), yet still within the 68% confidence interval.

Zoom In Zoom Out Reset image size

Our derived spectral slope is very similar to the $\beta ={2.3}_{-0.6}^{+0.8}$ measured at 1.3 mm by Dexter et al. (2014). Those authors also noted a break in the power spectrum at ${8}_{-4}^{+3}$ hr, which is a longer timescale than we have access to in our 12.75 hr intervals. Meyer et al. (2009) also found a slope of $\beta =2.1\pm 0.5$ at infrared wavelengths (mostly 2.2 μm), but with a spectral break around 2.5 hr, and Hora et al. (2014) found consistent characteristics at 4.5 μm. The consistency in slope from 1.3 mm, through the SPIRE bands, out to the IR is not unexpected, as emission at all of these wavelengths is expected to arise very close to the black hole, and therefore to be subject to the same variations in the accretion process.