THE EFFECTS OF MAGNETIC FIELDS AND INHOMOGENEITIES ON ACCRETION DISK SPECTRA AND POLARIZATION

The most well-understood accretion flow state onto a black hole is that of a geometrically thin, optically thick disk. Considerable theoretical effort has been devoted to calculating spectral models of such disks for comparison with observations (e.g., Kolykhalov & Sunyaev 1984; Laor & Netzer 1989; Shimura & Takahara 1995; Sincell & Krolik 1998; Hubeny et al. 2001; Davis et al. 2005; Hui et al. 2005; Davis & Hubeny 2006). These models have been compared with observed colors and spectra of active galactic nuclei (e.g., Laor 1990; Bonning et al. 2007; Davis et al. 2007), ultraluminous X-ray sources (Hui & Krolik 2008), and black hole X-ray binaries (Davis et al. 2006; Done & Davis 2008). In the last case, where the models appear to perform quite well for the high/soft state, attempts have been made to use them to measure the spins of black holes from observed X-ray continuum spectra (Shafee et al. 2006; Davis et al. 2006; McClintock et al. 2006; Middleton et al. 2006; Liu et al. 2008; Gou et al. 2009).

Note, however, that all such models are based on certain ad hoc assumptions. First, they generally adopt radial profiles of surface mass density and emissivity based on some form of the alpha stress prescription of Shakura & Sunyaev (1973) and some inner boundary condition on the stress (usually zero). These models also make various assumptions about the vertical structure of the disk at each radius. In particular, at a given radius the structure depends only on height above the midplane, and inhomogeneities generated by instabilities and turbulence are neglected. The disk is usually assumed to be vertically supported against the tidal gravity of the black hole by gas and radiation pressure alone. Some prescription for the vertical distribution of dissipation is adopted, e.g., that the dissipation rate per unit mass is constant. Finally, some mechanism for vertical heat transport is assumed, e.g., radiative diffusion and/or convection.

Significant progress in understanding the vertical structure of accretion disks has been made with vertically stratified shearing box simulations of magnetorotational turbulence (Balbus & Hawley 1998). In particular, radiation magnetohydrodynamical simulations of the vertical structure have now been done for a broad range of midplane radiation to gas pressure ratios of relevance to black hole accretion disks (Turner 2004; Hirose et al. 2006, 2009; Krolik et al. 2007; Blaes et al. 2007). These simulations have demonstrated that magnetic forces contribute significantly to the vertical hydrostatic balance of the disk, and can be dominant over gas and radiation pressure gradient forces in the surface layers where the spectrum is formed. At the same time, these magnetically dominated regions are Parker unstable, and as a result, very large density inhomogeneities form near the thermalization and scattering photospheres. The mass density falls off steeply with height above the midplane, and most of the accretion power is (numerically) dissipated at high optical depth near the midplane regions. Radiative diffusion completely dominates the vertical heat transport, and the fraction of accretion power carried out to the photospheres by Poynting or mechanical fluxes is negligible. Finally, the overall turbulent stress scales approximately with the total thermal pressure near the midplane. Nonetheless, the inner radiation pressure-dominated regions of the disk are thermally stable (Turner 2004; Hirose et al. 2009). They may, however, be subject to inflow instability (Lightman & Eardley 1974).

In previous work (Davis et al. 2005; Blaes et al. 2006), we incorporated the numerical dissipation profiles computed in the simulations of Turner (2004) and Hirose et al. (2006) in one-dimensional non-LTE vertical structure models of disk annuli. For models that neglect magnetic support, we found very little difference in the emergent spectra compared with models based on the assumption of constant dissipation per unit mass, the common prescription that gives constant density in radiation-dominated regions. Even though the numerical dissipation profiles have the dissipation per unit mass, generally increasing outward near the photosphere, most of the accretion power is still dissipated and thermalized at high optical depth, resulting in very similar spectra.

In contrast, the addition of substantial magnetic support of the outer layers of the disk results in significant differences in the emergent spectra (Blaes et al. 2006). Compared with models that neglect magnetic forces (e.g., Davis et al. 2005), magnetically supported atmospheres have larger density scale heights, resulting in smaller densities at the photosphere. This increases non-LTE effects and the ionization state of the gas, and also increases the ratio of scattering to absorption opacity. All of these effects result in a hardening of the disk spectrum (Blaes et al. 2006; Begelman & Pringle 2007).

However, the strong density inhomogeneities that are seen in the surface layers may act to soften the spectrum, both by increasing the thermalization of the radiation through the enhanced absorption opacity in the denser regions and by enhancing the rate of cooling as photons diffuse faster through the low density regions (Davis et al. 2004; Begelman 2006).

The fact that accretion disks in black hole X-ray binaries in the high/soft state are expected to be electron scattering dominated implies that their thermal spectra should be significantly polarized (Chandrasekhar 1960). However, the complex photospheric geometry that results from the strong density inhomogeneities seen in the simulations may reduce the degree of polarization (Coleman & Shields 1990). Even if the disk surface is flat, relativistic effects can dilute the polarization (Stark & Connors 1977). On the other hand, relativistic effects can also cause photons emitted by one part of the disk to scatter off a different part of the disk, enhancing the polarization and rotating its angle (Agol & Krolik 2000; Schnittman & Krolik 2009). Strong photospheric magnetic fields can also Faraday depolarize the emergent radiation field (Gnedin & Silant'ev 1978), and the fields seen in the simulations are marginally strong enough to do just that (Blaes et al. 2007). Future X-ray polarimetry satellites may therefore provide valuable diagnostics of surface magnetic fields in black hole accretion disks (Gnedin et al. 2006; Silant'ev 2007), providing a good observational test of the simulation results.

In this paper, we present the results of spectropolarimetric radiative transfer calculations through representative simulation domains that cover the range of radiation to gas pressure regimes that are relevant to black hole X-ray binaries. We quantify the degree of spectral hardening due to magnetic support and the degree of spectral softening due to density inhomogeneities. We also calculate the expected polarization signatures of both these effects and confirm that X-ray polarimetry will indeed be able to provide interesting diagnostics on black hole accretion disk magnetic fields.

This paper is organized as follows. In Section 2, we briefly describe the numerical simulation data that we use in our radiative transfer calculations. In Section 3, we compute one-dimensional non-LTE atmosphere models based on horizontally averaged dissipation and magnetic force profiles taken from the simulations. In Section 4, we incorporate the fully three-dimensional structure of the simulations by using Monte Carlo radiative transfer calculations to compute the local emergent spectrum. We keep track of the polarization in these calculations and discuss the emergent polarization spectrum in Section 5. Using certain scalings that we derive from comparing the different simulations, we present an illustrative full disk spectrum of the polarization in Section 6. In Section 7, we discuss our results, in particular quantifying the uncertainties in spectral color correction factors and discussing how X-ray polarimetry could be used as a diagnostic of disk magnetic fields. We summarize our conclusions in Section 8. Our Monte Carlo calculations fully incorporate polarized Compton scattering, and we summarize the relevant equations and methods in the Appendix. Readers interested mainly in observational applications of our work may begin by reading Sections 7 and 8.

We use data from three vertically stratified shearing box simulations in this paper: a simulation of a gas pressure dominated box (0326c, Hirose et al. 2006), a box with comparable midplane gas and radiation pressures (0528a, Krolik et al. 2007; Blaes et al. 2007), and a radiation pressure-dominated simulation (1112a, Hirose et al. 2009). The reader should consult these papers for details of these simulations. We summarize the most relevant parameters for spectral modeling in Table 1. In particular, M is the mass of the black hole, assumed to be non-spinning; r is the radius in the disk where the box would be located, in units of the gravitational radius r_G ≡ GM/c²; and Ω is the local angular velocity of the flow, used to approximate the local tidal gravitational acceleration through |g| = Ω²|z|, where z is the height above the midplane. Because we are interested in the full three-dimensional structure of the medium, we do not use time averages but instead select three representative epochs in each simulation. Table 1 also shows the column mass density m₀ to the midplane and the effective temperature T_eff of the local emergent flux at that epoch as computed by the simulation.

Figure 1 shows the distribution of column mass times the dissipation per unit mass as a function of column mass m for each of the three simulation data sets. We measure this in two ways. For the upper panel, we difference the vertical radiative flux with respect to column mass, and for the lower panel we show the actual numerical dissipation rate per unit mass . If local radiative equilibrium held exactly, then these two quantities should have identical distributions. This is not exactly true in the simulations, because a small fraction of the dissipated heat is transported vertically outward through the advection of internal energy (Hirose et al. 2009). However, this is most important near the midplane, and is completely negligible in the low column mass density regions where the emergent spectrum is formed.

Zoom In Zoom Out Reset image size

We have scaled the data shown in Figure 1 with the midplane column mass density m₀ and the emergent radiative flux F₀ for each simulation. With this normalization, it is remarkable how similar all the distributions are, despite the different physical conditions in each simulation and the different numerical resolutions in the vertical direction. Within the noise, all the dissipation profiles are well fit by a power law of the form

$$\begin{eqnarray} &&m\epsilon =-m\frac{dF}{dm}=\frac{F_0}{2}\left(\frac{m}{m_0}\right)^{1/2}. \end{eqnarray} \tag{ 1 }$$

Note that this differs slightly from the broken power-law profile used by Blaes et al. (2006, Equation (1)), which was based on a fit to the dissipation profile computed from the 0326c simulation, but is equivalent to the scaling derived by Krolik et al. (2007) for the 0528a simulation and used for the initial condition of the radiation dominated run 1112a (Hirose et al. 2009). Blaes et al. (2006) found mdF/dm ∝ m^0.4 at high column mass densities and ∝m at low column mass densities. Since this is only an approximate relation and subject to differences in averaging of the simulation profiles, the small difference in the exponents is not significant. We adopt the 1/2 exponent because it is the best match to the high column density behavior in all three simulation, to within the uncertainties. It is important to recognize that all these dissipation profiles have most of the dissipation occurring at high column mass densities, much deeper than the spectral forming regions of the annuli.

Figure 2 shows the horizontally averaged profiles of magnetic acceleration a_mag, scaled with the local gravitational acceleration g, as a function of scaled column mass density from the three simulation data sets. The magnetic acceleration is the net result of magnetic pressure and tension forces, which are comparable in magnitude in the low column mass density regions. Near the surface, the magnetic acceleration can be larger than the gravitational acceleration (even after averaging) due to deviations from hydrostatic equilibrium, and even be negative (inward).

Zoom In Zoom Out Reset image size

Such excursions may give rise to interesting variability, but we are primarily interested in the time-averaged behavior. In an attempt to enforce hydrostatic equilibrium better, we rescale the magnetic acceleration via

$$\begin{eqnarray} &&a^{\prime }_{\rm mag}=\frac{a_{\rm mag}}{a_{\rm tot}}g, \end{eqnarray} \tag{ 2 }$$

where a_tot is the total acceleration due to radiation, magnetic fields, and gas pressure gradients in the simulation. This rescaling is also shown in Figure 2. It significantly improves the agreement between the vertical density profiles of the simulation data and those produced in our one-dimensional atmosphere models to which we now turn.

We use the TLUSTY code (Hubeny & Lanz 1995) to compute one-dimensional vertical structure models and emergent spectra from annuli whose input parameters (m₀, T_eff, and Ω) correspond to the three snapshots summarized in Table 1. We also incorporate the vertical dissipation and magnetic support profiles measured from the three simulation epochs.

The effects of the vertical support and dissipation profile are qualitatively similar in all three cases, so we only plot the model corresponding to 1112a in Figure 3. (A similar plot for the 0326c snapshot can be found in Figure 6 of Blaes et al. 2006.)

Zoom In Zoom Out Reset image size

The solid curve shows the "standard" annulus with constant dissipation rate per unit mass (i.e., constant dF/dm) and no magnetic support. Since the annulus is radiation pressure dominated throughout most of its extent, the density is very nearly constant, as is commonly assumed. However, the effective optical depth of unity (marked by the squares) shows that the spectral forming layer occurs close to the surface where gas pressure gradients dominate and balance gravity to ensure hydrostatic equilibrium. As discussed in previous work (Davis et al. 2005, 2006; Davis & Hubeny 2006; Blaes et al. 2006; Done & Davis 2008), this means that the spectral properties are sensitive to the gas pressure scale height near the surface, but rather insensitive to midplane properties.

A comparison of Figures 3 and 4 demonstrates this lack of sensitivity. The dashed curve corresponds to a model that includes the modified dissipation profile, but still neglects magnetic support. Figure 3 shows a significant change in the density profile, which is no longer constant at large depths and yields a higher central density. Also, the temperature is generally lower at a given height beneath the effective photosphere relative to the standard model. However, the observed spectra in Figure 4 remain very similar with only a small amount of hardening for the model with enhanced surface dissipation. This can be understood by noting that both the temperature and density at the effective photosphere (denoted by the squares) occur near the surface where the density declines rapidly. Because the surface pressure scale heights are still very similar for the two models, the spectra are relatively unchanged even though the midplane properties differ. This is generically the case as long as only a small fraction (≲10%) of the dissipation occurs above the effective photosphere (see Done & Davis 2008, for further discussion). All the three snapshot models obey this constraint, and inspection of Figures 4–6 shows that modified dissipation implied by the simulations has only a small effect on the spectrum.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The model including magnetic support is shown as the dot-dashed curve in Figure 3. The addition of magnetic support greatly increases the pressure scale height near the surface. The resulting profiles are in qualitative agreement with the horizontally averaged profiles taken directly from the simulation (dotted curves), although they are not as smooth as the simulation results. Also, the density is not strictly monotonic decreasing as required for stable hydrostatic equilibrium. These inconsistencies result from the assumption of hydrostatic equilibrium in the atmosphere calculation, even though the simulations can have significant, non-hydrostatic dynamics near the surface. Also, the temperature differences arise partly due to differences in radiative equilibrium in the two cases. The radiation in the simulations is handled through flux-limited diffusion (Levermore & Pomraning 1981) and gray opacity which only accounts for free–free emission and absorption. The one-dimensional calculation solves full radiative transfer with bound-free opacities, but neglects the effects of inhomogeneities in the full three-dimensional geometry of the simulations. This point is expanded in Section 4.2.

Due to the much larger pressure scale height at the surface, the density (temperature) at the effective photosphere is much lower (higher). As discussed by Blaes et al. (2006), this combination of higher temperature and lower density leads to harder spectra, due to both the effects on absorption features at lower temperatures and the effects of electron scattering at higher temperatures. Indeed, Figures 4–6 show that magnetically supported annuli all give harder spectra (dot-dashed curves), demonstrating that the qualitative result (harder spectra) is quite general. We can make this more quantitative by comparing the color-corrected blackbody models that best match the magnetized and standard atmospheres. This estimate is most robust for the model corresponding to the 1112a snapshot since the strong absorption edges in the other models prevent the color-corrected blackbody from providing a good fit, making the color correction (or hardening factor) somewhat ambiguous. For the 1112a model, which lacks significant edge features, we find the color correction increases from about 1.67 to 1.84 when incorporating magnetic support, a roughly 10% increase.

4.1. Monte Carlo Radiative Transfer

4.2. The Specific Intensity: Effects of Inhomogeneities

Figures 7–9 show the results of the MC calculations described above for the 0326c, 0528a, and 1112a snapshots, respectively. For the moment, we concentrate on the top panels of each figure, which show the specific intensity at an inclination of 79°. As one would expect on physical grounds, the specific intensity shows little azimuthal variation, so we plot the azimuthally averaged output of the MC calculation to improve the statistics. The solid black curves correspond to the results of MC calculations through the full three-dimensional grid and the red curves are the MC spectra from horizontally averaged one-dimensional grids. In each case both sets of calculations yield rather similar spectral shapes, and the spectral peaks occur at nearly identical energies. The main difference is that the three-dimensional model intensity is greater at all but the lowest energies. Although we have only shown the spectra at 79°, this result qualitatively holds at all inclinations, implying that a larger flux is being radiated in the three-dimensional calculations, in good agreement with previous results (Davis et al. 2004).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This flux enhancement is primarily the result of the density inhomogeneities. Although there are some horizontal variations in the temperature, they are generally smaller than the variations in the density. The weaker variation in temperature results from the relatively short photon diffusion timescale, which efficiently smooths out temperature fluctuations (see e.g., Turner et al. 2005). The large density inhomogeneities are the result of highly variable magnetic forces that dominate the surface layers. Densities range over factors of roughly ten above and below the mean at the horizontally averaged effective photospheres in all the simulations (Blaes et al. 2007; Hirose et al. 2009). It is possible that photon bubbles, if they exist, may make the density even more inhomogeneous (Turner et al. 2005) in the radiation pressure-dominated regime.

The density inhomogeneities enhance the total emission because free–free emission (and thermal emission processes in general) scales as the square of the density. Therefore, even though average density is identical in the two calculations, the excess emission from high density regions more than compensates for the deficit in low density regions. Of course, free–free absorption also scales as the square of the density and if it were the only opacity source, increased absorption would offset the increased emissivity. However, since electron scattering is proportional to one power of density and dominates the opacity at most frequencies, photons emitted in high density regions scatter out to lower densities before they are absorbed, increasing their likelihood of escape.

Ultimately, we would like to address whether or not these inhomogeneities are essential to obtaining the correct specific intensity. The fact that both the one-dimensional averaged domain and the three-dimensional domain yield the same spectral shape suggests that the spectra may be insensitive to the inhomogeneities. However, the flux enhancements suggest that we might be able to get an equivalent flux with a lower surface temperature, possibly lowering the spectral hardening. The problem with comparing the one-dimensional horizontally averaged and full three-dimensional MC calculations is that they do not represent the same radiative equilibrium. In a real disk the flux of radiation (and therefore the radiative equilibrium) is fixed by the overall amount and location of turbulent energy dissipation. One needs to compare the three-dimensional calculation against a one-dimensional calculation with the same radiative equilibrium. Fortunately, we have already described such calculations in Section 3.

Figure 10 shows a comparison of the three-dimensional MC spectrum (bins) from Figure 9 with an equivalent TLUSTY atmosphere spectrum (solid curve). The TLUSTY model includes the simulation derived dissipation profile and magnetic support, but differs from the one shown in Figure 4 (the dot-dashed curve) in that we have dropped the bound-free opacity and used a slightly lower effective temperature to better match the three-dimensional MC calculation. The lower effective temperature is needed because the total flux in the MC calculation is slightly lower than that found directly by the simulation. We find that the three-dimensional MC spectrum is significantly softer than the one-dimensional TLUSTY spectrum. In fact, the peak is lower by about 12%, almost completely offsetting the hardening due to the magnetic support discussed in Section 3.

Zoom In Zoom Out Reset image size

In fact, a similar cancellation also occurs in the 0326c and 0528a snapshots. Again, we repeated the one-dimensional TLUSTY calculations from Section 3, but neglecting bound-free opacities and fixing the effective temperature to match the MC results. The spectral hardening due to magnetic support was about the same level as the softening due to inhomogeneities in both cases: ∼10% for 0326c and ≲5% for 0528a. Since the magnetic fields providing the support are also the primary source of the inhomogeneities, a correlation between the effects is not surprising. However, it is remarkable that the magnitude of the effects is such that they cancel out in all three sets of calculations.

We plot a comparison of the temperature profiles from the two calculations in Figure 11. The dotted curve shows the horizontally averaged temperature from the 1112a simulation snapshot and the solid curve is the one-dimensional atmosphere model. Overall, the one-dimensional profile is a better match to the simulation than the equivalent model in Figure 3 (dot-dashed curve) since we have modified the input parameters to better match the MC calculation. Nevertheless, the temperature gradients differ near the surface and the effective photosphere (square symbol) is ∼15% higher than the simulation at the same height above the midplane. We generated similar plots for the other two snapshots and they are qualitatively consistent with Figure 11.

Zoom In Zoom Out Reset image size

This result, along with the flux enhancements discussed above, strongly suggests that the inhomogeneities allow the same flux to be radiated with average lower surface temperatures, leading to a softer spectrum. However, one must note that a number of other differences between the two calculations may also contribute to the resulting discrepancies in the temperature profiles and spectra. The treatment of radiation transport in the simulations (gray opacity, flux limited diffusion) differs from both the MC and the TLUSTY calculations. Also, non-hydrostatic motions are significant in the surface layers of the simulations while the one-dimensional TLUSTY calculation strictly enforces hydrostatic equilibrium. Distinguishing between these effects is difficult, but it is clear that the three-dimensional and dynamic nature of the simulations plays a non-negligible role in determining the disk spectrum. Ignoring these effects will lead to an underestimate of the spectral hardening.

For viewing angles inclined to the surface normal, homogeneous atmospheres will appear moderately linearly polarized (up to 12%) at photon energies for which electron scattering is nearly elastic and dominates the opacity (see, e.g., Chandrasekhar 1960). Since electron scattering opacity is typically dominant for the energy ranges of interest, the MC calculations do produce significant linear polarization at moderate to high inclinations. This can be seen in Figures 7–9, which plot the polarization viewed from an inclination of 79°.

These polarization results are most easily discussed in terms of the Stokes parameters corresponding to the total specific intensity I and the linearly polarized intensities Q and U. (Electron scattering does not impart circular polarization, so the Stokes parameter V is identically zero in our calculations.) We will also find it convenient to define the total polarization P ≡ (Q² + U²)^1/2. In this paper we always plot the normalized Stokes parameters (q ≡ Q/I, u ≡ U/I) and the degree of polarization p ≡ (q² + u²)^1/2.

We define Q such that positive and negative Q correspond (respectively) to polarization perpendicular and parallel to the surface normal. U corresponds to polarization defined relative to axes rotated by 45°. For a homogeneous domain, U would also be identically zero by symmetry. Although the simulation domains are not homogeneous, they are still highly symmetric, and U is nearly zero in all three snapshots for all viewing angles. Therefore, we only plot Q in Figures 7–9. Like the intensity, we do not observe significant azimuthal variation in the polarization, so we plot only the azimuthal average. The only exception is when the effects of Faraday rotation are included, as discussed in Section 5.3.

The direction of the polarization is generally determined by the angular dependence of the emergent radiation field. Limb darkened radiation yields polarization in the plane of the disk (positive Q), isotropic radiation is unpolarized, and limb-brightened radiation provides polarization parallel to the surface normal (negative Q). Since Thomson scattering dominated atmospheres are limb darkened, one nominally expects polarization with positive Q (Chandrasekhar 1960). However, a number of other processes are relevant: absorption at low energies, changes in photon energy due to inelastic electron scattering, Faraday rotation by the magnetic fields, and geometric effects due to the inhomogeneities.

Absorption opacity changes the polarization both by competing with scattering opacity and by altering the anisotropy of the radiation field. Depending on the vertical gradient of the thermal source function, this can increase the polarization compared to a pure scattering atmosphere (Harrington 1969; Loskutov & Sobolev 1979; Bochkarev et al. 1985), or decrease it and even turn it negative (the "Nagirner effect," c.f. Gnedin & Silant'ev 1978). For the free–free opacity assumed here, absorption dominates at lower photon energies, and the Nagirner effect is clearly evident in the bottom panel of Figure 7 at energies below 0.03 keV. It is also present in the other two snapshots, but at lower energies than those plotted in Figures 8 and 9.

5.1. Compton Scattering

The effects of inelastic (Compton) electron scattering manifest themselves as a turnover in Q at photon energies well above the spectral peak. If electron scattering is treated in the elastic (Thomson) limit the polarization would remain positive. We have deliberately neglected to plot this turnover in Figures 7–9 due to the poor statistics at these energies. But, the effect can be seen in Figure 12 where we compare a full three-dimensional model where electron scattering is treated in the Compton limit (black curve) with one in which the it is computed in the Thomson limit (red curve).

Zoom In Zoom Out Reset image size

At high energies above or near the peak, the predictions for polarization differ significantly. The result is not due to the differences in the scattering matrices, since the full Compton-scattering matrix is well approximated by the Rayleigh matrix at the energies of interest; see the Appendix for further details. Instead, it is related to the different frequency dependence of the angular distribution of the radiation in the two cases. Since the photon energies are fixed in the Thomson model, high energy photons tend to be preferentially emitted deeper in the atmosphere where temperatures are greater. Since photons are initially emitted at large scattering depths, the radiation field becomes limb darkened, the spectrum becomes a modified blackbody (Rybicki & Lightman 1979), and the polarization is well approximated by the predictions for a homogeneous, semi-infinite, scattering-dominated atmosphere (see, e.g., Chandrasekhar 1960).

However, when Compton effects are included these high energy photons emitted at large depth can now efficiently exchange energy with the electrons near the surface. Since the surface is cooler, these electrons are typically lower in energy and photons will tend (on average) to lose energy until they have a mean energy similar to the photospheric electrons. This leads to the Wien tail at high photon energies seen in the specific intensity for the Compton spectra. The photon angular distribution and the energy distribution are now coupled by the scattering history, and the polarization can therefore differ from that of an optically thick Thomson scattering atmosphere. Indeed, the polarization of 2 keV photons is slightly enhanced by Compton scattering in Figure 12.

Near the surface, photon diffusion smooths out most, but not all, temperature inhomogeneities. In a few localized regions the gas temperature significantly exceeds the effective temperature of the radiation, though this is likely artificial as simulation 0528a neglected Compton energy exchange between the radiation and the gas (Blaes et al. 2007). In these pockets, photons are upscattered to higher energies, producing the break seen at high energy in the top panel of Figure 12. Due to the small optical depth of the hot material, these photons only overcome the emission in the Wien tail at high energies. At these energies the up scattered emission is modestly limb brightened, due to the slightly increased probability of laterally propagating photons to be scattered before escape. This is completely analogous to the polarization flip produced by inverse Compton scattering in an optically thin, hot corona in plane parallel geometry (Haardt & Matt 1993; Hsu & Blaes 1998), and is a reminder that such coronal geometries, if they exist in nature, will produce high energy radiation that is polarized parallel to the projected disk rotation axis.

5.2. Inhomogeneities

The effects of inhomogeneities can be seen in the bottom panels of Figures 7–9. The black curves show the polarization of the full three-dimensional calculations (excluding the effects of Faraday rotation) while the red curves correspond to the one-dimensional, horizontally averaged domain. We find that the three-dimensional polarization curve is generally less polarized than the one-dimensional horizontally averaged curve, although the effect is rather small.

This slight reduction in polarization is in qualitative agreement with the expectations of Coleman & Shields (1990), who argued that deviations from a planar photosphere (e.g., a corrugated surface) might produce a more isotropic, or even limb brightened, radiation field. They postulated that such a photospheric geometry could, in principle, explain the apparent results that AGNs are at most weakly polarized (≲2%) and that the polarization is perpendicular to the inferred disk plane, contradicting the standard expectations of scattering dominated atmospheres. However, the effect observed in our calculations is much smaller than suggested by Coleman & Shields (1990) and may partly be the result of the increased ratio of absorption to scattering opacity at some energies.

The model of Coleman & Shields (1990) implicitly assumes that the photosphere geometry is complex, but that characteristic length scales of the horizontal variations is significantly greater than the photon mean-free path. (If the horizontal length scale is smaller or comparable with a mean-free path, their assumption of a Chandrasekhar (1960) limb-darkening profile relative to the local surface normal would be invalid.)

To compare with their model, we examined the shape of the photosphere in our three-dimensional calculations. In one-dimensional atmospheres, the photosphere corresponds to a surface of constant escape fraction, a concept which is easily generalized to three-dimensional. Therefore, we computed in each grid zone the fraction of photon packets escaping after each scattering or emission event. In Figure 13 we plot the contours of constant escape fraction as white curves on top of the density in a two-dimensional slice of the 0528a snapshot. Although deviations from a planar surface are significant, the length of horizontal variations is comparable with or smaller than a typical mean-free path to electron scattering. Therefore, the geometric effect discussed by Coleman & Shields (1990) is nearly negligible for these simulation domains.

Zoom In Zoom Out Reset image size

Nevertheless, these results do not invalidate the Coleman & Shields (1990) hypothesis for real systems. Since periodic boundary conditions are assumed, larger length scale horizontal variations are not possible in the current simulations. The surface inhomogeneities are primarily the result of Parker instability (Blaes et al. 2007), so larger wavelength variations may be possible in actual accretion flows, and this hypothesis should be reexamined when larger domains become computationally feasible.

5.3. Magnetic Fields

The bottom panels of Figures 7–9 also show the effects of magnetic fields on the polarization. Both the blue and black curves are three-dimensional calculations, but with and without (respectively) the effects of Faraday rotation included. Each time a photon packet is propagated through a Thomson optical depth τ_T, the polarization vector is rotated by an amount

$$\begin{eqnarray} &&\chi _{\rm F} = \frac{3 \lambda ^2 \tau _{\rm T}}{16\pi ^2 e} {\bf B}\, {\bm \cdot}\,{\bf \hat{k}}, \end{eqnarray} \tag{ 3 }$$

where B is the magnetic field in the current grid zone, $\bf {\hat{k}}$ is a unit vector in the direction of photon propagation, e is the electron charge, and λ is the wavelength. For the observed polarization, the rotation that occurs after the last scattering is typically most important, so τ ∼ 1. Since $\bf {\hat{k}}$ is a unit vector, the typical rotation angle only depends on λ, B, and the direction. Due to the λ² dependence in Equation (3), one is normally not concerned with Faraday rotation in X-ray sources. However, the near-equipartition strengths of magnetic fields in the simulation (∼10⁵–10⁶ G) are large enough that Faraday rotation can be significant even at these short wavelengths.

For any single photon, Faraday rotation can only change the direction of the polarization, not its magnitude. However, due to departures from uniformity in the magnetic field and differences in trajectories, different photons will experience rotation through different angles. This is particularly significant for low energies where χ_F ≫ 1 and small differences along nearby trajectories can lead to rather different polarization angles. As a result, photons are imparted with a nearly uniform distribution in polarization, yielding zero net polarization after integration. Therefore, Faraday rotation tends to completely depolarize at low energies, gives modest reductions at intermediate energies where χ_F ≳ 1, and has a negligible effect for higher energies where χ_F ≪ 1.

For photon energies with χ_F ∼ 1, there is an azimuthally dependent rotation, due to the nearly toroidal net magnetic field in the snapshot. Figure 14 shows the polarization angle

$$\begin{eqnarray} &&\psi \equiv 0.5\tan ^{-1}(U/Q) \end{eqnarray} \tag{ 4 }$$

plotted as a function of azimuthal angle ϕ for the 0528a snapshot. A ψ of zero corresponds to polarization parallel to the simulation domain surface. This plot shows the polarization angle for an inclination of 79°, averaged over photon energies from 0.3–1 keV, where χ_F ∼ 1. For ϕ ∼ 110° and 290° the escaping photons are traveling parallel and antiparallel to the net toroidal field, obtaining non-zero polarization angles. When averaged over all azimuth, ψ is consistent with zero.

Zoom In Zoom Out Reset image size

The observed spectrum from an unresolved disk is a weighted average over all azimuth. It is "weighted" because relativistic beaming can enhance emission from the approaching side and diminish emission from the receding side of a disk annulus. Therefore, cancellation will not be complete deep in the relativistic potential well of a black hole, so some net rotation may be observed. Nevertheless, we have chosen to plot the uniformly azimuthally averaged q in Figures 7–9, since the details of the relativistic beaming are complicated to model, and depend precisely upon where the simulation would be located in a real accretion flow. For higher or lower photon energies, there is very little net rotation, and U remains consistent with zero, even before azimuthally averaging the output.

In the bottom panels of Figures 7–9, the dashed curve corresponds to the simple analytic relation

$$\begin{eqnarray} &&q=q_0(\mu)[1+(\lambda /\lambda _B)^2]^{-1}, \end{eqnarray} \tag{ 5 }$$

where q₀(μ) is the non-magnetized polarization value appropriate for the observed inclination in a scattering dominated atmosphere (Chandrasekhar 1960), λ²_B ≡ 8π³e/(3B_ph), B_ph is the average magnitude of B at the photosphere of the snapshot, and μ is the cosine of the inclination angle. The quantity (λ/λ_B)² is roughly the Faraday rotation angle for unit optical depth. For energies near or below the peak, this functional form provides a rather good approximation to the three-dimensional MC results when Faraday rotation is included. It breaks down at high energies because it does not account for the Compton scattering turnover, but is a nearly perfect match to the Thomson scattering calculations.

Figures 7–9 are all plotted for the same inclination of 79°, a relatively high value at the edge of the observed distribution for black hole X-ray binaries (see, e.g., Narayan & McClintock 2005). Since polarization increases with inclination, it will generally be lower in most sources and the effects discussed here will be more difficult to measure. The inclination dependence of Q is shown in Figure 15. The black, blue, and dashed curves have the same meaning as in the bottom panels of Figures 7–9.

Zoom In Zoom Out Reset image size

These results can be compared with previous work on one-dimensional atmospheres with uniform (Silant'ev 2002; Shternin et al. 2003; Gnedin et al. 2006) and turbulent (Silant'ev 2007) magnetic fields. Our ad hoc relation (5) is qualitatively similar to the polarization dependence derived in these works. They find depolarization at photon energies with large Faraday rotation angles and make detailed predictions for the angular and energy dependence of the polarization. Their analytic formulas are most precise in the asymptotic limit of large Faraday rotation angles, but it is difficult for us to probe the asymptotic dependence due to limited photon statistics at low energies. For a uniform magnetic field, they find p ∝ δ⁻¹ for δ ≫ 1, where δ is half the Faraday rotation angle for optical depth of unity. This is the same dependence as in Equation (5) for λ ≫ λ_B. For isotropic turbulent fields, Silant'ev (2007) predicts p ∝ λ⁻⁴ for large λ, a somewhat stronger dependence than found here. Although the simulation magnetic fields do have turbulent fluctuations, there is a mean toroidal field that may explain our consistency with their uniform field relation.

The effects discussed in Section 5 indicate that the assumption of Chandrasekhar (1960) polarization for a semi-infinite atmosphere is not correct at either high or low photon energies. However, note that our simulations only represent a small patch of the disk and (at best) only approximate the emission from a single, narrow annulus in the accretion flow. In order to assess the possibilities for obtaining constraints from actual observations, we need to compute full disk models. Given the uncertainties and complexities involved, our goal is only an approximate model, which is presented to motivate future calculations and observations.

Following previous work (Hubeny et al. 2000, and references therein), we assumed a fully relativistic thin disk model (Shakura & Sunyaev 1973; Novikov & Thorne 1973) with zero torque inner boundary condition at the innermost stable circular orbit. We also use the KERRTRANS code (Agol 1997) to compute the relativistic effects on photon geodesics and parallel transport of the polarization vectors as they travel from the disk surface to the observer at infinity.

For simplicity, we will assume that the local spectra are color-corrected blackbodies with a color correction of 1.7, which is a qualitatively reasonable approximation at the accretion rates of interest (Shimura & Takahara 1995; Davis et al. 2005). Of course, the results presented in Section 3 (and previous work) indicate that this will not be quantitatively correct in detail. The main obstacle to a more precise calculation is that we do not have simulations for the innermost rings which account for most of the radiation, so we cannot robustly model the effects of dissipation and magnetic support for the whole disk. Based on the simulations presented in this paper (which, however, may not generalize to hotter annuli), the overall shift in the spectrum may be very small, and the numerous approximations (outlined below) ultimately make this a qualitative enterprise, anyway.

We additionally assume that the local spectra have an angular distribution that approximately matches the Chandrasekhar (1960) limb-darkening profile, although this only has a small effect on the output spectrum. Note that we are also assuming that the local rest-frame polarization is azimuthally symmetric and either parallel or perpendicular to the disk plane. We neglect the azimuthal dependence of polarization position angle discussed in Section 5.3, due to the much larger uncertainties inherent in this illustrative full disk model.

For the energy dependence of the polarization, we use Equation (5) as an analytic approximation for each individual annulus. This accounts for the effects of Faraday rotation, but omits two other polarization-changing effects that could be important. If there is a hot, optically thin corona above the disk surface, Compton scattering within the corona could rotate the polarization angle to be nearly perpendicular to the surface at energies above the thermal peak (Haardt & Matt 1993; Hsu & Blaes 1998). In the simulation data we study, the gas in the region above the thermalization photosphere is hotter than the effective temperature, but the inverse Compton scattering occurring there is much too weak to explain the coronal emission typically observed; consequently, our calculations likely underestimate the importance of this effect, i.e., overestimate the energy at which the polarization rotation takes place. Secondly, we do not consider the scattering of returning radiation (Agol & Krolik 2000; Schnittman & Krolik 2009), a global general relativistic effect that also rotates the polarization (and strengthens it) at energies above the maximum effective temperature found in the disk.

The last requirement is a specification for B_ph as a function of radius in the disk model. In the three shearing box simulations, B_ph correlates well with the equipartition magnetic field strength B_eq, the magnetic field strength for which magnetic pressure equals the total midplane pressure. Although there is some scatter in the correlation, we find that B_ph ≃ B_eq/40. Therefore, we estimate B_ph by first computing B_eq as a function of radius in our thin disk models and then use this relation to obtain B_ph.

This estimate is particularly uncertain, principally because the total midplane pressure depends on α in the thin disk model, so that our estimate of B_ph inherits this dependence: B_ph ∝ α^−1/2. We assume α = 0.01 in our models, a value which is ≃1/2 the typical long-term time average found in our stratified shearing box simulations, and roughly an order of magnitude smaller than both the typical numbers found in the disk body in global disk simulations and some observationally based estimates (e.g., King et al. 2007). Using a constant α for this purpose may create a further problem in that global disk simulations often show a sizable increase in the time average of this quantity in the region just outside the ISCO, where a significant part of the luminosity is created. An additional uncertainty is introduced by the fact that what matters for Faraday rotation is ${\bf B}\, {\bm \cdot}\, {\bf \hat{k}}$, and the relative magnitudes of the different components of B may change systematically with radius.

After combining these prescriptions, we compute full disk models, the results of which are shown in Figures 16 and 17. The top panels show the total specific intensity, the middle panels show the degree of polarization, and the bottom panels show the polarization angle ψ. Polarization parallel to the plane of the disk corresponds to ψ = 0°, while ψ = 90° corresponds to polarization parallel to the projected rotation axis of the disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It has been well established by previous work that the polarization is sensitive to the properties of the spacetime and observer inclination (see, e.g., Connors & Stark 1977; Stark & Connors 1977; Connors et al. 1980; Laor et al. 1990; Agol 1997; Agol & Krolik 2000; Dovčiak et al. 2004, 2008; Li et al. 2009; Schnittman & Krolik 2009). Here we briefly review the effect of spin on the polarization results and refer the reader to previous work for further discussion. The solid and dashed curves in Figure 16 are computed for a_* = 0 and a_* = 0.99, respectively. Increasing the spin parameter has the well-known effect of shifting the spectral peak to higher energies.

The thin curves show the effects of the spacetime on polarization most clearly. These models assume that the polarization at the disk surface is everywhere identical to the Chandrasekhar (1960) value. Because spacetime near the black hole is strongly curved, photons traveling at an oblique direction at infinity are in some cases emitted at fluid-frame directions rather closer to the polar axis. The parallel transport of polarization then leads to an effective dilution of the observed polarization. This effect is strongest at small radii, where the highest-energy photons are predominantly created, so the polarization diminishes somewhat toward higher energy. There is also a net rotation of the polarization angle which can be seen in the bottom panel. For observer-viewing angles as shown, which are less than 90° to the disk rotation axis (defined in a right-handed sense), the sense of the rotation of the polarization position angle is clockwise as one moves to higher photon energies (Connors et al. 1980). Note, however, that all our calculations omit another relevant relativistic effect: the scattering of returning radiation. As shown by Agol & Krolik (2000) and in greater detail by Schnittman & Krolik (2009), this effect can lead both to a rise in polarization at high energies and a nearly 90° rotation of its direction.

The thick curves in the bottom two panels of Figure 16 show the results from disk models for which Equation (5) specifies the polarization at the disk surface. The low energy polarization is significantly reduced by Faraday depolarization and is negligible at or below 0.1 keV. The degree of polarization increases with photon energy, but even near the peak, the polarization is reduced significantly from the non-magnetized Thomson scattering models.

From Figure 16 it is clear that Faraday rotation has a significant effect for the magnetic field strengths adopted in our model. However, it may be the case that our assumed B_ph ≃ B_eq/40 ≡ B₀ does not generally hold or that our estimate of B_eq from the thin disk model may be either higher or lower than in real flows. Therefore, we explore the sensitivity to this assumption in Figure 17 by plotting several spectra with identical spin, mass, accretion rate, and inclination, but with different assumptions for magnetic field strength.

At the high energy end, all the models asymptote to a similar polarization set by the Chandrasekhar value for this inclination, but with a slight reduction due to relativistic effects. The effects of varying the magnetic field are seen at energies near the spectral peak (in νI_ν) and at low energies. Quite generally, the level of polarization decreases at all energies as we strengthen the magnetic field. The wavelength dependence is also sensitive to field strength. The model with B_ph = 0 (triple-dot-dashed curve) shows a monotonic decline in polarization as photon energy increases. In the models with lower but non-zero magnetic field (B_ph = B₀/10, B₀/3), the degree of polarization generally increases as photon energy increases and reaches a maximum or begins to level off at or just below the spectral peak. For even higher fields (B_ph = B₀, 3B₀), the polarization continues to rise and only levels off out in the Wien tail.

We can compare our results with those of Gnedin et al. (2006), who also examined the wavelength dependence of polarization in accretion disks subject to Faraday rotation. They computed the polarization from accretion disks with purely vertical magnetic field using the analytic models of Silant'ev (2002), but neglecting the effects of relativistic transfer. They reported the long wavelength asymptotic dependence of the polarization for various assumptions about radial dependence of the magnetic field. We find an approximate asymptotic dependence of p ∝ λ^−s with s ∼ 0.3. Long wavelength photons are predominantly emitted at larger radii where relativistic effects are negligible and P_gas > P_rad, so B_ph ∝ B_eq ∝ r^−5/4 (Shakura & Sunyaev 1973). For this radial dependence, Gnedin et al. (2006) find s = 1/3 in reasonable agreement with our result. The agreement is not surprising since the depolarization parameter (their δ) is not very sensitive to the field geometry, and the asymptotic dependence found by Silant'ev (2002) is very similar to that of Equation (5).

7.1. Implications for Spectral Hardening and Continuum Spin Estimation

7.2. Polarization Predictions and Potential Magnetic Field Measurements

7.3. Generality of these Predictions

If our assumed scaling of B_ph with B_eq generally holds, Faraday rotation is likely to be important in most luminous accretion regimes. The innermost regions of near Eddington accretion disks are radiation dominated for ∼10 M_☉ black holes, and are even more radiation dominated for supermassive black holes. Using the radiation-dominated relations of Shakura & Sunyaev (1973), we can compute the characteristic rotation angle for optical depth unity near the spectral peak using Equation (3)

$$\begin{eqnarray} &&\chi _{\rm F,p} \simeq \frac{3 \lambda _{\rm p}^2}{16\pi ^2 e} {\bf B_{\rm ph} \cdot {\hat{k}}}, \end{eqnarray} \tag{ 6 }$$

where λ_p ≃ ch/(4k_BT_eff). The effective temperature T_eff is computed from the flux

$$\begin{eqnarray} &&T_{\rm eff} = \left(\frac{3 c^3 \dot{m} R_{\rm R}}{2R_g \kappa _{\rm es} \eta \sigma r^3}\right)^{1/4}, \end{eqnarray} \tag{ 7 }$$

where R_g = GM/c² is the gravitational radius, κ_es is the electron scattering opacity, σ is the Stefan–Boltzmann constant, r is the radius in gravitational radii, η is the radiative efficiency, $\dot{m}$ is the accretion rate normalized to the Eddington rate $\dot{M}_{\rm Edd}=4\pi R_g c/(\kappa _{\rm es} \eta)$, and R_R is a function encapsulating both relativistic corrections and the ISCO torque boundary condition (notation following Krolik 1999). Defining B_eq = (8πP_rad)^1/2 and using the results of Shakura & Sunyaev (1973), we have

$$\begin{eqnarray} &&B_{\rm ph} = \frac{B_{\rm eq}}{40} = \left(\frac{8\pi }{\alpha \kappa _{\rm es}R_g}\right)^{1/2} \frac{c}{40 r^{3/4}}\left(\frac{R_{\rm z}R_{\rm T}}{R_R}\right)^{1/2}, \end{eqnarray} \tag{ 8 }$$

where R_T, like R_R, summarizes both relativistic and boundary condition corrections to the stress profile, and R_z is the relativistic correction to the vertical gravity. Combining Equations (6)–(8) with the definition of λ_p, we arrive at

$$\begin{eqnarray} &&\chi _{\rm F,p} \sim 1 \left(\frac{\eta }{\dot{m}}\right)^{1/2} \frac{r_{10}^{3/4}}{\alpha _{0.01}^{1/2}q_4^2 Q_{40}} \frac{(R_{\rm z} R_{\rm T})^{1/2}}{R_{\rm R}} {\bf \hat{b}}\, {\bm \cdot}\, {\bf \hat{k}}, \end{eqnarray} \tag{ 9 }$$

where ${\bf \hat{b}}$ is a unit vector in the direction of the mean field, q = hc/(k_BT_effλ_p), and Q = B_eq/B_ph.

Remarkably, the characteristic Faraday rotation angle is independent of black hole mass wherever the disk is dominated by radiation pressure. Moreover, at least to the accuracy of this estimate, it is always ∼O(1) for wavelengths near the thermal peak for the radii where most of the light is produced when $\dot{m} \sim \eta$. Although the scaling we have derived should be robust, the actual magnitude of the effect is subject to significant uncertainty. As we have already discussed, the magnitude of α could likely be somewhat larger than the fiducial value we have chosen for it. Color corrections could alter q. Although simulations at a variety of radii have all pointed toward Q ≃ 40, we cannot say for certain whether this ratio applies in the innermost disk. The particular radius dominating the light output depends on black hole spin and the ISCO boundary condition. Because r ≃ 10 is appropriate for non-rotating black holes with little torque at the ISCO, spin and extra stress in the inner disk would diminish the characteristic rotation angle, both by decreasing the characteristic radius of emission and by increasing R_R at small radii (typically R_T ≃ R_R, while R_z ≳ 1).

We therefore predict that Faraday rotation could substantially affect the emergent polarization both for high-accretion rate Galactic black holes (in the thermal or steep power-law states) and for many AGNs (as has been previously pointed out by, e.g., Blandford 1990). Prospects for using this effect to constrain disk field properties are better in the stellar-mass case than in the supermassive case, however, because in AGNs it appears that the observed polarization does not arise from the accretion disk itself (e.g., Antonucci 1988), but from scattering by material further out. Accretion disks in cataclysmic variables also emit the bulk of their light in the optical and ultraviolet, but these disks are not scattering dominated. The theoretically predicted continuum polarization in the outburst phase is very small (Cheng et al. 1988), and has in fact defied observational attempts to separate it out from interstellar polarization (Naylor et al. 1996; Moffat et al. 2001). X-ray polarimetry of X-ray binaries might therefore be an optimal way of constraining accretion disk magnetic fields.

Using snapshots of local shearing box simulations of accretion disks with a broad range of radiation to gas pressure ratios, we have calculated how a realistic vertical structure established by MRI turbulence affects the emergent photon spectrum and polarization. Most of the dissipation within the turbulence occurs at high effective optical depth in all the physical regimes that have been simulated thus far. As a result, the spectra are very insensitive to the details of that dissipation profile and are extremely close to standard model spectra that assume that the vertical profile of dissipation per unit mass is constant. On the other hand, standard models usually neglect the fact that the photospheric regions are supported against gravity by magnetic forces. These forces reduce the horizontally averaged densities at the thermalization photosphere, resulting in harder one-dimensional model spectra compared with standard models.

In addition to lifting the atmosphere, those same magnetic forces also produce complex three-dimensional density inhomogeneities in the photospheric regions. Through a comparison between three-dimensional Monte Carlo calculations and one-dimensional calculations that were designed to have the same emergent flux, we demonstrated that these inhomogeneities act to soften the spectrum. Somewhat surprisingly, this softening largely cancels the hardening due to vertical magnetic support in all three snapshots, even though the magnitude of each effect differs between the three simulations. Therefore, it may be that the color corrections derived from, e.g., BHSPEC models (Davis & Hubeny 2006) will be approximately correct. However, we cannot rule out the possibility that this cancellation was fortuitous in the three snapshots we considered here. Given the magnitude of the effects we have found, we caution observers that color correction factors derived from atmosphere models that neglect magnetic support and three-dimensional inhomogeneities may be incorrect by approximately 10% in either direction.

Because the density inhomogeneities produce a complex, structured photosphere, the degree of polarization of the emerging radiation field is reduced compared with a plane-parallel atmosphere, but only slightly. Faraday rotation by the strong magnetic fields in the atmosphere produces a much stronger effect, producing a significant reduction in polarization for photon energies near or below the peak in the local spectrum. This fortunate coincidence that Faraday rotation is strong, but not too strong, means that future X-ray polarimetry measurements of the polarization of the thermal component of black hole X-ray binaries could be used as a diagnostic of magnetic fields in disks. More extensive theoretical calculations will be necessary, however, before this can be done quantitatively.

We thank J. Schnittman for useful discussions. This work was supported in part by NSF grant AST-0707624. S.D. is supported by NASA grant PF6-70045, awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060.

Berestetskii et al. (1982) provide an incisive summary and first-principles derivation of the cross sections for both polarized and unpolarized Compton scattering. Polarized radiation can be completely described by a Stokes vector I ≡ (I, Q, U, V)^T, where I is the total intensity, Q and U are intensities of linearly polarized radiation, and V is the intensity of circular polarization (Chandrasekhar 1960). We only consider polarization arising from electron scattering, so V = 0 and we will drop it in what follows. Since we only consider scattering of individual photon packets, it will be convenient to work in terms of normalized Stokes parameters q ≡ Q/I and u ≡ U/I. Note that these quantities correspond to ξ₃ and ξ₁, respectively, in the notation of Berestetskii et al. (1982).

We consider the general problem of modeling electron scattering in the comoving frame of the fluid (as opposed to the electron rest frame). We define photon and electron momentum-four vectors

$$\begin{eqnarray} K_\mu &=& \omega m_{\rm e} c (1,{\bf k}) \nonumber \\ P_\mu &=& \gamma m_{\rm e} c (1,{\bbeta }), \end{eqnarray} \tag{ A1 }$$

where ω ≡ hν/(m_ec²) and all other variables have their usual meanings. It is convenient to work with relativistic invariants x ≡ 2P^μK_μ/(m_ec)² and x' ≡ 2P^μK'_μ/(m_ec)², where primes denote the scattered quantities. Conservation of momentum then requires

$$\begin{eqnarray} &&x^{\prime }=x - 2 \omega \omega ^{\prime } (1- \cos \theta), \end{eqnarray} \tag{ A2 }$$

where cos θ ≡ k · k' is the fluid frame scattering angle.

The differential scattering cross section can be summed over final photon polarizations to yield

$$\begin{eqnarray} &&\frac{d \sigma }{d \Omega ^{\prime }}=\frac{3 \sigma _{\rm T}}{16 \pi }\left(\frac{x^{\prime }}{x}\right)^2 \left[\left(\frac{x^{\prime }}{x}+\frac{x}{x^{\prime }}\right)-4 \delta (\delta +1)(1-q_{\rm in}) \right],\qquad \end{eqnarray} \tag{ A3 }$$

where δ ≡ 1/x − 1/x' and q_in is the Stokes parameter evaluated in the "internal" basis defined by the scattering plane. In general, this differs from the "external" basis defined relative to the fluid frame. They are related to each other through rotation matrices L(χ), which transform I to the internal basis before scattering, and L(−χ'), which transforms I' back to the external basis after scattering. The matrices are given by Chandrasekhar (1960):

$$\begin{eqnarray} L(\chi)=\left(\begin{array}{@{}ccc@{}} 1 \;& 0 \;& 0 \\ 0 \;& \;\;\cos {2\chi } \;& \sin {2\chi } \\ 0 \;& -\!\sin {2\chi } \; & \cos {2\chi }\end{array}\right) \end{eqnarray} \tag{ A4 }$$

Using these rotation matrices, one finds q_in = qcos 2χ − usin 2χ

Nagirner & Poutanen (1993) provide a thorough discussion of polarization bases and derive expressions for the angles χ and χ' in their appendices. Here, we just quote their results

$${\fontsize{9.0}{11.0}\selectfont{ \begin{eqnarray} \cos \chi & = & \frac{\gamma (k_z^{\prime }-\cos \theta k_z)-\beta \gamma [(1-\cos \theta) \hat{\beta }_z+ {\bf k}\, {\bm \cdot }\, \hat{\bbeta } (k_z^{\prime }-k_z)]}{\sqrt{1-k_z^2} \sqrt{(\delta -1)/\delta } (1-\cos \theta)}, \nonumber \\ \sin \chi & = & \frac{\gamma (1-\bbeta\, {\bm \cdot}\,{\bf k})(k_y k_x^{\prime }-k_x k_y^{\prime }) -\bbeta \gamma (1-\cos \theta)(\hat{\beta }_x k_y-\hat{\beta }_y k_x)}{\sqrt{1-k_z^2} \sqrt{(\delta -1)/\delta } (1-\cos \theta)},\nonumber\\ \cos \chi ^{\prime } & = & \frac{\gamma (k_z-\cos \theta k_z^{\prime })-\beta \gamma [(1-\cos \theta) \hat{\beta }_z+ {\bf k^{\prime }}\,{\bm \cdot }\, \hat{\bbeta } (k_z-k_z^{\prime })]}{\sqrt{1-k_z^{\prime 2}} \sqrt{(\delta -1)/\delta } (1-\cos \theta)}, \nonumber\\ \end{eqnarray}}} \tag{ A5 }$$

$${\fontsize{9.0}{11.0}\selectfont{ \begin{eqnarray} \sin \chi ^{\prime } & = & \frac{\gamma (1-\bbeta\, {\bm \cdot}\, {\bf k^{\prime }})(k_y^{\prime } k_x-k_x^{\prime } k_y) -\beta \gamma (1-\cos \theta)(\hat{\beta }_x k_y^{\prime }-\hat{\beta }_y k_x^{\prime })}{\sqrt{1-k_z^{\prime 2}} \sqrt{(\delta -1)/\delta } (1-\cos \theta)}.\nonumber\\ \end{eqnarray}}} \tag{ A6 }$$

Here, $\hat{\beta } \equiv \bbeta /\beta$ is a unit vector in the electron momentum direction.

The effects of scattering on the Stokes vector are accounted for by the matrix R = L(χ')SL(−χ) where F is the transformation induced by Compton scattering in the internal basis (see, e.g., Nagirner & Poutanen 1993; Berestetskii et al. 1982):

$$\begin{eqnarray} \bf {S}=\left(\begin{array}{@{}ccc@{}} S_a+S_b & S_c & 0 \\ S_c & S_a & 0 \\ 0 & 0 & S_d\end{array}\right), \end{eqnarray} \tag{ A7 }$$

where S_a = 4δ(δ + 1) + 2, S_b = (x'/x + x/x') − 2, S_c = S_a − 2, and S_d = 4δ + 2. In the β, β' → 0 limit, x' → 2ω', x → 2ω, and δ → (cos θ − 1)/2 so that S_a → cos²θ − 1, S_b → 0, S_c → cos²θ + 1, S_d → 2cos θ. Therefore, S corresponds to the Rayleigh matrix for Thomson scattering in the appropriate limit. Carrying out the matrix multiplication, we find

$$\begin{eqnarray} \bf {R}&=&\left(\begin{array}{@{}cc@{}} S_a+S_b & S_c \cos 2\chi\\ S_c \cos 2\chi ^{\prime } & S_a \cos 2\chi \cos 2\chi ^{\prime }+ S_d \sin 2\chi \sin 2 \chi ^{\prime }\\ -S_c \sin 2\chi ^{\prime } & -S_a \cos 2\chi \sin 2\chi ^{\prime }+ S_d \sin 2\chi \cos 2 \chi ^{\prime } \end{array}\right.\nonumber\\ &&\left. \begin{array}{@{}c@{}} -S_c \sin \chi \\ -S_a \sin 2\chi \cos 2\chi ^{\prime }+ S_d \cos 2\chi \sin 2 \chi ^{\prime }\\ S_a \sin 2\chi \sin 2\chi ^{\prime }+ S_d \cos 2\chi \cos 2 \chi ^{\prime } \end{array}\right) \end{eqnarray} \tag{ A8 }$$

A.1. Monte Carlo Implementation

Pozdniakov et al. (1983; hereafter PSS83) describe and evaluate Monte Carlo methods for Compton scattering of unpolarized radiation. Using the above relations, it is straightforward to generalize the methods presented in PSS83 to include polarization. For the sake of brevity, we focus only on the additional complexity offered by the polarization physics and refer the reader to PSS83 for further details.

The angle-integrated cross section σ(x) is independent of polarization, and equivalent to Equation (2.10) of PSS83. Therefore, the mean-free path and selection of the electron momentum can be evaluated using the same method as presented in Sections 9.4 and 9.5 of PSS83. If the photon momentum is accepted, the scattered photon direction k' is drawn randomly, and cos θ and x' are evaluated as before. However, the scattering probability (Equation (9.8) of PSS83) now depends on the Stokes parameters through Equation (A3). The quantity X (defined in Section 2.1.2 of PSS83) is replaced by the quantity in square brackets in Equation (A3). Therefore, we must evaluate q_in, which requires calculation of χ via Equations (A5).

If accepted, the photon direction and energy are updated as in PSS83. Additionally, the normalized Stokes parameters after scattering must be calculated. Using I' = R · I and Equation (A8) we find

$$\begin{equation} \begin{array}{c} \dsty q^{\prime } = \frac{R_{21} + R_{22} q + R_{23} u}{R_{11} + R_{12} q +R_{13} u}, \\ \ms \dsty u^{\prime } = \frac{R_{31} + R_{32} q + R_{33} u}{R_{11} + R_{12} q +R_{13} u}.\end{array} \end{equation} \tag{ A9 }$$