THE X-RAY FLUX DISTRIBUTION OF SAGITTARIUS A* AS SEEN BY CHANDRA

After over a decade of multiwavelength monitoring of Sgr A*, the supermassive black hole at the center of the Galaxy, significant uncertainties regarding the physics of its X-ray emission persist. Observationally, the X-ray emission seen by Chandra, XMM-Newton, NuSTAR, and Swift appears to be composed of a steady quiescent background interrupted roughly daily by flares (e.g., Baganoff et al. 2001, 2003; Goldwurm et al. 2003; Porquet et al. 2003; Eckart et al. 2004; Bélanger et al. 2005; Eckart et al. 2006; Porquet et al. 2008; Trap et al. 2011; Nowak et al. 2012; Neilsen et al. 2013; Degenaar et al. 2013; Barrière et al. 2014). The quiescent flux is exceedingly faint, corresponding to a 2–10 keV luminosity L_X ≈ 3.5 × 10³³ erg s⁻¹ ≲ 10⁻¹¹L_Edd, where the Eddington luminosity L_Edd is the canonical maximum luminosity from an accreting 4.1 × 10⁶ M_☉ black hole. Presently (for a recent review, see Yuan & Narayan 2014), the quiescent X-ray source is understood as thermal plasma emission from an accretion flow extending out to the Bondi radius (Quataert 2002; Baganoff et al. 2003; Yuan et al. 2003; Liu et al. 2004; Xu et al. 2006; Wang et al. 2013). Sgr A* appears to owe its low luminosity to a combination of extremely radiatively inefficient accretion (Narayan et al. 1995; Quataert & Gruzinov 2000; Yuan et al. 2003, 2012; Blandford & Begelman 1999) and an outflow that removes ∼99% of the inflowing matter before it reaches the event horizon (see Wang et al. 2013, and references therein).

The origin of the X-ray flares is somewhat more elusive, a fact best illustrated by reference to the panoply of models suggested to produce them, which include particle acceleration and/or heating via magnetic reconnection or other stochastic processes in the inner accretion flow or base of an outflow (e.g., Markoff et al. 2001; Liu & Melia 2002; Yuan et al. 2003; Liu et al. 2004), as well as the tidal vaporization of asteroids (Čadež et al. 2008; Kostić et al. 2009; Zubovas et al. 2012). The X-rays themselves have been attributed to direct synchrotron, synchrotron self-Compton (SSC; in which infrared-emitting electrons inverse Compton scatter the infrared synchrotron photons), or other inverse Compton scenarios, in which either the IR photons are upscattered by submillimeter (sub-millimeter) emitting electrons or the sub-millimeter photons are upscattered by the IR-emitting electrons (Markoff et al. 2001; Liu & Melia 2002; Liu et al. 2004; Yuan et al. 2003, 2004, 2009; Eckart et al. 2004, 2006, 2009; Marrone et al. 2008; Dodds-Eden et al. 2009; Yusef-Zadeh et al. 2009, 2012; Witzel et al. 2012; Nowak et al. 2012; Barrière et al. 2014).

Part of the reason for the model degeneracy is the significant interstellar absorption toward the Galactic center, which makes it difficult to confirm or exclude curvature in the X-ray spectrum over the narrow bandpasses of Chandra and XMM-Newton. In addition, it has proven difficult to acquire truly simultaneous sub-millimeter, IR, and X-ray observations of a large number of flares. Presently, however, synchrotron models with cooling breaks appear to be slightly favored over inverse Compton and SSC scenarios, given the absence of significant curvature out to ∼60 keV in a NuSTAR spectrum of a bright flare (Barrière et al. 2014; see also Dodds-Eden et al. 2009; Dibi et al. 2014). But even if the model degeneracy remains, the abundance of models reveals these flares to be astrophysically significant, by far the closest example of supermassive black hole variability. With connections to the fundamental plane (and thus active galactic nuclei (AGNs), blazars, and even X-ray binaries; Merloni et al. 2003; Falcke et al. 2004; Markoff 2005; Plotkin et al. 2012), these flares may be our best opportunity to understand Sgr A* in the context of other accreting systems.

Thus an alternative approach may be required to provide deeper insight into the radiation physics of Sgr A*. For example, a multiwavelength statistical analysis that takes advantage of the wealth of available data may be a more powerful way to reveal correlations between wavelength bands and physical processes. A great deal has already been learned from statistical investigations of the Galactic center in the near infrared (NIR), where flares are at least four times more common (e.g., Genzel et al. 2003; Eckart et al. 2006). For example, Meyer et al. (2009) reported a break in the NIR power spectrum at a timescale of ∼150 minutes, which is consistent with a linear scaling between break frequency and black hole mass in relation to the AGN studied by Uttley & McHardy (2005; see also McHardy et al. 2006; Meyer et al. 2008, 2014; Do et al. 2009). Dodds-Eden et al. (2011, hereafter D11) analyzed five years of Very Large Telescope (VLT) Ks-band observations of Sgr A* and reported that the distribution of infrared flux F_IR could be described by a log-normal distribution (Equation (4)) with a median flux of ∼1 mJy, breaking to a power-law tail at high flux (≳ 5 mJy, $\sim F_{\rm IR}^{-2.7}$). Mirroring the X-ray band, they argued that this implies two states of infrared emission (quiescent and flaring). However, in a subsequent analysis of a larger VLT data set, Witzel et al. (2012, hereafter W12) found that the flux distribution was consistent with a pure power law ($\sim F_{\rm IR}^{-4.2}$). In either case, it appears that Sgr A* is continuously variable over long and short timescales in the NIR.

Thanks to the 2012 Chandra X-ray Visionary Project¹⁸ on Sgr A* (hereafter the XVP, a 3 Ms campaign to observe the Galactic center at high spectral/spatial resolution and high cadence), we have made a number of recent advances in understanding the black hole's X-ray emission, most recently demonstrating that the quiescent emission cannot be produced by a cluster of coronally active stars (Wang et al. 2013, and references therein). The Chandra data include the brightest known X-ray flare from Sgr A* (Nowak et al. 2012), as well as dozens of others (nearly tripling the number of observed flares). Based on our analysis of the luminosity and fluence distributions of these flares (Neilsen et al. 2013), we estimated that (1) detectable flares constitute approximately one-third of the total X-ray emission on timescales of 3 Ms, and (2) undetected flares likely contribute ≲ 10% of the quiescent emission. It can also be demonstrated that when extrapolated to very high luminosity (L_X ≳ 10³⁹ erg s⁻¹), our flare luminosity distribution is broadly consistent with the historical frequency of extreme outbursts of Sgr A* (e.g., Ponti et al. 2010; Clavel et al. 2013).

In addition, in the XVP data we now have a statistical sample of X-ray fluxes that approaches what is available at NIR wavelengths. We can therefore apply similar statistical techniques as a complementary approach to observationally expensive multiwavelength campaigns; by comparing statistics between wavebands, we can search for connections and new insights into the radiation mechanisms of the black hole. In this context, we have undertaken an analysis of the X-ray flux distribution of Sgr A*. In this paper, we address the X-ray light curves from a new perspective: rather than dividing the data into "quiescent" and "flare" states, we suppose that the emission in each time bin is the sum of two continuous processes, one associated with the quiescent accretion flow and one associated with the flaring, active (variable) region close to the black hole. We briefly describe the data in Section 2 and detail our analysis in Section 3. Flux distribution in hand, we present preliminary considerations of radiation mechanisms in light of our results and the NIR results in Section 4, and summarize in Section 5.

For the details of the Chandra observations and our data reduction and analysis methods, we refer the reader to Neilsen et al. (2013). Briefly, our data consist of 38 high spatial/spectral resolution Chandra HETGS observations of Sgr A* from 2012. We restrict our analysis to the 2012 data set, with its unique combination of the HETGS spectral resolution, observing cadence, and uniform calibration. From each observation, we extract 2–8 keV light curves in 300 s bins (taking detected counts from the inner 125 of the zeroth order image and the ±1st diffraction orders; see Figure 1 of Neilsen et al. 2013). The typical X-ray flare duration is ∼3 ks, so most detectable flares are well resolved by these data. However, as discussed in Nowak et al. (2012) and Neilsen et al. (2013), the X-ray light curve of Sgr A* may exhibit substructure on timescales shorter than 300 s (see also D. Haggard et al., in preparation), so the count rates reported here should be considered averages over 300 s, analogous to the (shorter) integration times in the NIR. We base our flux calculations (2–8 keV) on the results of Nowak et al. (2012): the quiescent unabsorbed flux is 0.45 × 10⁻¹² erg s⁻¹ cm⁻² and the average unabsorbed flux of the brightest flare is 21.6 × 10⁻¹² erg s⁻¹ cm⁻². We perform all our analysis in the Interactive Spectral Interpretation System (ISIS; Houck & Denicola 2000; Houck 2002).

In order to determine the NIR flux distribution of Sgr A*, W12 used Kolmogorov–Smirnov (K-S) tests to evaluate the goodness of fit for a range of power-law indices. Here, because we are interested in the fluxes in excess of the quiescent background emission from the Galactic center, we adopt a modified version of this technique. Since we incorporate every 300 s time bin from the 3 Ms XVP campaign, our present analysis samples all flares, even those too short or too faint to be identified by our detection algorithms (Neilsen et al. 2013).

3.1. Observed Flux Distribution

Following W12, we analyze a (complementary) cumulative distribution function (CDF), i.e., the fraction of count rates greater than or equal to a given rate:

$$\begin{equation} f_{\ge }(r)=\frac{1}{M}\sum _{i=1}^{M}{\tt if}(r_i\ge r), \end{equation} \tag{ 1 }$$

where r_i are the count rates and M = 9964 is the number of data points (a data point being the count rate in one 300 s time bin). X-ray count rates are preferable to flux as an observable because they are not model-dependent. The resulting empirical CDF of Sgr A* from 2012 is shown in black in Figure 1. For reference, the blue dotted line is the CDF for a Poisson distribution with a rate Q ∼ 5.24 × 10⁻³ counts s⁻¹ (see below). It is apparent that the lowest count rates (r ≲ 0.02 counts s⁻¹) are dominated by this Poisson component. Although we cannot define quiescence based on the count rate alone, this figure supports our previous estimate that ≳ 90% of the quiescent emission can be described as a Poisson process. At higher count rates, there is also a clear excess above the quiescent Poisson noise. The sharp break at the bright end of the CDF is likely a feature of the empirical CDF itself, which must converge to 1/M at the highest count rate bin.

Zoom In Zoom Out Reset image size

To study the intrinsic flux distribution, we must be able to predict count rates given an intrinsic flux. Since there is no conclusive evidence¹⁹ for significant variations in the flare spectrum (Porquet et al. 2008; Nowak et al. 2012; Neilsen et al. 2013; Degenaar et al. 2013; Barrière et al. 2014), we can use the light curve analysis from Neilsen et al. (2013) and the spectral analysis of the brightest flare from Nowak et al. (2012) to convert any unabsorbed flare flux F from our models into a model flare count rate r (i.e., above the quiescent background). Specifically,

$$\begin{equation} r(F)\approx \frac{1}{\Delta t}{\tt prand}\left[\mathcal {P}\left(F \frac{\mathcal {P}^{-1}(r_{\rm max})}{F_{\rm max}}\right)\Delta t\right], \end{equation} \tag{ 2 }$$

where Δt = 300 s is the bin time, F_max = 21.6 × 10⁻¹² erg s⁻¹ cm⁻² is the mean unabsorbed flux of the brightest flare, r_max ∼ 0.13 counts s⁻¹ is the mean observed count rate of the brightest flare, prand is a Poisson random number generator to incorporate counting noise, and $\mathcal {P}$ represents the suppression of the count rate due to photon pileup (Davis 2001; Neilsen et al. 2013). The result is approximate because we have neglected the contribution of the quiescent background to pileup. For the very low count rates involved here, pileup is still fairly linear in the count rate, and we estimate that the approximation is good to within 1%, so that any errors are at a level well below the counting noise.

3.2. Power-law Excess

To model the excess flare emission, we first consider a power law model, since a scale-free representation simplifies our investigation of the radiation mechanism (Section 4). For an excess flux F in units of 10⁻¹² erg s⁻¹ cm⁻²:

$$\begin{equation} P(F) = \left\lbrace \begin{array}{@{} ll}k F^{-\xi } & : F_1 \le F\le F_2 \\ 0 & : \mathrm{otherwise} \end{array}, \right. \end{equation} \tag{ 3 }$$

where P(F) is the probability of an excess flux F, the lower bound F₁ = 0.04 is well below our detection limit (F ∼ 0.25), the upper bound F₂ = 41 corresponds to the maximum observed count rate of ∼0.23 counts s⁻¹, and $k=(1-\xi)/(F_1^{1-\xi }-F_2^{1-\xi })$ is a normalization constant. F₁ = 0.04 leads to the best match to our data (compared to other values between 0.01 and 0.1), and should be far enough below our detection limit and the upper limit F₂ that the bounded power law model is effectively scale-free in practice. We consider two ways to test this model: statistical tests using synthetic (simulated) data and maximum likelihood methods, specifically Markov Chain Monte Carlo (MCMC).

3.2.1. Synthetic Data

To use synthetic data, we modify the technique used by W12: as they did, we generate a power law random variable according to the probability distribution in Equation (3). Here, however, we also apply Equation (2) to convert fluxes to counts, then compare the sum of the power-law process and a Poisson process (representing the quiescent emission, which is not seen in the NIR) to the observed data. Instead of a K-S test, we use a two-sample Anderson–Darling test (A-D; Anderson & Darling 1945; Pettitt 1976; Scholz & Stephens 1987), which is more sensitive to the tails of the tested distributions (a valuable property given the proportion of values at high fluxes in our data). We refer the reader to Scholz & Stephens (1987), particularly their example in Section 7, for the details of calculating the test statistic in the presence of ties, as well as a table of percentile points and significance levels p. Outside the range of this table, we use a linear extrapolation of log (p/(1 − p)), similar to the implementation in the R statistical package.²⁰

We repeat this test for a grid of 101 evenly spaced values of ξ between 1.75 and 2.2 and 101 evenly spaced values of Q_pl between 0.0048 and 0.0058 counts s⁻¹, inclusive (the subscript pl indicates quantities found in the power-law analysis). Outside these ranges, additional tests indicate that the match probabilities are negligible. For each pair (ξ, Q_pl), we generate 1000 simulated data sets including counting noise and record the average p-value from our A-D tests with the observed data. The contours of constant p-value are shown as dashed red lines in Figure 2. The highest p-values (∼0.5) are returned in the vicinity of Q_pl ∼ 5.26 × 10⁻³ counts s⁻¹ and ξ ∼ 1.93. Next, we confirm this model using a more detailed MCMC analysis.

Zoom In Zoom Out Reset image size

3.2.2. Maximum Likelihood

For the models considered here, in which the data are treated as the sum of two random processes with significant counting noise, the likelihood function is not trivial to write down. For details, see Appendix A, but briefly the probability of observing n counts is a convolution of the Poisson distribution (representing the quiescent emission) and a piled-up power law with counting noise (here representing the flare emission). The power-law term can be calculated analytically, so the log-likelihood easily calculated. This makes the problem ideal for MCMC. We have implemented the parallel stretch from emcee (Foreman-Mackey et al. 2013) in ISIS. We use an ensemble of 300 walkers—150 per free parameter—initially distributed uniformly throughout a suitably wide region of parameter space (1.5 ⩽ ξ ⩽ 2.5, $0.004\le Q_{{\rm pl}}\ {\rm {\rm counts}}^{-1}\ {\rm s}\le 0.007$). The walkers are evolved for 3000 steps, and we discard the first 600 (we estimate that the average autocorrelation time is ≲ 50 steps). In order to keep the acceptance fraction below 0.5, we set the a-parameter, which determines the size of the stretch move at each step (see Foreman-Mackey et al. 2013; Goodman & Weare 2010), to 4. We assume uniform priors for all parameters.

In Figure 2, we show the contours (black) in the ξ − Q_pl plane that contain 68%, 90%, 95%, and 99.7% of the resulting ensemble of walkers, which agrees well with the p-value contours from the A-D analysis. There is clearly a small correlation between ξ and Q_pl, but both parameters are well constrained and the correlation is smaller than in the A-D analysis. The probability density functions (PDFs) for ξ and Q_pl are presented in Figure 3. We find $\xi =1.92_{-0.02}^{+0.03}$ and Q_pl = (5.24  ±  0.08) × 10⁻³ counts s⁻¹ (these values correspond to the 16th, 50th, and 84th percentiles of the marginalized distributions²¹). To demonstrate that the ensemble is able to reproduce the observed CDF, we draw 100 random (ξ, Q_pl) pairs from the probability distribution shown in Figure 2 and use them to generate simulated Poisson-plus-power-law data (converting power law fluxes to counts using Equation (2)). The associated CDFs for the simulated data sets, shown in gray in the top panel of Figure 1, match the data very well at all count rates. Note that if we allow F₁ to vary as well, the result is F₁ = 0.04 ± 0.02 and there is no significant change in the other parameters, although the error bars are somewhat larger.

Zoom In Zoom Out Reset image size

3.3. Log-normal Excess

It is worth exploring alternatives to the power law model, and the log-normal is a natural choice. Tilted accretion disk models of Sgr A* tend to produce NIR flux distributions that are well described by a log-normal model (Dexter & Fragile 2013). A log-normal flux distribution is expected for an exponentiated damped random walk process (Meyer et al. 2014), and damped random walk models can successfully describe Sgr A*'s sub-millimeter variability (see Dexter et al. 2014, and references therein).

Ideally, we would compare the power law model to a model of the same form as used by Dodds-Eden et al. (2011), i.e., a log-normal distribution at low flux transitioning to a power law at high flux. But since the power law already provides a satisfactory description of the bright end of the X-ray flux distribution and the faint end is dominated by the quiescent Poisson process, and the more complicated function form would have an additional three free parameters, we focus on a pure log-normal process for the time being.

We proceed as in Section 3.2, replacing the power law random distribution with a log-normal probability distribution:

$$\begin{equation} P(F) = \frac{1}{\sqrt{2\pi }\sigma F}\exp \left(-\frac{(\ln {F}-\mu)^2}{2\sigma ^2}\right). \end{equation} \tag{ 4 }$$

With F scaled to 10⁻¹² erg s⁻¹ cm⁻², μ and σ are dimensionless. As explained in Dodds-Eden et al. (2011), this log-normal distribution has a median flux of exp (μ) × 10⁻¹² erg s⁻¹ cm⁻² and a multiplicative standard deviation of exp (σ).

3.3.1. Synthetic Data

For the synthetic log-normal analysis, we calculate A-D test statistics over a 101 × 101 grid covering −6 ⩽ μ ⩽ −2 (corresponding to median fluxes for the log-normal component from ∼2.5 × 10⁻¹⁵ erg s⁻¹ cm⁻² to ∼1.4 × 10⁻¹³ erg s⁻¹ cm⁻²), and 1 ⩽ σ ⩽ 5; we fix Q at Q_ln = 0.00526 counts s⁻¹, which corresponds to the peak p-value in the synthetic power-law analysis. All the μ values tested here correspond to median fluxes for the log-normal process that lie below our detection limit. This is not a weakness of the log-normal analysis, rather it is a generic property of the high flux tail: the observed excess represents only a few percent of the data bins and only a fraction of the variable process (regardless of the model; see Figure 1). Indeed, while all the simulated fluxes in both models are positive, most are so small that the most likely discrete count rate for both the power law and the log-normal components is 0 counts s⁻¹.

We present the resulting p-value contours in red in Figure 4. μ and σ are strongly correlated, with the peak match probability around μ ∼ −2.7 ≡ 6.6 × 10⁻¹⁴ erg s⁻¹ cm⁻² and σ ∼ 1.7. The p-values are somewhat lower than in the power-law analysis (this can be seen from the contours themselves, which contain very little probability above the 68% confidence line). While this result is suggestive, our likelihood analysis is better suited to determining the relative utility of the power law and log-normal models.

Zoom In Zoom Out Reset image size

3.3.2. Maximum Likelihood

The MCMC analysis for the log-normal model is very similar to that of the power-law model, although the likelihood function is somewhat more difficult to write down (Appendix B). Here, a = 3 is sufficient to keep the acceptance fraction below 0.5. Once again, we use 150 walkers per free parameter, this time initializing them to be uniformly distributed between 0.005 ⩽ Q_ln counts⁻¹ s ⩽0.007, −6 ⩽ μ ⩽ −2 and 1 ⩽ σ ⩽ 4.

As the log-normal process typically contributes somewhat less flux (see Section 3.5) than the power-law process, the quiescent count rate Q_ln is slightly higher than Q_pl:Q_ln = (5.7 ± 0.1) × 10⁻³ counts s⁻¹. The log-normal component has a median flux $1.8_{-0.6}^{+0.8}\times 10^{-14}$ erg s⁻¹ cm⁻² (μ = −4.0  ±  0.4) and a multiplicative standard deviation σ = 2.4 ± 0.2. As noted above, most of the probability mass of the log-normal distribution lies below our detection limit.

The log-normal contours and PDFs are shown in black in Figures 4 and 5. Although there was good qualitative and quantitative agreement between the MCMC and A-D results in the power-law analysis, it does appear that these methods have different sensitivities: here there is much less overlap in the A-D contours and the MCMC contours: the synthetic data analysis prefers a higher median flux and a smaller multiplicative standard deviation, perhaps to reduce the number of simulated bins brighter than the highest-flux data point. If we fix Q_ln = 0.00526 counts s⁻¹, the new MCMC contours lie between the A-D contours and the original MCMC contours, so it seems that the discrepancy is due in part to our choice to fix Q_ln in the synthetic analysis.

Zoom In Zoom Out Reset image size

As in Section 3.2.2, we draw random samples from the (Q_ln, μ, σ) ensemble, use them to generate simulated data sets, and plot the CDFs in the bottom panel of Figure 1. Again, the simulated data appear to be in fairly good agreement with the observed CDF over most of the observed count rates. The only noticeable difference with respect to the power-law process is at the high flux end, where the data and the log-normal component diverge slightly. This is not entirely unexpected: the power law has an advantage over the log-normal in matching the turnover seen in the data at high flux (the former is bounded by definition, while the latter is not). A power law not bounded at the high flux end would behave comparably to the log-normal component here. Ultimately, however, this difference is significant enough to select the power law as a formally superior model: via the Akaike information criterion, due to the difference in log-likelihoods (∼6) and the number of free parameters, the likelihood of the power-law model relative to the log-normal model is ≳ 1300. At face value, it seems that the log-normal model can effectively be ruled out for failing to perform as well as the power law at capturing the behavior of the very brightest flares from Sgr A* (but see Section 3.4).

3.4. A Note on Correlations in the Data

3.5. Undetected Flares in Quiescence

We can also use our decomposition of the flux distribution into a steady component and a variable component to address the contribution of undetected flares to the quiescent emission from the inner 125 of the Galaxy. Here, we use the same random (Q_pl, ξ) and (Q_ln, μ, σ) samples as in Figure 1 to create synthetic data sets. For each data set, we calculate the fraction of the observed counts contributed by the variable process in each simulated time bin. In Figure 6, we show the cumulative fractional contributions of the power-law and log-normal processes, i.e., the fraction of the counts at or below a given count rate that can be attributed to the variable component.

Zoom In Zoom Out Reset image size

Unsurprisingly, the low flux bins are dominated by the Poisson process. The power-law process contributes a total of ∼12% of the counts in bins with one count and less than 20% of the counts out to count rates of ∼0.02 counts s⁻¹ (corresponding to nearly 4 × the baseline quiescent count rate Q_pl). In general, the log-normal component contributes somewhat less than the power law, but the dependence on count rate is very similar. We may also note that the value of Q_pl is only ∼84% of the average count rate outside the flares (Neilsen et al. 2013), so that the remaining flux must come from flares in our analysis. This accounting (i.e., a ∼10%–15% contribution of undetected flares to the steady background emission) is consistent with the extrapolation of the observed flare distribution to low fluence and the power spectrum and the distribution of waiting times between photons in quiescence (Neilsen et al. 2013), as well as the quiescent spectrum (Wang et al. 2013) and the surface brightness distribution of Sgr A* (Shcherbakov & Baganoff 2010).

In the preceding section, we demonstrate that the X-ray flux distribution of Sgr A* is consistent with a Poisson process, representing the steady thermal emission from the quiescent accretion flow, plus a power-law process representing flare emission from within ∼tens of gravitational radii from the black hole (e.g., Barrière et al. 2014). The flaring may also be consistent with a log-normal variability process, but the power law provides a superior description. As previous analyses have found evidence of a power-law flux distribution in the NIR (D11, W12), here we make a preliminary consideration of the extent to which a statistical comparison of the X-ray and NIR flux distributions can provide insight into the X-ray radiation mechanism. We proceed with the following formalism. The flux distribution (specifically, the differential flux distribution) is defined:

$$\begin{equation} P(F_\lambda) \equiv \frac{dN}{dF_\lambda }, \end{equation} \tag{ 5 }$$

and the cumulative flux distribution (i.e., the continuous extension of Equation (1)) is given by:

$$\begin{equation} N_{\ge }(F_\lambda) \equiv \int _{F_\lambda }^{\infty }P(f_\lambda)df_\lambda, \end{equation} \tag{ 6 }$$

where F_λ is the flux at wavelength λ and dN/dF_λ is the number of times per unit flux that Sgr A* is observed with flux F_λ at wavelength λ. Given an infrared flux distribution P(F_IR), we may ask what any given radiation mechanism predicts for the X-ray flux distribution P(F_X). The IR is dominated by optically thin synchrotron emission (e.g., Hornstein et al. 2007); for a one-zone model, the flux is:

$$\begin{equation} F_\nu ^{S}\propto n_eR^3B^{1+\alpha }\nu ^{-\alpha }, \end{equation} \tag{ 7 }$$

where n_e is the IR-emitting electron density, R is the size of the emitting region, B is the magnetic field, ν is the frequency, and α is the spectral index in the optically thin regime (e.g., Rybicki & Lightman 1979). For fluxes associated with detectable X-ray emission, we suppose the NIR flux distribution can be described as a power law with index q, i.e., $P(F_{\rm IR})\propto F_{\rm IR}^{-q}.$ A parallel analysis could be performed to explore other functional forms for the flux distributions.

First, we derive a general result in the case of a power law relationship between the X-ray and IR, i.e., $F_{\rm X}\propto F_{\rm IR}^\beta,$ where β ≠ 0 is some constant. Under this condition, we have

$$\begin{eqnarray} P(F_{\rm X}) &\equiv & \frac{dN}{dF_{\rm X}} \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} &{\propto} & F_{\rm IR}^{1-\beta }\frac{dN}{dF_{\rm IR}} \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} &{\propto} & F_{\rm IR}^{1-\beta -q} \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} &{\propto} & F_{\rm X}^{-(q+\beta -1)/\beta }. \end{eqnarray} \tag{ 11 }$$

In other words, if the infrared flux distribution is a power law with index q and $F_{\rm X}\propto F_{\rm IR}^\beta$, the X-ray flux distribution will be a power law with index ξ = (q + β − 1)/β. It can also be seen that the index of the cumulative distribution in the X-ray is 1/β times that in the NIR.

In general, β is not trivial to derive from first principles, even in a one-zone model. For instance, the synchrotron flux $F_{\rm IR}^{S}$ (see Equation (7) above) depends on several physical parameters, and the response of the X-ray to any variations in F_IR will be sensitive to which of the synchrotron parameters vary, as well as their time dependence and the X-ray emission mechanism (see, e.g., Dodds-Eden et al. 2010; Yusef-Zadeh et al. 2009; Eckart et al. 2012). This sort of detailed time-dependent modeling of individual flares is beyond the scope of this paper, and is furthermore not easily applied to a statistical analysis of the flux distribution.

However, it can still be enlightening to consider the multiwavelength flux distribution in this context. For a simple, illustrative example, let us consider a scenario in which the X-ray emission is the high-energy tail of the synchrotron radiation observed in the IR. Since the synchrotron cooling time at X-ray frequencies is much shorter than the duration of the flares, radiative losses are likely very important (Kardashev et al. 1962; Yuan et al. 2003). Indeed, modeling the spectral energy distribution of a single bright flare, Dodds-Eden et al. (2009) found that synchrotron models required a cooling break (at an intermediate frequency ν_c ∼ 10¹⁵ Hz) in order to avoid violating their upper limits on the mid-IR flux. Such a break could also explain why the flare spectra may appear to be slightly steeper in the X-ray than in the NIR (see Hornstein et al. 2007; Barrière et al. 2014). In this model, assuming that ν_c is in the optical or UV, the X-ray flux is given by

$$\begin{eqnarray} F_{\rm X} &\approx & F_{\rm IR} \left(\frac{\nu _c}{\nu _{\rm IR}}\right)^{-\alpha }\left(\frac{\nu _{\rm X}}{\nu _c}\right)^{-(\alpha +1/2)} \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} &{\propto} & F_{\rm IR} \nu _c^{1/2}. \end{eqnarray} \tag{ 13 }$$

If ν_c has a power-law dependence on the IR flux, then there is some β for which $\nu _c\sim F_{\rm IR}^{2\beta -2}$, so that $F_{\rm X}\, {\propto}\, F_{\rm IR}^\beta$ and Equation (11) applies.

Observationally, since ξ ≲ q, we use Equation (11) to infer β ≳ 1. For the specific case of the cooling break in Equation (13), β ≳ 1 implies that ν_c is either independent of or increases with F_IR. This is counterintuitive because ν_c and the synchrotron flux (Equation (7)) decrease and increase, respectively, with both B and R (for ν_c; see, e.g., Equation (3) in Dodds-Eden et al. 2009). However, β ≳ 1 can be understood in the context of more sophisticated models of the time-dependent emission. For instance, Dodds-Eden et al. (2010) modeled the light curve of a single flare with a temporary increase in n_e, anticorrelated with changes in the magnetic field (producing both an increase in flux and β > 1). Such correlations might be expected for magnetic reconnection events, and can be accounted for in self-consistent simulations of the plasma parameters during flares (e.g., Dibi et al. 2014; S. Dibi et al., in preparation).

It is less straightforward to apply similar analytical considerations to the other proposed X-ray emission mechanisms (i.e., SSC and external Comptonization scenarios). Consider the SSC case. In the Thomson scattering limit for a homogeneous sphere of low optical depth (see Bloom & Marscher 1996; Marrone et al. 2008, and references therein), the SSC flux from first-order scattering is proportional to the synchrotron flux times the Thomson optical depth, which is proportional to the product n_eR:

$$\begin{equation} F_\nu ^{\rm SSC}\propto n_e^2R^4B^{1+\alpha }\nu ^{-\alpha }. \end{equation} \tag{ 14 }$$

If the IR variability could be attributed to power-law variations in either n_e or R in this scenario, one could write down β and predict the X-ray flux distribution. However, we have already established that it is likely that flares involve possibly correlated variations of several parameters. Furthermore, in the model of Dodds-Eden et al. (2009), the SSC flux peaked in or near the X-ray, so that the assumptions used to derive Equation (14) may not apply. A complete treatment, beyond the scope of this paper, would require integrating over the distributions of seed photons and electrons and allowing for variations in physical parameters.

The situation for external Comptonization scenarios is similar. If the X-ray emission is produced by IR photons scattering off sub-millimeter-emitting electrons, we might expect the X-ray and IR to be related more or less linearly. If instead the sub-millimeter photons scatter off IR-emitting electrons, there need not be any correlation between the X-ray and infrared. However, the accuracy of these conclusions depends on the similarity of the plasma properties and evolution at sub-millimeter and IR wavelengths.

Perhaps it should come as no surprise that in the 15 yr since Baganoff et al. (2001, 2003) discovered the incredibly faint counterpart of Sgr A* with Chandra, it has proven difficult to pin down the underlying radiation physics of this elusive X-ray source. To date, there is still no consensus as to the X-ray emission mechanism or the underlying physical process(es) responsible for the variability. However, there appears to be some grounds for optimism, as the statistical behavior of the X-ray emission appears to mirror that of the NIR, at least qualitatively. For instance, at both wavelengths, the source is continuously "on" and its flux distribution is partially or fully consistent with a power law distribution (e.g., Dodds-Eden et al. 2011; Witzel et al. 2012; Meyer et al. 2014, and references therein for the NIR).

Although no clear correlation has been found between the peak flux of flares observed simultaneously in NIR and X-rays (Trap et al. 2011), the flare peaks are typically simultaneous within observational uncertainties (see Yusef-Zadeh et al. 2012, and references therein). Furthermore, flares at both wavelengths exhibit significant fast variations on short timescales (NIR: 47 s, Dodds-Eden et al. 2009; X-ray: ∼100 s, Nowak et al. 2012; Barrière et al. 2014), indicating that excursions to high flux come from compact regions (within ≈10 Schwarzschild radii of the black hole in the case of the X-ray flares; Barrière et al. 2014). The emission at different wavelengths need not be produced by the same electrons (e.g., Hornstein et al. 2007), but it is plausible that these variable processes share a common physical origin.

Thus, the excellent data now available present a new and exciting opportunity to gain insight into the origin of the multiwavelength variability of Sgr A*. In this paper, we have explored the total X-ray flux distribution from the inner 125 of the Galaxy using statistical methods previously applied to Sgr A*'s NIR emission. Rather than focusing on individual flares, we approach the radiation physics of Sgr A* from a different perspective, looking for connections between emission at different wavebands once we understand the statistical behavior of each waveband separately.

As is evident from Figure 1, the X-ray emission is dominated²² by the faint steady quiescent source and diffuse background. Our analysis of the flux distribution confirms the steadiness of the quiescent emission: roughly 85%–90% of it is consistent with a Poisson process, i.e., it exhibits no detectable intrinsic variability at all (see also Neilsen et al. 2013). To the extent that this thermal plasma emission from large scales (e.g., Wang et al. 2013) is constant in time, all the X-ray variability from Sgr A* can be attributed to non-thermal flaring.

The variable source, likely located in the inner accretion flow, is typically an order of magnitude fainter than the thermal plasma, but our statistical analysis indicates that it can be described by a single variability process over a dynamic range of three orders of magnitude in flux. At the faint end (below our flare detection limit), the variable source contributes the remaining ∼10%–15% of the flux during quiescent/steady intervals; at the bright end, flares can exceed the steady background by factors of 100 or more (Nowak et al. 2012; Porquet et al. 2003, 2008). Both a power-law process (Equation (3)) with index $\xi =1.92_{-0.02}^{+0.03}$ and a log-normal process (Equation (4)) with location parameter μ = −4.0 ± 0.4 and shape parameter σ = 2.4 ± 0.2 are consistent with the observed excess above the quiescent emission, but the power law provides a significantly better description of the data. Notably, ξ appears to be the same as the power law index of the distribution of flare luminosities ($1.9^{+0.3}_{-0.4}$; Neilsen et al. 2013). This result makes sense if the variable process is the sum of many discrete flares, but more work is required to know if this is a unique interpretation.

Our results place a strong observational constraint on theoretical models of the multiwavelength variability of Sgr A*. In addition to the flare SED and the variable X-ray/NIR ratio, any viable model must be able to reproduce the distribution of fluxes in both the NIR and the X-ray independently. Even the one-zone models considered in Section 4 can begin to inform our understanding of the relationship between the NIR and X-ray emission, but it is instructive to note why (in addition to the limited power of one-zone models) a complete answer may require more than simple analytical considerations. The primary reasons are physical:

• 1.  
  Even given an X-ray emission mechanism, the X-ray variability can only be predicted analytically from the NIR variability if the origin of the NIR variability is known.
• 2.  
  The NIR power law is probably not attributable to a single parameter with a power law distribution. Instead, the synchrotron flux likely varies due to non-power-law variations of multiple parameters (as in the time-dependent analysis of Dodds-Eden et al. 2010).

Self-consistent calculations of the particle distribution during flares (Dibi et al. 2014) can overcome these difficulties, and we will analyze the multiwavelength flux distribution in the context of these calculations in future work (S. Dibi et al., in preparation). It may also be possible to understand the variability in light of other simulations of the accretion flow (e.g., Mościbrodzka & Falcke 2013).

In addition, several factors complicate multiwavelength statistical analyses of Sgr A*:

• 1.  
  There is disagreement in the literature over the functional form of the NIR flux distribution, which leads to different estimates of the index q of the power-law portion (i.e., q = 2.7, D11; q = 4.2, W12). This is particularly relevant in the context of Equation (11).
• 2.  
  Any deviations from a pure power-law distribution break the scale-free nature of our calculations, so that in order to validate the comparison performed in Section 4, we would need to know which NIR fluxes correspond to which X-ray fluxes. This may not be problematic, as deviations from a power law (Dodds-Eden et al. 2011) have only been reported at NIR fluxes that are not typically associated with detectable X-ray flares (≲ 5 mJy; see also Marrone et al. 2008).
• 3.  
  The flux distribution may not be completely stationary on timescales of individual observations. W12 argue that the two-state flux distribution of D11 (i.e., lognormal at low flux, power law at high flux) is biased by a single bright flare. Such behavior might be expected if the power law is an aggregate effect of isolated flares in distinct regions of the accretion flow. Caution is therefore merited when comparing X-ray and NIR flux measurements taken over different time intervals and at different epochs (particularly because the X-ray flux distribution is consistent with several functional forms).
• 4.  
  A related issue is that the characteristic timescales of the X-ray variability process are unknown. Although the near simultaneity of X-ray and NIR flares suggests that the two are fundamentally related, an empirical confirmation of similar variability timescales would lend significant credence to any multiwavelength statistical analysis of the radiation mechanisms at work in Sgr A*. We intend to perform this computationally demanding measurement in the near future.

In short, there is still more work to be done to understand the shape of the flux distributions at different wavelengths and their relationships to each other. These points highlight the need for continued coordinated multiwavelength campaigns on Sgr A* (specifically, simultaneous observations). By tracking simultaneous variations in the X-ray and IR, we can make a direct measurement of the IR/X-ray flux scaling, free of any lingering concerns about observations taken at different epochs. Strictly simultaneous observations would also break the degeneracy related to the functional form of the NIR/X-ray flux distributions, since the same model or parameterization can be applied to both wavebands. Although complicated somewhat by frequency-dependent lags (e.g., Marrone et al. 2008; Falcke et al. 2009, and references therein), the addition of radio/sub-millimeter variability statistical analyses (e.g., Dexter et al. 2014) will paint a more comprehensive picture of variability in the Galactic center.

Although we must leave the question of the primary X-ray radiation mechanism of Sgr A* for future work, a statistical approach is still promising. Our clear and precise measurement of the X-ray flux distribution is a critical first step, and our analysis provides a new and strong observational constraint for theoretical models of Sgr A*. With the advances described above, the prospects for revealing the physics of variable X-ray emission in Sgr A* in the near future remain strong.

We thank all the members of the Sgr A* Chandra XVP collaboration, and we are immensely grateful to Chandra Mission Planning for their support during our 2012 campaign. We thank Daniel Wang, Salomé Dibi, Feng Yuan, and Stefan Gillessen for very useful discussions concerning the analysis and interpretation of our data. We acknowledge the role of the Lorentz Center, Leiden, and the Netherlands Organization for Scientific Research Vidi Fellowship 639.042.711 (S.M.). J.N. gratefully acknowledges funding support from NASA through Hubble Postdoctoral Fellowship grant number HST-HF2-51343.001-A and through the Einstein Postdoctoral Fellowship, grant PF2-130097, awarded by the Chandra X-Ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060, and from NASA through the Smithsonian Astrophysical Observatory contract SV3-73016 to MIT for support of the Chandra X-Ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060. F.K.B. acknowledges support from the Chandra grant G02-13110A under contract NAS8-03060. N.D. is supported by NASA through Hubble Postdoctoral Fellowship grant number HST-HF-51287.01-A.

In order to use maximum likelihood methods, we naturally must be able to write down the likelihood function. In general, for the sum of two random variables X and Y, the probability of observing a value of n given parameters θ_X and θ_Y is a convolution:

$$\begin{equation} \mathcal {P}(n|\theta _X,\theta _Y)=\sum _{j=0}^{n}\mathcal {P}_X(j|\theta _X)\mathcal {P}_{Y}(n-j|\theta _Y), \end{equation} \tag{ A1 }$$

and for a data set D with M measurements n_i, the likelihood is:

$$\begin{equation} \mathcal {P}(D | \theta _X,\theta _Y)=\prod _{i=1}^{M}\mathcal {P}(n_i|\theta _X,\theta _Y). \end{equation} \tag{ A2 }$$

Below, we calculate the likelihood function using the models described in the main text.

As described in Section 3.2, we first model the light curve as the sum of a Poisson process and a power-law process. The Poisson process has parameter Λ ≡ (300 s)Q, where Q is the quiescent count rate, and the well-known probability of observing n counts is:

$$\begin{equation} \mathcal {P}_p(n|Q)=\frac{\Lambda ^{n}e^{-\Lambda }}{n!}. \end{equation} \tag{ A3 }$$

The power-law probability is given by Equation (3), with P(F)∝F^−ξ on the interval [F₁, F₂]. In order to use this power law in Equation (A1), we must convert to counts, include the effects of pileup, and add counting noise. The latter two requirements make this a non-trivial likelihood function to write down.

First, we set out (and remind the reader of) some variables. As before, n shall be the number of observed counts and F shall be a power-law flux. We shall use c to denote the counts implied by applying the flux/counts scaling from the brightest flare (Section 3; see also Nowak et al. 2012) to F, and m shall be the result of applying pileup effects to c. In this scheme, F, c, and m are continuous variables, while n is discrete. To add counting noise, we need to find the probability of observing n counts given m counts in the piled-up power-law process and then integrate over all values of m:

$$\begin{equation} \mathcal {P}_{{\rm pl}}(n|\xi)=\int _{m_1}^{m_2}P(n|m)P(m|\xi)dm. \end{equation} \tag{ A4 }$$

Here m₁ and m₂ are the upper and lower limits of the integral calculated by scaling the bounds for the power law F₁ and F₂ to counts and then applying our pileup model. The complication is that $\mathcal {P}(m|\xi)$ is not a power law. Based on our analysis of pileup in these data (Nowak et al. 2012; Neilsen et al. 2013), we have determined that c(m) can be approximated to within 0.5% over the observed range of count rates as:

$$\begin{equation} c(m) = A_1 m^{A_2} e^{A_3 m}, \end{equation} \tag{ A5 }$$

with A₁ = 0.9761415, A₂ = 1.004234, and A₃ = 0.002401043. Then

$$\begin{eqnarray} \mathcal {P}_{{\rm pl}}(m|\xi) &=& \frac{dN}{dm} \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} &=& \frac{dN}{dc}\frac{dc}{dm} \end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray} &=&K_{{\rm pl}}(A_1 m^{A_2} e^{A_3 m})^{1-\xi }(A_3 + A_2/m). \end{eqnarray} \tag{ A8 }$$

The normalization constant K_pl is easily found to be:

$$\begin{equation} K_{{\rm pl}} = \frac{(\xi -1)}{\left(A_1 m_1^{A_2} e^{A_3 m_1} \right)^{1-\xi }- \left(A_1 m_2^{A_2} e^{A_3 m_2} \right)^{1-\xi }}. \end{equation} \tag{ A9 }$$

Finally, we can write the probability of observing n counts in the power law as

$$\begin{eqnarray} \mathcal {P}_{{\rm pl}}(n|\xi) &\propto & K_{{\rm pl}} \int _{m_1}^{m_2}\frac{m^n e^{-m}}{n!}(A_1 m^{A_2} e^{A_3 m})^{1-\xi }(A_3 + A_2/m) dm \nonumber\\ \end{eqnarray} \tag{ A10 }$$

$$\begin{eqnarray} &\equiv &\frac{K_{{\rm pl}} A_1^{1-\xi }}{n!}\int _{m_1}^{m_2}m^{A-1} e^{-B m}(A_3 + A_2/m) dm \end{eqnarray} \tag{ A11 }$$

$$\begin{eqnarray} &=&\frac{K_{{\rm pl}} A_1^{1-\xi }}{n!} [A_3 B^{-A}[\Gamma (A,m_1B)-\Gamma (A,m_2 B)] \nonumber\\ &&+ A_2 B^{-A+1}[\Gamma (A-1,m_1B)-\Gamma (A-1,m_2 B)]] \qquad \end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray} &=&\frac{K_{{\rm pl}} A_1^{1-\xi }B^{-A}}{n!} [A_3 [\Gamma (A,m_1B)-\Gamma (A,m_2 B)] \nonumber\\ &&+A_2 B[\Gamma (A-1,m_1B)-\Gamma (A-1,m_2 B)]]. \end{eqnarray} \tag{ A13 }$$

For simplicity, we have made the substitutions A ≡ n + A₂(1 − ξ) + 1 and B = 1 − A₃(1 − ξ). Γ is the upper incomplete gamma function. It turns out that Equation (A13) itself is not normalized, so we need divide by a sum over all plausible values of n,

$$\begin{equation} \mathcal {P}_{{\rm pl}}(n|\xi)\leftarrow \frac{\mathcal {P}_{{\rm pl}}(n|\xi)}{\sum ^{i=150}_{i=0}\mathcal {P}_{{\rm pl}}(i|\xi)}. \end{equation} \tag{ A14 }$$

To test this formula, we generated random numbers to represent ξ, F₁, and F₂ (4.34, 0.76, and 88.24, respectively), drew 10 million random numbers from the associated power-law distribution, and incorporated pileup and counting noise. The distribution of the resulting counts is shown in black in Figure 7, and is very well matched by Equation (A14) (red). With confidence in the power law probability, we can easily use Equations (A1) and (A2) to calculate the logarithm of the likelihood given any particular set of parameters.

Zoom In Zoom Out Reset image size

For the log-normal model, the situation is very similar, except for the obvious substitution in Equation (A8). Here, with the definition κ = F/c, we have:

$$\begin{eqnarray} \mathcal {P}_{{\rm ln}}(m|\mu,\sigma)&=&\frac{dN}{dF}\frac{dF}{dc}\frac{dc}{dm} \end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray} &=&\frac{K_{{\rm ln}}}{\sqrt{2\pi }\sigma F}\exp \left[-\frac{(\ln {F}-\mu)^2}{2\sigma ^2}\right]\kappa A_1 m^{A_2}e^{A_3 m}\left(A_3+\frac{A_2}{m}\right) \nonumber\\ \end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray} &=&\frac{K_{{\rm ln}}}{\sqrt{2\pi }\sigma }\exp \left[-\frac{(\ln {(\kappa A_1 m^{A_2}e^{A_3 m})}-\mu)^2}{2\sigma ^2}\right]\left(A_3+\frac{A_2}{m}\right). \nonumber\\ \end{eqnarray} \tag{ B3 }$$

The normalization K_ln can be found for a finite interval [m₁, m₂]:

$$\begin{eqnarray} && 2\left[{\rm erf}\left(\frac{\ln {(\kappa A_1 m_2^{A_2}e^{A_3 m_2})}-\mu }{\sqrt{2}\sigma }\right) \right. \nonumber\\ &&\quad - \left.{\rm erf}\left(\frac{\ln {(\kappa A_1 m_1^{A_2}e^{A_3 m_1})}-\mu }{\sqrt{2}\sigma }\right)\right]^{-1}. \end{eqnarray} \tag{ B4 }$$

In the limit a → 0, b → ∞, K_ln ≈ 1. Unfortunately, the probability of observing n counts in the log-normal component, i.e.,

$$\begin{eqnarray} && \mathcal {P}_{{\rm ln}}(n|\mu,\sigma) \propto \int _{m_1}^{m_2}\frac{m^n e^{-m}}{n!}\frac{K_{{\rm ln}}}{\sqrt{2\pi }\sigma } \nonumber\\ && \times \exp \left[-\frac{(\ln {(\kappa A_1 m^{A_2}e^{A_3 m})}-\mu)^2}{2\sigma ^2}\right]\left(A_3+\frac{A_2}{m}\right) dm, \qquad \end{eqnarray} \tag{ B5 }$$

cannot be calculated analytically. Instead, we estimate $\mathcal {P}_{{\rm ln}}(n|\mu,\sigma)$ empirically by generating 2.5 million log-normal random fluxes for each (μ, σ) pair, converting to piled-up counts, and adding counting noise. It is only necessary to estimate this quantity for n between 0 and 69. The rest of the likelihood calculation proceeds as in Appendix A.