text:
# DYNAMICS OF FE-NI BUBBLES IN YOUNG SUPERNOVA REMNANTS

## 1 INTRODUCTION

The presence of large-scale mixing of ejecta in core-collapse supernovae (SNe) has been well established on both observational and theoretical grounds. The early, unexpected emergence of X-rays and $`\gamma `$-rays shortly after the explosion of SN 1987A provided dramatic evidence for mixing of radioactive Ni throughout the He core and H-rich envelope. Extensive mixing of Ni is also required to explain the light curve and spectral observations of SN 1987A (see McCray 1993 for a review of SN 1987A). From studies of light curves and spectra of numerous core-collapse SNe, it is now known that such large-scale mixing is nearly always present. Another line of evidence is provided by studies of meteoritic graphite grains which condensed in SN ejecta (Travaglio et al. 1999 and references therein). Such grains apparently contained radioactive, short-lived isotopes such as 44Ti at the time of their formation, which implies extensive mixing of freshly-synthesized material deep in the SN ejecta with the C-rich layer at the bottom of the He core. This overwhelming observational evidence for mixing is consistent with our present understanding of the collapse and explosion of massive stars. The neutrino-driven Rayleigh-Taylor instability just outside the proto-neutron star leads to convective motions, strengthening the post-bounce shock and maybe even making the explosion possible in many cases (e. g., Mezzacappa et al. 1998; Kifonidis et al. 2000). Further mixing occurs later during the SN explosion because of Rayleigh-Taylor instabilities generated at interfaces between stellar layers with different chemical composition (e.g., Fryxell, Arnett, & Mueller 1991; Hachisu et al. 1992; Kifonidis et al. 2000).

Large-scale turbulence generated during the explosion appears not to be sufficient to mix SN ejecta completely, but only macroscopically. In Cas A, optical observations revealed the presence of ejecta with very different chemical abundances, such as O-, S-, Ar-, and Ca-rich ejecta knots, which were macroscopically mixed during the SN explosion as evidenced by the lack of spatial stratification expected in the absence of mixing (e. g., Fesen & Gunderson 1996). The most recent IR and X-ray observations of Cas A provide further evidence for this macroscopic but not microscopic mixing (Arendt et al. 1999; Douvion et al. 1999; Hughes et al. 2000; Hwang et al. 2000). Spectroscopic observations of numerous SNe also confirm the presence of clumpy ejecta (Spyromilio 1994). Detailed studies of SN 1987A lead to the same conclusion.

In SN 1987A, the inhomogeneities take a particularly interesting form: Fe-Ni bubbles (Li et al. 1993), inferred from observations of Fe, Co, and Ni lines. Ni-rich ejecta clumps, transported outwards by turbulent motions, are heated by their own radioactive energy input, and expand in the ambient substrate of other heavy elements, forming low-density Fe bubbles. The SN structure is then expected to be like a Swiss cheese, with Fe bubbles occupying a substantial ($`0.5`$) fraction of the ejecta volume (Li et al. 1993; Basko 1994). These bubbles should lead to vigorous turbulence and mixing in young supernova remnants (SNRs) in general.

Most of ejecta in young SNRs can be seen in X-rays, where modern X-ray observatories such as Chandra and XMM-Newton have begun to provide high-quality information. Interpretation of these observations must be done in the framework of appropriate hydrodynamical models. But realistic simulations for young SNRs, which include the presence of Fe bubbles in SN ejecta, are lacking. We report here on such 3-dimensional hydrodynamical simulations (preliminary results based on 2-D simulations were reported by Borkowski et al. 2000). Our goal is to study the basic hydrodynamics of the interaction by investigating a single Fe bubble, and to explore how the presence of multiple Fe bubbles changes the global dynamics of young SNRs.

## 2 HYDRODYNAMIC SIMULATIONS

We use a parallel version (implemented with Message Passing Interface) of the Virginia Hydrodynamics (VH-1) numerical code to study the complex dynamics of SNRs with Fe bubbles. The simulations were computed on a spherical grid of $`400^3`$ zones covering an angular span of 1.6 steradians. This resolution is comparable to the simulations by Chevalier at al. (1992) that used a grid of $`256^2`$. A higher resolution simulation would show more small scale mixing, but the large-scale dynamics of the problem would not change Chevalier et al. (1992). Periodic boundary conditions were applied in both the $`\theta `$ and $`\varphi `$ directions. This choice for $`\theta `$ is somewhat unorthodox, but we wished to avoid the numerical artifacts associated with the coordinate singularity at $`\theta =0`$, so we centered our grid about the equator ($`0.3\pi <\theta <0.7\pi `$). We experimented with other boundary conditions (e.g., reflecting or zero gradients), but they proved unsatisfactory for various reasons. The numerical grid was expanded to follow the forward blastwave so the evolution could be tracked for many expansion times.

The simulations were initialized with a self-similar driven wave (SSDW) solution Chevalier (1982) with an ejecta density power law of $`n=9`$ expanding into either a uniform ambient density (an ambient density power law of $`s=0`$) or a relic stellar wind (an ambient density power law of $`s=2`$). The radial boundary conditions were set to match these power laws. This self-similar structure has the following components: an outer blast wave, a wavy (dynamically unstable) contact discontinuity between the shocked ambient medium and the shocked SN ejecta, and an inner reverse shock. For these parameters, the forward shock decelerates according to $`R_st^{6/7}`$ ($`t^{6/9}`$ for $`s=0`$), and the reverse shock is located at a radius of $`0.785R_s`$ ($`0.842R_s`$ for $`s=0`$). The initial conditions are spherically symmetric, so the narrow shell of shocked ejecta starts off smooth but eventually shows signs of convective instability Chevalier et al. (1992).

## 3 DYNAMICS WITH A SINGLE BUBBLE

Hydrodynamical simulations are shown in Figure 1 of a SSDW ($`n=9`$, $`s=2`$) with a single Fe bubble in the supersonic ejecta. The bubble, with a fractional radius of 0.4 (the radius of the bubble is 0.4 times the distance of the center of the bubble from the center of the SN explosion), is 100 times less dense than the ambient SN ejecta at a radius corresponding to the center of the bubble. The pressure inside the bubble is equal to the pressure of the surrounding ejecta in order to minimize motions, but the Mach number of expansion is still relatively high within the bubble.

The initial contact between the Fe bubble and the reverse shock is illustrated in Figure 1 (left). When the outer edge of the bubble reaches the reverse shock, the reverse shock accelerates into the bubble because of its low density, while a rarefaction wave propagates back into the shocked ejecta. The pressure behind the reverse shock propagating through the bubble drops to more than an order of magnitude smaller than the nominal pressure between the forward and reverse shocks in the spherical SSDW. There is only a modest drop in pressure (compared to the spherical SSWD) in the shocked ambient medium ahead of the bubble because the SNR expansion timescale is comparable to the sound crossing time across the bubble. The high pressure of the relatively tenuous ambient gas then drives dense shocked ejecta deep into the bubble’s interior. This situation is Rayleigh-Taylor unstable, which leads to vigorous mixing of dense shocked ejecta with the low-density shocked Fe from the bubble and the shocked ambient gas. This enhanced mixing is evident in the middle-left frame of Figure 1. Note that although the reverse shock is moving inwards in this expanding frame, at no time in the simulation is gas ever moving radially inwards in the fixed frame of the explosion.

Eventually the reverse shock reflects off the bottom of the bubble (middle-right panel of Figure 1) and propagates back through the layer of mixed gas with a shock Mach number of $`1.5`$. The reverse shock is also transmitted into the dense ejecta, reforming the thin layer of shocked ejecta present in the initial SSDW. The bottom edge of the bubble is driven back out to the original SSDW by the high density of the ejecta, and in the process the extended layer of mixed gas inside the bubble is swept up and compressed into a thin layer.

The shock reflected from the bubble’s bottom eventually overtakes the blast wave and the layer of shocked ejecta is pushed back to the original radius of the reverse shock. This marks the transition back to the original self-similar stage, and the end of the transient stage associated with the bubble’s presence.

What is less apparent, but perhaps most important, in this simulation is the role of the oblique shocks along the bubble walls in generating vorticity. As the reverse shock traverses the low-density bubble, it drives an oblique shock into the wall of dense ejecta surrounding the bubble. This oblique shock generates vorticity at the bubble’s boundary, and the thin layer of shocked ejecta along the bubble’s wall quickly becomes unstable. This effect is most pronounced along the sides of the bubble where the transmitted shock is very oblique. By the time the bubble has been completely shocked, this excess vorticity has piled up in a thick ring with a radius slightly smaller than the original Fe bubble, as shown in Figure 2. It is this torus of high vorticity where the most vigorous mixing occurs. This mixing can lead to relatively large velocities, with some shocked ejecta traveling outwards 50% faster than the gas immediately behind the forward shock.

A hydrodynamical problem of a plane shock interacting with a low density spherical (or cylindrical) bubble has been of interest to the fluid dynamics community. Experiments (e. g., Haas & Sturtevant 1987) and hydrodynamical simulations (e. g., Quirk & Karni 1996) revealed how vortex lines (for cylindrical bubbles) and vortex rings (for spherical bubbles) form in this interaction process, although a complete understanding of this complex interaction (including details of mixing of bubble’s material with the ambient medium) is still lacking. Although geometry and dynamics in our simulations are more complex than in these idealized experiments and simulations, the emergence of vortex rings in our simulations is in qualitative agreement with results reported in the literature. We refer the reader to a recent review by Zabusky (1999) for further details about generation of vorticity in shock-accelerated inhomogeneous flows.

This hydrodynamical problem of a shock interacting with a low density bubble can be compared with the more often studied problem of a shock interacting with an overdense cloud (e. g., Klein et al. 1994; Xu & Stone 1995; see also Zabusky 1999 and references therein). In both cases the important long-term result of the shock interaction is the generation of vorticity and a rapid mixing of the gas inside and outside the bubble/cloud. This is demonstrated in a dramatic fashion below when we consider SN ejecta with multiple bubbles instead of a single bubble.

## 4 DYNAMICS WITH MULTIPLE BUBBLES

A more realistic model of Fe bubbles in SN ejecta should include multiple bubbles occupying a substantial fraction of the ejecta volume. We have evolved such models for both $`s=0`$ (Figure 3) and $`s=2`$ (Figure 4). In these examples, bubbles with a fractional radius of 0.2 and with a uniform density 100 times smaller than the ejecta density are distributed randomly in the ejecta, with an average volume filling fraction of $`0.5`$. We tracked the filling fraction by computing the fractional area occupied by bubbles at a radius of $`0.5R_s`$. Due to the clustering of bubbles created by the random number generator, the filling fraction went as high as 0.7 and as low as 0.3, but averaged reasonably close to 0.5.

The outer shock wave is relatively unaffected by the presence of the Fe bubbles in these simulations; the shape of the outer shock remains relatively spherical and the deceleration parameter ($`Vt/R_s`$) is within a few percent of the analytic solution of 6/7 for a spherical SSDW with $`n=9`$ and $`s=2`$ (6/9 for $`s=0`$). The time evolution of the deceleration parameter from both runs is shown in Figure 5. The SSDW has been evolved more than six orders of magnitude in radius, requiring 30,000 timesteps, in order to allow the blastwave to reach a quasi self-similar form (see the time-dependence of turbulent energy density in Figure 9).

In contrast to the forward shock, the presence of multiple bubbles dramatically alters the reverse shock and makes the interaction region very turbulent and inhomogeneous. From the results of the single bubble simulation, we expect the position of the reverse shock to vary randomly on a scale of order the radius of the bubbles. This variation is illustrated in the 2D slices of the multiple bubble simulation shown in Figures 3 and 4. The deviations can in fact be much larger than a bubble diameter because of the overlap of bubbles when the filling factor is large. In the $`s=2`$ model the reverse shock came close to the inner radial boundary at $`0.5R_s`$ at several times during the simulation. In the $`s=0`$ model we were forced to extend the radial grid beyond $`0.5R_s`$ in order to keep the reverse shock contained within the simulation domain. In addition to a strong variation in the radius of the reverse shock, the average radius of the reverse shock is significantly smaller than in the SSDW solution. The average width of the interaction region grows by $`50\%`$ during the first half of the simulations, from the analytic value of $`0.21`$ to an average value of $`0.31`$ for $`s=2`$ and from $`0.16`$ to $`0.26`$ for $`s=0`$.

Angle-averaged radial profiles of density, turbulent energy, velocity, and pressure for $`s=0`$ (Fig. 6) and $`s=2`$ (Fig. 7) quantitatively demonstrate how the presence of Fe bubbles affects the SNR structure. Density and velocity profiles are remarkably smooth when compared with analytic and numerical solutions for homogeneous ejecta. A spatially-distinct shell of shocked ejecta visible in these solutions is no longer discernible in simulations with bubbles. This nearly complete obliteration of radial structure is not surprising in view of the irregular shape of the reverse shock just discussed. The increase in the width of the interaction region mentioned above is most clearly demonstrated by the broad pressure profiles. Angle-averaged pressures are generally lower than for homogeneous ejecta because of the presence of unshocked ejecta with negligible pressure throughout the interaction region. In the vicinity of the blast wave angle-averaged pressure is also lower because of the presence of the unshocked ambient gas located between the blast wave protrusions.

The irregular shape of the reverse shock has an important effect on the dynamics of the SNR. If a region of interbubble ejecta encounters a relatively oblique reverse shock, it will be decelerated less than it would in a spherically symmetric SSDW. This dense region would then propagate through the intershock region to produce a protrusion of the forward shock. In effect, the geometry of the interbubble ejecta produced small overdense regions that behaved more like clumps within low-density ejecta than bubbles within high-density ejecta. The regions of interbubble ejecta that appeared to have the most effect on the intershock region were comparable in size to the bubbles themselves. Larger regions of ejecta without bubbles resulted in a relatively planar reverse shock, while smaller regions of ejecta between bubbles did not significantly perturb the intershock region.

However, even in the most favorable scenario, regions of dense shocked ejecta traveling through the intershock region were quickly disrupted by fluid instabilities and ultimately had a minimal impact on the forward shock. Dense clumps being slowed by the more tenuous shocked ambient gas are Rayleigh-Taylor unstable, while shear flow around the clumps creates Kelvin-Helmholtz instability. Both processes act to spread out the clump laterally and increase the drag. The result is that any dense clump creating a protrusion of the forward shock is almost immediately sheared apart, with the remaining pieces quickly advected back into the intershock region. Note that these same processes are responsible for shaping the shell of shocked ejecta in the case of spherically-symmetric ejecta Chevalier et al. (1992). In the absence of Fe bubbles, however, the drag within the intershock region prevented the shocked ejecta from getting anywhere close to the shock front.

The impact of these ejecta clumps on the shape of the forward shock is shown in Figure 8. The most pronounced deformations of the forward shock were relatively small wavelength (less than a bubble radius) protrusions driven by clumps of shocked ejecta. Some of these protrusions pushed the forward shock out an extra 15%, while longer wavelength variations in the radius of the forward shock were limited to only a few percent. Angle-averaged profiles (Figs. 6 and 7) show that $`5`$% excursions of the blast wave are common.

Wang & Chevalier (2000) arrived at similar conclusions regarding clumpy ejecta, namely that even extremely overdense clouds have relatively little effect on the forward shock. They used 2D simulations to study the evolution of a small, dense cloud placed in the ejecta of a Type Ia SN model in which the nominal ejecta are described by an exponential density profile. As long as sufficient spatial resolution was provided, their overdense clumps were flattened out and broken up by fluid instabilities before they reached the shock front. In 3D and with higher numerical resolution, this process should happen even faster.

The dynamics in these multiple-bubble simulations appear much richer than isolated bubbles or clouds alone. In the multiple bubble simulations the propagation of the reverse shock through bubbles and around interbubble regions leads to a substantial amount of turbulence in the intershock region. We tracked the turbulent energy density in the intershock region by summing up the transverse (non-radial) kinetic energy in each zone between the reverse and forward shocks. We adjusted this value by $`3/2`$ to account for the radial component of the turbulent flow (i.e., we assume all three components of the turbulent flow have comparable magnitude: $`v_r^2v_\theta ^2v_\varphi ^2`$). This average was then normalized by the kinetic energy density flowing through the shock front, $`\rho _oV_s^2`$. The results are plotted as a function of radius in Figures 6 and 7, and as a function of time in Figure 9. In contrast to the case of homogeneous ejecta, the turbulent region is no longer restricted to the shocked ejecta and the adjacent shocked ambient gas. The turbulence extends across the whole interaction region, as implied by broad, slowly varying turbulent energy profiles. The average turbulent energy density grows relatively slowly, but eventually reaches a significant fraction ($`>10\%`$) of the energy density associated with the forward shock.

## 5 OBSERVATIONAL CONSEQUENCES

### 5.1 Gas Temperature and Velocity

The irregular shape of the reverse shock has important implications regarding the X-ray spectrum produced by these SNR models. Any deviation from a spherical shock will decrease the shock velocity and hence the postshock temperature. From Figures 3 and 4 we see that this will be mildly important for the forward shock, but can be expected to have dramatic effects for the reverse shock. Furthermore, as the reverse shock moves backward through a bubble (in the expanding frame), the shock velocity can be higher than in the spherical case, resulting in higher postshock temperatures (but with small emission measure because of the low density).

To estimate these effects without undertaking the complicated effort of calculating X-ray spectra, we have summed up the emission measure (EM) for each zone in the simulation, and plotted this as a function of temperature and radial velocity (scaled to the postshock temperature and shock velocity). (This is only a very approximate procedure as we neglected variations in the mean molecular weight between Fe bubbles and the ambient ejecta. For detailed comparisons with X-ray observations one also needs to consider electron temperature which is generally lower than the mean gas temperature discussed here.)

The results for both the $`s=2`$ and $`s=0`$ multiple bubble simulations are shown in Figures 10 and 11. Superimposed on these plots are the curves corresponding to emission from spherically symmetric SSDWs. In both cases these curves are composed of two pieces, one for the shocked ejecta and one for the shocked ambient gas. Shocked gas immediately behind the forward shock is located at $`\mathrm{log}(T)=0`$ and $`V=0.75`$. For the $`s=2`$ case the gas temperature in the SSDW decreases with increasing distance from the forward shock, while the opposite holds for $`s=0`$. Shocked ejecta immediately behind the reverse shock are located at $`\mathrm{log}(T)=1.76`$ and $`V=0.82`$ and decrease in temperature away from the shock in the case of $`s=2`$. For $`s=0`$, the shocked ejecta start off at $`\mathrm{log}(T)=0.75`$ and $`V=0.95`$ and increase in temperature away from the shock. Despite these differences in the spherically symmetric models, the emission maps from the 3D simulations with multiple bubbles appear relatively similar. In both cases X-ray emitting gas is significantly spread out in temperature. This is most dramatic in the $`s=0`$ case, for which there is no cool gas in the spherical solution.

There is also a large spread in the radial velocity of the X-ray emitting gas, particularly for the case of $`s=0`$. In both simulations the radial velocity of the shocked ambient gas is spread to lower velocities as a result of protrusions in the forward shock, while the radial velocity of the shocked ejecta is spread to higher velocities. In the case of $`s=0`$, the bulk of the X-ray emitting ejecta is traveling faster than the forward shock.

### 5.2 SNR Morphology

To provide a qualitative estimate of how Fe bubbles might affect the observed morphology of SNRs, we created a volume rendering of the emission measure of shocked gas computed from the last frame in our multiple-bubble simulations. While the actual emission observed from a SNR will also depend on the local temperature, abundance, and ionization, the emission measure provides a simple, convenient means for exploring the X-ray morphology implied by a given hydrodynamical model. As such, these images are meant only to show an overall qualitative agreement with the observed X-ray morphology of SNRs like Cas A.

The rendered images shown in Figure 12 are dominated by narrow filaments of shocked ejecta. These filaments are barely resolved in our numerical simulations, and they occupy a very small fraction of the volume of the SNR. In the $`s=2`$ model, 90% of the EM comes from less than 2% of the volume of shocked gas. These filaments often show up as rings or partial rings representing the circumference of Fe bubbles as they pass through the reverse shock. These rings show up more prominently in the $`s=0`$ simulation.

## 6 SUMMARY

Our primary conclusion is that large-scale inhomogeneities associated with the presence of Fe bubbles in heavy-element ejecta lead to vigorous turbulence and mixing in young SNRs. The turbulent energy density may be increased by an order of magnitude when compared with homogeneous ejecta. This turbulence radically changes the spatial structure of the interaction region between the blast wave and the reverse shock. While the time-averaged dynamics and the propagation of the forward shock are well described by a self-similar solution Chevalier (1982); Chevalier et al. (1992), the reverse shock geometry and the structure of the shocked region are very different in the presence of the Fe bubbles. The reverse shock is no longer approximately spherically symmetric, and jet-like fingers of ejecta can even affect the location of the blast wave. Angle-averaged density profiles rise smoothly toward the remnant’s interior, with no signs for the presence of a spatially distinct shell of ejecta. Pressure and temperature variations of orders of magnitude are also present, unlike the case of homogeneous ejecta which can be described by well-defined density, pressure, and temperature profiles. Such large variations should dramatically influence X-ray emission, which depends strongly on temperature and density. In particular, low temperature, dense ejecta are particularly efficient in producing X-ray emission, and may dominate X-ray spectra. But this emission is generally produced by ejecta which suffered the least amount of deceleration, so that their radial velocities and proper motions may be even higher than that of the blast wave. This example demonstrates that standard SSDW solutions should be used with an extreme caution for young SNRs with inhomogeneous ejecta. The amount of turbulence generated by the presence of Fe bubbles in SNRs most likely varies from remnant to remnant because of large ($`100`$) variations in the radioactive Ni yields in core-collapse SNe. In addition, the amount of turbulence might be lower for those SNRs which did not have enough time to enter the self-similar turbulent regime, because of long timescales necessary to achieve it. But unless mixing of heavy elements is somehow inhibited during a SN explosion, it should be clear that one-dimensional models are clearly not acceptable for remnants of core-collapse SNe, and that multidimensional hydrodynamical modeling is essential for understanding dynamics of heavy-element ejecta in such remnants. It is also possible that Ni bubbles can be formed in Type Ia supernova explosions (Wang & Chevalier 2000), so that our results may be applicable there as well.

We expect that small-scale turbulence will have a strong effect on the synchrotron radio morphology of a young remnant. In 2-D MHD simulations of the SSDW phase of a young Type Ia ($`s=0`$) supernova remnant, Jun & Norman (1996) found that turbulent energy eventually rose to about 0.6% of the kinetic energy in the remnant, and that random magnetic energy rose as well, to a level of about 0.3% of the turbulent energy. The random magnetic-field energy did not track the total turbulent energy precisely, but spatial locations of strong turbulence (Kelvin-Helmholtz unstable edges of Rayleigh-Taylor fingers) were also the locations of strong magnetic-field amplification. In synchrotron visualizations, Jun & Norman (1996) also found that regions of strong magnetic field produced strong synchrotron emission, but their simulation used a relatively simple description of relativistic-electron acceleration and transport. We expect that magnetic fields in a bubble simulation would also track locations of turbulence and vorticity, that is, should indicate bubble walls, and ought to be noticeable in synchrotron images. In addition, the far higher levels of turbulent energy we find compared to smooth ejecta models such as Jun and Norman’s should result in considerably higher magnetic energy densities and consequent higher synchrotron emissivities. However, hydrodynamic simulations including magnetic-field tracking and particle acceleration and subsequent evolution will be necessary to make more definite predictions.

Our results are in qualitative agreement with the observed morphology of Cas A, undoubtedly the best example of a young SNR with ample evidence for vigorous mixing. Our simulations are able to reproduce filamentary ejecta emission seen in Chandra images. The ragged appearance of Cas A, including its “jet” feature, may also be a consequence of an interaction of bubbly ejecta with the ambient medium, and not necessarily the result of a strongly asymmetric explosion. It is also tempting to identify its rings of optical knots (Reed et al. 1995) with the most dense ejecta clumps at the boundaries of Fe bubbles. Our current simulations are however not suitable for a detailed modeling of the Cas A dynamics because of the evidence for the dynamically important shell of circumstellar matter in Cas A Chevalier & Liang (1989). The circumstellar interaction, studied by us through one-dimensional hydrodynamical simulations Borkowski et al. (1996), must be now simulated in 3-D in the framework of bubbly ejecta. A quantitative study of spatial morphology of Cas A, coupled with a more detailed examination of spatial structures in hydrodynamical data sets, would also be useful.

X-ray observations should provide a most complete test of the two-phase model, because they probe the bulk of the shocked material. While modeling of X-ray emission based on multidimensional calculations clearly demands a separate effort, we briefly outline what major effects are expected. If most Fe initially resides in low-density shocked bubbles, while other heavy elements are located in the dense phase of SN ejecta, then Fe lines should generally be much weaker than lines from other abundant heavy elements. Fe lines indeed appear to be weaker than expected in Cas A Borkowski et al. (1996); Vink et al. (1996); Favata et al. (1997) when compared with strong Si, S, Ar, and Ca lines, based on nucleosynthetic yields of core-collapse SNe. Unlike in the optical and IR, freshly synthesized Fe was detected by Chandra Hughes et al. (2000); Hwang et al. (2000), but its spatial distribution is different than that of Si- and S-rich ejecta. This means that not all Fe resided in low-density bubbles, either because of microscopic mixing of Fe at the bubbles’ boundaries or because radioactive Ni was mixed to sufficiently high velocities where the bubble growth was inhibited by escape of $`\gamma `$-rays from Ni clumps. Quantitative analysis of Chandra Cas A data, coupled with hydrodynamical simulations and X-ray modeling, should provide a stringent test of the two-phase ejecta model with Fe bubbles.

As remnants become older, turbulent mixing will lead to a gradual strengthening of Fe lines with respect to other heavy elements. Even if Fe is microscopically mixed with other elements during the SNR evolution, its ionization age should lag behind that of other elements. In the next few years, a combination of new X-ray observations and sophisticated hydrodynamical and X-ray modeling should provide us with a vastly better understanding of how chemical elements are ejected in explosions of massive stars and how they are dispersed into ambient interstellar medium.

We thank Dick McCray for discussions about Fe-Ni bubbles in SN 1987A. The three-dimensional simulations reported here were performed at the North Carolina Supercomputing Center using 100 processors of an IBM SP2. We thank NCSC and IBM for their generous support of computing resources. Support for this work was provided by NASA under grant NAG-7153.


format: text

text:
# Rare Kaon Decays

## Introduction

The study of rare kaon decays has a glorious history. At the level of $`10^1`$ in branching ratio, parity violation was first “discovered” as the $`\theta `$-$`\tau `$ puzzleleeyang , CP violation was discovered at the level of $`10^3`$firstcp , and the GIM mechanismgim  was suggested by the strong suppression of the decay $`K_L^0\mu ^+\mu ^{}`$, or the absence of the Flavor Changing Neutral Currents (FCNC)— the decay mode was subsequently found at the level of $`10^8`$. Further searches for rare kaon decays at lower branching ratio undoubtably involve rich physics and may yet offer another surprise— e.g. they are likely to elucidate the origin of CP violation, and provide some clues for physics beyond the Standard Model (SM). In this talk, the decay modes via FCNC with the expected branching ratios of less than $`10^9`$ are discussed. Since exotic decays violating the lepton flavor conservation law have been discussed by the previous speakerklaus , the focus of this talk is on the rare kaon decays which are allowed in the SM.

Typical leading Feynman diagrams of these decays are shown in Fig. 1, where $`d`$ and $`u`$ at the top lines are for $`K^0`$ and $`K^+`$ decays, respectively, and $`\mathrm{}`$ indicates the electron or the muon. Because of the GIM mechanism, the lowest order diagrams come from second order weak interactions with a $`u`$-type quark in the loop diagrams, in which the top-quark contributions dominate because of the mass. This makes these decay modes very sensitive to $`V_{td}`$, the least constrained coupling-constant of the top and down quarks in the Cabibbo-Kobayashi-Maskawa (CKM) matrix. Due to one of the unitarity conditions of the CKM matrix $`V_{ud}V_{ub}^{}+V_{cd}V_{cb}^{}+V_{td}V_{tb}^{}=0`$ and with the approximation $`V_{ud}V_{tb}1`$, $`V_{td}`$ and $`V_{ub}^{}`$ form two sides of a unitarity triangle as shown in Fig. 2. The height of the triangle $`Im(V_{td})`$ is an indication of “direct” CP violation, or CP violation through the decay amplitude. For the decay $`K^+\pi ^+\nu \overline{\nu }`$, the branching ratio is roughly proportional to $`|V_{td}|^2`$ and for $`K_L^0\pi ^0\nu \overline{\nu }`$ to $`Im^2(V_{td})`$. These measurements alone can determine the unitarity triangle, and when they are combined with measurements from $`B`$-decays it will provide over-constrained information that is sensitive to a presence of supersymmetry and other physics beyond the SMbb ; nir .

## Decays $`K_L^0\pi ^0\mathrm{}^+\mathrm{}^{}`$ and $`K_L^0\mathrm{}^+\mathrm{}^{}`$

If the final states include a charged lepton pair, there are additional contributions from diagrams with virtual photons that usually prevent clear interpretations of measurementsvalencia .

The decays $`K_L^0\pi ^0\mathrm{}^+\mathrm{}^{}`$ have a “direct” CP violating component, which is expected to occur, if it were the only component, at a branching ratio of $`5\times 10^{12}`$. There are two other contributions, however, at the same level to this process; one arises from the mixing of the CP even state in $`K_L^0`$, and the other from two-virtual-photon intermediate states that conserve CP. At present, there are theoretical ambiguities in the estimations of these contributions, which may eventually be sorted out by measurements of $`K_S^0\pi ^0\mathrm{}^+\mathrm{}^{}`$ and other radiative decays. To make the matter a bit more complicated, there is a physical background coming from the radiative decay $`K_L^0\mathrm{}^+\mathrm{}^{}\gamma \gamma `$, which has the identical event topology. With the tightest cuts, the background level is estimated to be still at a $`10^{11}`$ level for $`\mathrm{}=e`$greenlee . The situation is similar in the case of $`\mathrm{}=\mu `$. The present upper limits (90 % c.l.) are $`B(K_L^0\pi ^0e^+e^{})5.6\times 10^{10}`$pee  and $`B(K_L^0\pi ^0\mu ^+\mu ^{})3.8\times 10^{10}`$pmm .

The decay $`K_L^0e^+e^{}`$ is very similar to $`K_L^0\mu ^+\mu ^{}`$, which is sensitive to the $`\rho `$ parameterlw , but it is further suppressed by the helicity mechanism. This decay is sensitive to pseudo-scalar interactions coming from physics beyond the SM. Four events from the decay $`K_L^0e^+e^{}`$ have been reported recently which correspond to a branching ratio, $`B(K_L^0e^{}e^+)=8.7_{4.1}^{+5.7}\times 10^{12}`$ee , being consistent with the unitarity bound expected from the long-distance contributions.

## Decay $`K^+\pi ^+\nu \overline{\nu }`$

In the SM calculation of the decay $`K^+\pi ^+\nu \overline{\nu }`$, the dominant contribution comes from second order loop diagrams with a virtual top quark (Fig. 1). The hadronic matrix element can be extracted from the decay $`K\pi ^0e^+\nu `$ and theoretical uncertainties in the calculation due to long distance contributions and other effects are smallbb . The signature of the decay $`K^+\pi ^+\nu \overline{\nu }`$ is a single pion with no other observable decay-products. Definitive observation of this signal requires that all possible backgrounds are suppressed well below the signal level. Major background sources are: a muon from the copious decay $`K^+\mu ^+\nu `$ ($`K_{\mu 2}`$) which is misidentified as a pion, a pion from the decay $`K^+\pi ^+\pi ^0`$ ($`K_{\pi 2}`$) when two photons from the $`\pi ^0`$ decay are unobserved, a beam pion scattered by the target into the detector, and charge exchange reactions of $`K^+`$’s which result in decays $`K_L^0\pi ^+\mathrm{}^{}\overline{\nu }`$, where $`\mathrm{}`$ ($`e`$ or $`\mu `$) is undetected. In order to suppress the backgrounds, good particle identification, efficient photon veto, and detection of incoming particles are essential.

The E787 experiment at Brookhaven National Laboratory (BNL) as shown in Fig. 3 is designed to effectively distinguish these backgrounds from the signal. Kaons of about 700 MeV/c at a rate of $`(47)\times 10^6`$ per 1.6-s spill are detected and identified by a Čerenkov counter and hodoscopes, degraded by BeO and stopped in an active target, primarily consisting of 413 5-mm square scintillating fibers. The momentum ($`P`$), kinetic energy ($`E`$) and range ($`R`$) of decay products are measured using the target, a central drift chamber, 21 layers of 1.9-cm thick plastic scintillator (Range Stack) and two layers of straw chambers, all contained in a 1-T magnetic field. The $`\pi ^+\mu ^+e^+`$ decay sequence of the decay products in the Range Stack scintillator is observed by 500-MHz transient digitizers for particle identification. Photons are detected by a $`4\pi `$-sr calorimeter consisting of a 14-radiation-length-thick lead/scintillator barrel detector, 13.4-radiation-length-thick end caps of CsI crystals, and a 3.5-radiation-length-thick lead-glass Čerenkov counter which also works as an active beam degrader.

The E787 experiment reported an observation of one clean event at a branching ratio $`B(K^+\pi ^+\nu \overline{\nu })=4.2_{3.5}^{+9.7}\times 10^{10}`$ using the data sample taken in 1995e787a . Since then, additional data samples taken in 1996 and 1997, together with those taken in 1995, have been analyzed with improved algorithms, which have resulted in less non-Gaussian tails in $`P`$, $`R`$ and $`E`$ measurements, and a $``$30 % higher acceptance than that in Ref.e787a . Also, in the 1996–7 runs, lowering the incident $`K^+`$ beam momentum resulted in a larger fraction of kaons stopping in the target, which reduced accidental hits originated from nuclear reactions in the beam degrader. The higher proton intensity at the production target compensated the kaon yield loss at lower momentum.

In order to avoid a possible bias in the analysis, the signal region is kept untouched until the background estimations as well as optimization of acceptance have been done. In the background study, the data sample after applying all the final cuts except two orthogonal (uncorrelated) cut groups to be studied, e.g. the kinematical cuts and those related to the observation of the decay sequence $`\pi ^+\mu ^+e^+`$ in the stopping counter for the estimation of the $`K^+\mu ^+\nu `$ background, is used to obtain the suppression factor for each cut group. The correlation between the two cut groups, which may invalidate the method, is studied by varying the cuts being studied or by enhancing certain types of background events.

The total background for the entire 1995–1997 exposure with the final analysis cuts is estimated to be $`0.08\pm 0.02`$ events. The acceptance for $`K^+\pi ^+\nu \overline{\nu }`$, $`A=0.0021\pm 0.0001(stat)\pm 0.0002(syst)`$ is calculated based on data and Monte Carlo calculations. The largest uncertainty comes from the uncertainty in pion-nucleus interaction. The measurement of the branching ratio for $`K^+\pi ^+\pi ^0`$ within a few % of Ref.pdg  confirms the acceptance calculation. Analysis of the full data sample as shown in Fig. 4 has yielded only the single event previously reported. Based on the acceptance $`A`$ and the total exposure of $`N_{K^+}=3.2\times 10^{12}`$ kaons, the new branching ratio is $`B(K^+\pi ^+\nu \overline{\nu })=1.5_{1.2}^{+3.4}\times 10^{10}`$new . This provides a constraint, $`0.002<|V_{td}|<0.04`$.

The goal of a new experiment E949 at BNL is to improve the sensitivity by an order of magnitude. Since the AGS in the RHIC era will be used for two hours a day to feed heavy ions into the RHIC ring, the rest of 22 hours can be used for the high energy program. The operation is expected to provide a more stable and longer running period. Exploiting a higher beam intensity, the incident $`K^+`$ momentum can be further lowered to reduce accidental coincidence. The photon veto capability will be improved by additional lead/scintillator layers in the barrel region and by additional active degrader in the beam region. In the region below the $`K^+\pi ^+\pi ^0`$ peak at 205 MeV/c, where nuclear interactions result in a large momentum-tail, the additional photon veto capability may suppress the background by more than an order of magnitude as low as to the signal level, doubling the phase space in the search.

### Decay $`K_L^0\pi ^0\nu \overline{\nu }`$

The decay $`K_L^0\pi ^0\nu \overline{\nu }`$ violates the CP conservation law through decay amplitude. This process is expected to occur at $`B(K_L^0\pi ^0\nu \overline{\nu })3\times 10^{11}`$ bb . The contribution from CP mixing in $`K_L^0`$ is expected to be around $`10^{15}`$ll . Since no charged leptons are involved in this decay, the process is free from virtual photon contributions and clean for the study of the origin of CP violation. The theoretical ambiguity is only $`1\%`$ except those in the CKM matrix elements. Conversely, the observation of this decay mode uniquely determines $`Im(V_{td})`$ or $`\eta `$ in the Wolfenstein parametrization. The goal of the experiment is to determine $`Im(V_{td})`$ with a 10–15 % accuracy, which corresponds to a single event sensitivity of $`10^{12}`$. The present upper limit of this decay is $`5.7\times 10^7`$ktev .

The decay $`K_L^0\pi ^0\nu \overline{\nu }`$ is a three-body decay that involves only neutral particles. The signature of this decay is two photons from the $`\pi ^0`$ and no other activity in the detector. In this decay mode, available kinematical parameters are very limited; relatively easy ones to measure are the positions and energies of $`\gamma `$-rays. Also, the direction and position of the $`K_L^0`$ beam can be limited by tightly collimating the beam at a cost of beam intensity. The decay vertex can be reconstructed from these two constraints with the assumption of the pion mass. The KEK experimentinagaki  and the FNAL approachkami  are classified in this category.

The BNL experiment E926e926  attempts to measure more kinematic values; the directions of $`\gamma `$-rays and the momentum of the $`K_L^0`$. Measurements of $`\gamma `$-ray directions allow full reconstruction of the $`\pi ^0`$ kinematics without the assumption of the $`\pi ^0`$ mass, providing more constraints with redundancy, which is necessary to suppress the background to the level well below the signal. The time-of-flight (TOF) of a $`K_L^0`$ between the production target and the decay vertex can be measured for low momentum $`K_L^0`$’s if the incident proton beam is bunched. This measurement allows calculation of the missing mass and kinematical reconstruction in the center of mass system, which is effective to eliminate backgrounds from two-body decays. The major decay modes of $`K_L^0`$ are $`K_L^0\pi ^0\pi \pi `$ and $`K_L^0\pi ^\pm \mathrm{}^{}\nu `$. Suppression of most backgrounds is achieved by high-efficiency hermetic photon and charged-particle detector system surrounding the decay volume, and kinematical constraints.

In the BNL E926 experiment, a low energy $`K_L^0`$ beam with an average momentum around 650 MeV/c is produced by irradiating a target with a 24-GeV proton beam from the AGS and extracted at 40o with respect to the incident beam. The proton beam is bunched to form a $``$200 ps wide bucket at a rate of 25 MHz. About 16 % of $`K_L^0`$’s decay in the 4-m long decay volume, which is evacuated to a level of $`10^7`$ Torr to suppress the background from neutron-induced $`\pi ^0`$ production. The decay region is surrounded by a charged particle veto system and a photon veto system of 18-radiation-length-thick lead/scintillator sandwiches in the barrel region. The two photons from the $`\pi ^0`$ decay are converted into pairs of a positron and an electron in a 2-radiation-length preradiator next to the vacuum region for the measurement of the directions of the $`\gamma `$-rays. The preradiator consists of sandwiches of 2-mm thick scintillator, a copper plate as a mechanical support and radiator, and a tracking chamber. This is followed by an 18-radiation-length calorimeter for the measurement of $`\gamma `$-ray energies and for vetoing additional photons.

The estimates of sensitivity for $`K_L^0\pi ^0\nu \overline{\nu }`$ are tightly coupled to the cuts required for background suppression, particularly for the $`K_L^0\pi ^0\pi ^0`$ and $`K_L^0\pi ^0\pi ^+\pi ^{}`$ backgrounds. An acceptance of $`0.015`$ % for the case S/N=0.5 comes from the combination of factors; 0.58 for the fiducial region and usable kaon momentum region, 0.33 for the solid angle, 0.5 for the efficiency of the preradiator, and the remaining factor for the $`\pi ^0`$ mass cut and other cuts. Assuming three years of running with $`6\times 10^6`$ kaons per 2-s beam spill (50 % duty factor), the expected number of events is 65.

The toughest background is the CP violating $`K_L^0\pi ^0\pi ^0`$ decay when two photons are undetected. Backgrounds from $`K_L^0\pi ^0\pi ^0`$ arise when two photons are detected in the forward detector and the other two are undetected anywhere. These backgrounds can be classified into two categories; even pairing when two photons come from the same $`\pi ^0`$, and odd pairing when two photons come from different $`\pi ^0`$’s. Events in the even pairing category form a two-body-decay peak in the momentum spectrum of $`\pi ^0`$ in the CM system, while the two $`\gamma `$-rays in the odd pairing category do not reproduce the $`\pi ^0`$ mass. The above requirements and photon veto essentially suppress $`K_L^0\pi ^0\pi ^0`$ backgrounds to 0.2 of the signal level. Similarly, $`K_L^0\pi ^0\pi ^+\pi ^{}`$ backgrounds can be suppressed to 0.1 of the expected signal. Backgrounds from other $`K_L^0`$ decay modes are estimated to be less than 0.1 of the signal. Because of the large angle extraction, the cross sections for producing $`\mathrm{\Lambda }`$’s are small and they completely decay before reaching the decay volume. Backgrounds could arise from $`\mathrm{\Lambda }`$’s produced by halo neutrons and $`K_L^0`$’s, but they are estimated to be negligible. Neutrons with $`P_n800`$ MeV/c can react on the remaining atoms in the vacuum region to produce $`\pi ^0`$’s. This background level is again estimated to be 0.01. Accidental backgrounds are caused by beam halo neutrons and $`\gamma `$-rays which create a $`\pi ^0`$ signal in the preradiator. The background level is estimated to be 0.02 of the expected signal.

The BNL E926 was proposed in 1997 but the group is still waiting for full funding. In the present scenario, the detector construction is expected to start in 2002. The KEK experiment has partially been approved, but the sensitivity goal is only around the SM level—the new Japan Hadron Facility is expected to improve the sensitivity.

## Conclusion

Studies of rare kaon decays have contributed to the discoveries of several symmetry violations. The decays $`K^+\pi ^+\nu \overline{\nu }`$ and $`K_L^0\pi ^0\nu \overline{\nu }`$ are expected to “complete” the measurement of the CKM matrix in the next decade independently from the $`B`$-decay system, and to elucidate the origin of CP violation.


text:
# Monte Carlo Simulations of Thermal-Nonthermal Radiation from a Neutron Star Magnetospheric Accretion Shell

## 1 Introduction

The high energy radiation from compact astrophysical objects is emitted by relativistic or semi-relativistic thermal and nonthermal leptons (electrons and pairs) via synchrotron, bremsstrahlung, and Compton processes, plus bound-bound and bound-free transitions of high-Z elements. Since Compton scattering is a dominant radiation mechanism in this regime, the most efficient and accurate method to model the transport of high energy radiation is the Monte Carlo (MC) technique. During the past decade we have developed a versatile state-of-the-art space-and-time-dependent MC code to model the radiative output of compact objects (see e.g. Liang et al. 2000). Recently we have added the self-consistent evolution of the leptons using a Fokker-Planck scheme. The lepton evolution is then coupled to the photon transport. Since this is the first time that we report on results obtained with this code, the first part of the present paper is devoted to a detailed description of the capabilities of the code (§2) and its verification in comparison with previous work (§3).

In the second part of this paper (§4), we apply this coupled MC-FP code to the study of the reprocessing of blackbody radiation from a magnetized neutron star surface by a magnetized spherical shell at the Alfvén radius. Recent observations of weak-field neutron star binaries, such as low-luminosity, X-ray bursters (atoll sources, e.g., 4U 1608-522: Zhang et al. (1996), 4U 1705-44: Barret et al. (1996)), bursting soft X-ray transients (e.g., Aql X-1: Harmon et al. (1996)), or pulsar binary systems (e.g., PSR B1259-63: Tavani et al. (1996)) indicate that many of them exhibit soft (photon index $`2`$) power law tails extending beyond $`100`$ keV, at least episodically, in addition to the thermal component at temperatures of $``$ a few keV, which presumably originates from the stellar surface. The luminosity of this high-energy tail appears to be anti-correlated with the soft X-ray luminosity (Barret & Vedrenne (1994); Tavani & Liang (1996)).

The origin of this high energy tail is unexplained at present. It could be due to thermal Comptonization by a hot coronal plasma, or it could be due to nonthermal emission. Tavani and Liang (1996) examined systematically the possible sites of nonthermal emissions and concluded that the Alfvén surface is the most likely candidate since the dissipation of the rotation energy of the disk is strongest there, due to magnetic reconnection and wave turbulence generation. Here we first focus on particle acceleration by wave turbulence.

Because of the 1-D nature of our code at this stage, the neutron star is assumed to emit isotropic blackbody radiation, the reprocessing magnetospheric plasma is idealized as part of a spherical shell at the Alfvén radius, and the magnetic field is taken to be nondirectional, so that the radiation output is isotropic. Even though this is not a perfect representation of a quasi-dipolar magnetosphere and accretion flow, we believe that, except for very special viewing angles, our output spectra should be a reasonable approximation to the angle-averaged output of the reprocessing by a hybrid thermal/nonthermal magnetosphere. We assume that the leptons are energized by Coulomb collisions with virial ions and accelerated nonthermally by Alfvén and whistler wave turbulence, and cooled by cyclotron/synchrotron, bremsstrahlung, and inverse Comptonization of both internal soft photons and blackbody photons from the stellar surface.

The primary focus of the parameter study presented in §4 is the application to weakly magnetized neutron stars with surface magnetic fields of $`B_{\mathrm{surf}}10^{11}`$ G. In particular, we will show that the anti-correlation of the hardness and luminosity of the hard X-ray emission with the soft X-ray luminosity is a natural consequence of the energetics of particle acceleration and cooling near the Alfvén radius. We predict that the nonthermal tails in the hard X-ray spectra of accreting, weakly magnetized neutron stars may extend up to $`1`$ MeV. A solid detection and the measurement of the cutoff energy of these high-energy tails by the INTEGRAL mission, scheduled for launch in 2002, will provide important constraints on accretion-based models for the hard X-ray emission from accreting neutron stars.

In this context, it is interesting to note that Chatterjee, Hernquist, & Narayan (2000; see also Mereghetti & Stella (1995); Wang (1997); Chatterjee & Hernquist (2000)) have recently proposed a similar type of accretion-powered emission for anomalous X-ray pulsars (AXPs), as an alternative to models based on magnetic-field decay (Thompson & Duncan (1996)) or residual thermal energy (Heyl & Hernquist (1997)). According to Chatterjee et al. (2000) the X-ray emission from AXPs (which generally consists of a soft, thermal component with $`kT0.3`$ – 0.4 keV plus a hard X-ray tail with photon index $`\mathrm{\Gamma }3`$ – 4) is powered by accretion of material from the debris of the supernova which had formed the neutron star, onto a the surface of the neutron star, which possesses a typical pulsar magnetic field of $`B_{\mathrm{surf}}10^{12}`$ G. Therefore, we extend our parameter study to parameter values relevant to accreting pulsars. However, we point out that in the case of a magnetic field as high as $`B_{\mathrm{surf}}10^{12}`$ G, the assumed shell geometry and the quasi-isotropy of the emission from the neutron star surface may be a gross over-simplification. However, although consequently the precise parameter values used in this region of the parameter space should not be taken at face value, our parameter study might still provide interesting insight into the dependence of the equilibrium electron and photon spectra on the various input parameters in the high-magnetic-field case.

## 2 Physics of photon and lepton evolution

We use the Monte Carlo (MC) technique (Podznyakov, Sobol, & Sunyaev 1983; Canfield, Howard, & Liang 1987; Liang (1993); Hua, Kazanas, & Titarchuk 1997) to simulate relativistic photon transport. We include the full (energy-and-angle-dependent) Klein-Nishina cross section for Comptonization, relativistic bremsstrahlung from lepton-ion and lepton-lepton scattering (Dermer (1984)), and cyclo-synchrotron processes (Brainerd (1984)) in magnetic fields. The MC photon transport is fully space-and-time-dependent. Photons are born with a certain “weight” (which is basically the total energy in photons within a spatial region or emitted from a photon source, divided by the number of Monte-Carlo particles) which is diminished by absorption and escape, until the weight drops below a user-specified limit, at which point the photon is ”killed”. In the simulations shown later in the paper, we specify this energy weight cutoff as 1/100 of the initial statistical weight of the photon. The final results are insensitive to the choice of this energy weight cutoff, as long as it is $`1`$. Surviving photons are sampled at boundaries to provide time-and-frequency-dependent spectral output. In addition to self-emitted photons from the plasma, soft photons can be injected at zone boundaries and inside volume elements. Currently the code can handle 1-D spherical, cylindrical or slab geometries with an arbitrary number of spatial zones. However since the photon ray tracing is done with full angle informations, the generalization to 2-and-3-D transport is straight forward. The maximum number of photon frequency bins is 128. For more details of this code see Canfield et al. (1987) and Böttcher & Liang (1998). A typical MC run with a million particles at a Thomson depth of a few takes 10s of minutes on a DEC alpha server. Since the CPU time usage scales as the square of the Thomson depth ( number of scatterings), Thomson thick runs can be quite time consuming. We are currently developing a random walk scheme for Wien photons trapped in Thomson thick zones, which would save large amounts of CPU time without introducing too much error.

The lepton population is computed locally in each spatial zone using the Fokker-Planck approximation (Dermer, Miller, & Li 1996; Li, Kusunose & Liang 1996a), taking into account coulomb and Moller scattering, stochastic acceleration by Alfvén and whistler wave turbulence, and radiative cooling (plus pair processes if necessary). In general the lepton population consists of a low energy thermal population plus a high-energy tail truncated at the highest energies by radiative cooling. The photon and lepton evolutions are coupled to each other via a quasi-implicit time scheme in which we use an average of the photon-electron energy exchange rates between two subsequent time steps to compute the electron cooling rates, in particular due to Compton scattering. The Fokker-Planck equation governing the electron evolution is solved using a fully implicit scheme. In contrast to several other codes currently used in the literature (e.g., Stern et al. (1995); Li et al. 1996a , Malzac & Jourdain (2000)), we solve the evolution of the entire electron population with the FP scheme and do not introduce any artificial separation between thermal and non-thermal particles. Since the lepton distribution typically evolves much faster than the photon distribution, each photon cycle contains many lepton cycles. The user, however, can always turn off the nonthermal lepton acceleration and assume a strict thermal population whose temperature can be computed self-consistently from energy balance alone.

In oder to calculate the emissivities and opacities for thermal cyclotron, non-thermal synchrotron, and thermal bremsstrahlung emission (and only for this purpose), the electron distribution calculated with our Fokker-Planck scheme is decomposed into a thermal population plus a non-thermal tail. For the thermal bremsstrahlung and non-thermal synchrotron emission and absorption, the standard expressions from Rybicki & Lightman (1979) are used. For the thermal cyclotron emission, we add explicitly over the first 5 harmonics, beyond which we use the asymptotic continuum representation of Mahadevan, Narayan & Yi (1996).

Since in many cases of interest to the current investigation the moderately to strongly magnetized coronal plasma is optically thick at low frequencies, $`h\nu 1`$ keV, due to synchrotron-(self-)absorption, we use the following simplification in order to save CPU time: For any given frequency $`h\nu 1`$ keV, we compare the absorption length $`l_{\mathrm{abs}}^\nu `$ to the Compton mean free path, $`l_{\mathrm{Compt}}^\nu `$ and the radial extent $`\mathrm{\Delta }R`$ of the current zone of the Comptonizing region. If

$$l_{\mathrm{abs}}^\nu <0.1\mathrm{min}\{l_{\mathrm{Compt}}^\nu ,\mathrm{\Delta }R\}$$
(1)
any soft photon produced at this frequency has a very small probability of escaping the current zone or being up-scattered (to a frequency at which the absorption length will be different so that the above criterion has to be re-evaluated after the scattering event) before being re-absorbed. Thus, in the frequency ranges where Eq. (1) is fullfilled, the radiation escaping at the boundary of that zone will be approximately given by the respective section of the thermal blackbody spectrum. Accordingly, our code neglects the volume emissivities in those frequency ranges and, instead, produces thermal blackbody photons at the zone boundaries.

During each photon time step, the code keeps track of the energy transferred between photons and electrons due to Compton scattering and cyclotron/synchrotron and bremsstrahlung emission and absorption. The respective heating and cooling rates are used to scale the energy loss/gain coefficients of electrons at a given energy $`E_e=\gamma m_ec^2`$. For the bremsstrahlung energy loss rate, we use the approximate scaling law $`(d\gamma /dt)_{\mathrm{br}}\gamma ^{1.1}`$ (Böttcher, Pohl, & Schlickeiser 1999). Coulomb heating/cooling is included using the energy exchange and dispersion rates given in Dermer & Liang (1989). We assume that the ions have a pre-specified temperature $`kT_p1`$ MeV, which is not significantly affected by any changes of the electron temperature. For the energy exchange due to Møller scattering (elastic electron-electron scattering) we use the electron energy exchange and dispersion coefficients given in Nayakshin & Melia (1998).

In addition to Coulomb interactions and radiative losses, we also account for stochastic (2nd order Fermi) acceleration due to Alfvén and Whistler waves. For a given background magnetic field $`B`$, the magnitude and spectrum of hydromagnetic plasma wave turbulences are determined by the parameters $`\delta ^2(\mathrm{\Delta }B/B)^2`$, where $`\mathrm{\Delta }B`$ is the amplitude of the magnetic-field fluctuations, and $`q`$, the spectral index of the turbulence spectrum. We will generally use $`q=5/3`$, characteristic of Kolmogorov turbulence. The electron acceleration and energy dispersion rates are calculated using the formalism of Schlickeiser (1989). However, we have to take into account that those plasma waves interacting efficiently with the low-energy, quasi-thermal part of the electron spectrum will be efficiently damped and absorbed by the process of Landau damping (e.g., Schlickeiser, Fichtner & Kneller 1996). This leads to a strong truncation of the Kolmogorov wave spectrum above a critical wave number for which the equivalent absorption depth through the region occupied by the hot plasma exceeds unity. We take the effect of Landau damping into account by introducing an effective absorption depth $`\tau _k`$ at wave number $`k`$, and correcting the acceleration rate $`\dot{\gamma }_A^0`$ of Schlickeiser (1989) due to the “optically thin” plasma wave spectrum by an absorption term:

$$\dot{\gamma }_A=\dot{\gamma }_A^0\frac{\left(1e^{\tau _{k_{\mathrm{res}}}}\right)}{\tau _{k_{\mathrm{res}}}}.$$
(2)
Here $`\tau _{k_{\mathrm{res}}}=\mathrm{\Gamma }_{k_{\mathrm{res}}}t_\mathrm{A}`$, and $`\mathrm{\Gamma }_k`$ is the Landau damping rate at wave number $`k_{\mathrm{res}}`$ at which Alfvén waves are resonating preferentially with electrons of energy $`\gamma `$. $`t_A=\mathrm{\Delta }R/v_A`$ is the Alfvén crossing time of the zone. This modified acceleration term can be renormalized either to correspond to the pre-specified value of $`\delta ^2`$ for electrons resonating with waves in the weakly damped part of the Alfvén wave spectrum (i.e. at high electron energies), or by specifying a heat input rate into the electron ensemble due to resonant wave/particle interactions.

As an option, our code can account for pair production and annihilation and the corresponding photon absorption and emission. The pair annihilation rates and annihilation radiation emissivities are taken from Svensson (1982), and the $`\gamma \gamma `$ pair production rate of Böttcher & Schlickeiser (1997) is used. In addition, simple Compton reflection schemes, using the Green’s functions for reflection off neutral disk material of White, Lightman & Zdziarski (1988) and Lightman & White (1988), can be used to simulate a Compton reflection component either off the inner boundary (corresponding to a quasi-homogeneous slab geometry) or reflecting part of the radiation escaping at the outer boundary (e.g., corresponding to a cold outer disk).

Assuming local isotropy of the electron distribution, we solve the one-dimensional Fokker-Planck equation

$$\frac{n_e(\gamma ,t)}{t}=\frac{}{\gamma }\left[n_e(\gamma ,t)\frac{d\gamma }{dt}\right]+\frac{1}{2}\frac{^2}{\gamma ^2}\left[n_e(\gamma ,t)D(\gamma ,t)\right],$$
(3)
where $`D`$ is the energy dispersion coefficient. To solve Eq. 3, we use an implicit version of the discretization scheme proposed by Nayakshin & Melia (1998). We choose a logarithmic spacing in electron kinetic energy, $`x_i\gamma _i1`$. In the following the subscript $`i`$ refers to electron energy, while the superscript $`n`$ refers to time. We define $`f_i^nn_e(\gamma _i,t^n)/n_e`$ The discritization is then given by

$$f_i^n=f_i^{n+1}+\mathrm{\Delta }t\left(\frac{a_{i+1}^{n+1}f_{i+1}^{n+1}a_{i1}^{n+1}f_{i1}^{n+1}}{\mathrm{\Delta }x_i+\mathrm{\Delta }x_{i1}}\right)\mathrm{\Delta }t\left(\frac{\alpha D_{i+1}^{n+1}f_{i+1}^{n+1}2D_i^{n+1}f_i^{n+1}+\rho D_{i1}^{n+1}f_{i1}^{n+1}}{\alpha \mathrm{\Delta }x_i(\mathrm{\Delta }x_{i1}+\mathrm{\Delta }x_i)}\right),$$
(4)
where $`a_i=(d\gamma /dt)(x_i)`$, $`D_i=[d(\mathrm{\Delta }\gamma )^2/dt](x_i)`$, $`\alpha =2/(1+q)`$, $`\rho =2q/(1+q)`$, with $`q=x_{i+1}/x_i`$. The system of equations is supplemented by the boundary conditions specified in Appendix A of Nayakshin & Melia (1998). The system of equations (4) is in tridiagonal form and can be easily solved to find the electron distribution at time $`t^{n+1}`$. In order to evaluate the energy exchange and dispersion coefficients $`a_i^{n+1}`$ and $`D_i^{n+1}`$ the code performs an energy balance calculation to find the average electron energy in the ensemble after the current time step, which is then used to calculate the energy exchange and dispersion coefficients due to Møller scattering. All other coefficients evolve sufficiently slowly (i.e. on time scales much longer than the electron-evolution time scale) so that the coefficients evaluated under the conditions at the beginning of the current time step can be used.

The implicit scheme of Eq. 4 is known to approach the equilibrium solution for the electron distribution exactly, although the temporal evolution calculated this way becomes inaccurate on short time scales. Since in the problems of interest here, the electron distributions typically evolve on timescales much shorter than the photon distributions, the electrons always attain a distribution close to local equilibrium during each photon time step. Hence, the degree of accuracy provided by the implicit scheme is sufficient for our purposes. As an option, the code can perform the same simulations using an explicit scheme to solve the Fokker-Planck equation. In all test cases, both schemes produced virtually identical results, but the implicit scheme executes a factor of $`10`$ – 100 faster because of the larger individual Fokker-Planck time steps allowed in this scheme.

However, we did not find an appropriate method to calculate the coefficients for the pair production and annihilation rates at time step $`n+1`$. For this reason, we have to correct the electron and positron distribution functions for pair production and annihilation in an explicit manner. This becomes obviously inaccurate in the case of strongly pair producing model situations. Thus, our scheme is a significant improvement over existing schemes only for pair deficient situations with pair fractions $`f_{\mathrm{pair}}n_{e^+}/n_p`$ a few %. For strongly pair dominated situations, we would have to adjust the individual time steps of our simulations so far that such simulations would take a prohibitive amount of computing time.

In the present paper, we focus on results for equilibrium situations, for which we let the code evolve until a stable electron distribution and photon output spectrum is reached. The convergence of the electron and photon spectra can be greatly accelerated if one starts out with an appropriate first guess for the equilibrium electron temperature. For this reason, we have developed an analytical estimate of the equilibrium electron temperature, which is described in Appendix A. This analytical estimate is implemented in our code and can be used to determine the appropriate initial conditions for our equilibrium simulations. Detailed tests and applications of the time-dependent features of the code will be presented in future publications.

## 3 Tests of the numerical scheme and comparison with Previous Results

For verification of our code, we have first compared the individual energy exchange and diffusion rates, $`d\gamma /dt`$ and $`D(\gamma )`$, with those obtained in earlier work. Our numerical energy exchange and diffusion rates due to Coulomb scattering are in agreement with those of Dermer & Liang (1998)222note that in their Eq. (A2), the last $`\gamma ^{}`$ has to be replaced by $`(\gamma ^{})^2`$. The numerical values of the Møller scattering energy exchange and diffusion coefficients were found in good agreement with those of Nayakshin & Melia (1998)333note that there is a ’-’ sign missing in front of the term $`(\gamma \gamma _1)^2`$ in their Eq. (35).

We have verified that in cases in which radiative cooling and hydromagnetic acceleration are inefficient, our Fokker-Planck scheme correctly produces a thermal electron distribution in temperature equilibrium with the protons.

In the analysis of coronal energy and radiation transfer, it is customary to paramatrize energy input and dissipation rates $`\dot{u}`$ \[ergs cm-3 s-1\] by the dimensionless compactness, $`ł`$, which in spherical geometry, is defined by $`ł_{\mathrm{sph}}4\pi \sigma _\mathrm{T}R^2\dot{u}/(3m_ec^3)`$, where $`R`$ is the radius of the spherical volume. In slab geometry, this becomes $`ł_{\mathrm{slab}}\sigma _\mathrm{T}H^2\dot{u}/(m_ec^3)`$, where $`H`$ is the thickness of the slab. In our simulations, we specify the heating mechanisms to be Coulomb heating through the thermal protons and resonant wave-particle interactions. Thus, leaving the shape of the electron spectrum general, we calculate the respective energy dissipation rates as

$$\dot{u}_{\mathrm{Coul}/\mathrm{A}}=m_ec^2\underset{1}{\overset{\mathrm{}}{}}𝑑\gamma n_e(\gamma )\left(\frac{d\gamma }{dt}\right)_{\mathrm{Coul}/\mathrm{A}},$$
(5)
where the subscript ’Coul/wp’ stands for Coulomb scattering and wave-particle interaction, respectively. Obviously, the resulting compactness will depend not only on the energy density in protons and magnetic turbulence, but also on the current electron distribution (i.e. temperature, if the distribution is predominantly thermal), and can thus not be specified a priori as a free parameter. For comparison with previous work, we need to find appropriate values of the ion temperature and the wave turbulence amplitude in order to reproduce the dissipation compactness quoted in those papers. We define the temperature corresponding to the electron spectrum resulting from our simulations by requiring that the average electron Lorentz factor, $`\gamma =(1/n_e)𝑑\gamma \gamma n_e(\gamma )`$, be equal to the average Lorentz factor of a thermal electron population of the respective temperature,

$$\gamma _{\mathrm{th}}=\frac{K_3\left(\frac{1}{\mathrm{\Theta }_e}\right)}{K_2\left(\frac{1}{\mathrm{\Theta }_e}\right)}\mathrm{\Theta }_e.$$
(6)
To test our energy balance calculations, we have performed a series of simulations to reproduce the results of the slab corona model of Dove et al. (1997), in particular the inset of their Fig. 2, with which we find good agreement. For those simulations, we specified ion temperatures of $`80`$ MeV $`kT_i250`$ MeV, a magnetic field in equipartition with the ion population, and a negligibly low level of Alfvén wave turbulence to reproduce various values of the local heating compactness $`ł_c`$. As mentioned earlier, our code is rather inefficient in simulating strongly pair-producing, high-$`l_c`$ situations because of the explicit scheme to solve for pair balance. For this reason, we restrict the applications of the code in its present version to model situations with pair fractions $``$ a few %.

Li et al. (1996a, b) have developed a code solving simultaneously the Fokker-Planck equations for both the electron and the photon distribtions in a homogeneous medium, including Coulomb interactions, thermal bremsstrahlung, Compton scattering, resonant wave-particle interactions, and pair production and annihilation. There are two major differences between their approach and ours: (1) We use a Monte-Carlo method to solve the photon transport, and (2) we solve the Fokker-Planck equation for the entire electron spectrum, while Li et al. (1996a, b) split the electron distribution up into a thermal “bath” plus a non-thermal tail, assuming a priori that electrons of energies $`\gamma <\gamma _{\mathrm{thr}}=1+4\mathrm{\Theta }_e`$ attain a thermal distribution, and that acceleration due to wave/particle interactions affects only particles beyond $`\gamma _{\mathrm{thr}}`$. The latter simplification is justified by the argument that at low electron energies the thermalization time scale is much shorter than any other relevant time scale, and that long-wavelength plasma waves, with which low-energy electrons resonate preferentially, are strongly damped and in energy elequilibrium with the thermal pool of electrons. In our simulations, both effects are taken into account self-consistently without making the a-priori assumption of the existence or development of a non-thermal population.

Li et al. (1996b) present two model calculations to explain hard tails observed in the hard X-ray / soft $`\gamma `$-ray spectra of Cyg X-1 and GRO J0422+32. A spherical region of radius $`R=1.2\times 10^8`$ cm is assumed. For the case of Cyg X-1, they specify a total heating compactness of $`ł=4.5`$, and a non-thermal heat input into suprathermal particles ($`\gamma >\gamma _{\mathrm{thr}}`$) of $`ł_{\mathrm{st}}/ł=0.15`$, wich a turbulence amplitude of $`\delta ^2=0.059`$. A soft blackbody radiation component at $`kT_\mathrm{s}=0.1`$ keV from the outer boundary of the sphere is assumed to provide a soft radiation compactness of $`ł_\mathrm{s}/ł=0.07`$. The seed Thomson depth of the region is $`\tau _\mathrm{p}=0.7`$. In their simulations, Li et al. (1996b) find a significant suprathermal tail in the resulting electron distribution, which leads to an excess hard X-ray / soft $`\gamma `$-ray emission, consistent with the observed one. The temperature of the thermal part of their electron ensemble is found at $`kT_\mathrm{e}=139`$ keV, and the equilibrium pair fraction is $`f_{\mathrm{pair}}1.8`$ %. For the case of GRO J0422+32 they specify $`ł=0`$, $`ł_{\mathrm{st}}/ł=0.06`$, and $`ł_\mathrm{s}/ł=0.08`$. This resulted in a weaker nonthermal electron tail, an equilibrium temperature of $`kT_\mathrm{e}=133`$ keV, a wave amplitude of $`\delta ^2=0.065`$, and a pair fraction of $`f_{\mathrm{pair}}3`$ %.

For comparison with our code, we have run simulations with the same total compactness $`ł`$, soft compactness $`ł_\mathrm{s}`$, soft blackbody temperature $`kT_{\mathrm{BB}}`$, radius $`R`$, seed Thomson depth $`\tau _p`$, and the same values of the plasma wave amplitude normalization $`\delta ^2`$. A difficulty in comparing the two codes is that in Li et al. (1996a, b) the ion temperature and magnetic field are not specified explicitly. In our comparative simulations we assume a magnetic field in equipartition with the ion energy density. The results of our simulations are illustrated in Figs. 1 and 2. Our results are qualitatively similar to those of Li et al. (1996a, b). However, we find somewhat lower equilibrium temperatures and stronger nonthermal electron tails as well as stronger suprathermal acceleration compactnesses $`ł_{\mathrm{st}}`$. The lower temperatures may be attributed to a rather prominent cyclotron/synchrotron cooling (comparable to Compton cooling) in our simulations (see Figs. 1b and 2b). This, in combination with the larger $`ł_{\mathrm{st}}`$ values, seems to indicate that in our simulations we have used higher magnetic field values than Li et al. (1996a, b). However, even abandoning the equipartition assumption, we could not find self-consistent parameter values resulting in the same combination of input parameters used by Li et al. (1996a, b). Given this descrepancy in the way of input parameter specification, the qualitative agreement between the results of the two codes, using very different numerical methods, is encouraging.

## 4 Models of accretion onto a magnetized neutron star

While in the first part of this paper we were describing the general features and the verification of our MC/FP code, we are now applying this code to model the electron dynamics and photon transport arising from models of accretion onto a magnetized neutron star. In both weakly magnetized neutron stars (atoll sources and soft X-ray transient neutron star binary systems) with $`B_{\mathrm{surf}}10^{11}`$ G and X-ray pulsars (including, possibly, anomalous X-ray pulsars) with $`B_{\mathrm{surf}}10^{12}`$ G, energy dissipation will be most efficient at the Alfvén radius, where the optically thick, geometrically thin outer accretion disk is disrupted and the dynamics of the accretion flow becomes dominated by the magnetic field. We idealize this region of efficient energy dissipation at the Alfvén radius of disk accretion onto a magnetized neutron star as part of a spherical shell whose distance $`r_0`$, magnetic field, and column thickness are fixed by the accretion rate (Ghosh & Lamb 1979a ,b,1990). The distance $`r_0`$ is given by

$$r_0=2\times 10^8f\mu _{30}^{4/7}ł_{}^{2/7}M_{}^{1/7}R_6^{2/7}\mathrm{cm}$$
(7)
where $`\mu _{30}=(\mathrm{neutron}\mathrm{star}\mathrm{magnetic}\mathrm{moment})/(10^{30}\mathrm{G}\mathrm{cm}^3)`$, $`ł_{}=L/L_{\mathrm{Edd}}`$, $`M_{}=M_{\mathrm{NS}}/M_{}`$, and $`R_6=R_{\mathrm{NS}}/(10^6\mathrm{cm})`$. For the current simulations, for definiteness, we fix the Ghosh-Lamb fudge parameter $`f=0.3`$, and set $`M_{}=R_6=1`$. Hence, the dipole magnetic field at the Alfvén radius is

$$B_0=4.2\times 10^6ł_{}^{6/7}\mu _{30}^{5/7}M_{}^{3/7}R_6^{6/7}f_{0.3}^3\mathrm{G},$$
(8)
where $`f_{0.3}=f/0.3`$. The virial ion temperature at $`r_0`$ is

$$kT_i=\frac{2}{3}\frac{GMm_H}{r_0}9.3\left(\frac{r_0}{10^7\mathrm{cm}}\right)^1M_{}\mathrm{MeV}.$$
(9)
The column density of the shell can be estimated using the poloidal accretion rate $`\dot{M}4\pi r_0\mathrm{\Delta }r_0n_iv_pm_H`$, where we assume that the poloidal velocity $`v_pv_{\mathrm{ff}}/2`$ with $`v_{\mathrm{ff}}`$ being the free-fall velocity. Hence, the radial Thomson depth of the shell is approximately:

$$\tau _T=\mathrm{\Delta }r_0n_i\sigma _T\frac{\dot{M}\sigma _T}{2\pi r_0v_{\mathrm{ff}}m_H}0.97ł_{}^{8/7}\mu _{30}^{2/7}M_{}^{4/7}R_6^{1/7}f_{0.3}^{1/2}.$$
(10)
The neutron star (taken to be a 10 km spherical surface) is assumed to emit a blackbody luminosity at the temperature $`kT_{\mathrm{BB}}=1.78ł_{}^{1/4}`$ keV. In addition the level of wave turbulence is specified by the dimensionless amplitude $`\delta ^2=(\mathrm{\Delta }B/B)^2`$ and the spectral index q. The minimum wavevector $`k_{min}`$ is set to $`2\pi /(\mathrm{\Delta }r_0)`$, where $`\mathrm{\Delta }r_00.1r_0`$ is the shell thickness (Ghosh & Lamb 1979a ,b, 1990). Due to spherical symmetry of our simulations the magntic field is assumed to be nondirectional in the shell and the synchrotron emissivities and absorption coefficients are angle-averaged.

We have explored the parameter space by varying the accretion rate — corresponding to a variation of the parameter $`ł_{}`$ —, the magnetic field — corresponding to a variation of $`\mu _{30}`$ —, and the amplitude of wave turbulence, $`\delta ^2`$. q is set = 5/3 in all runs presented in this paper.

In Fig. 3, we demonstrate the effect of a varying accretion rate in the case of no turbulence, $`\delta ^2=0`$, and for a fixed magnetic moment of the neutron star, $`\mu _{30}=1`$, corresponding to a strong surface field of $`B_{\mathrm{surf}}10^{12}`$ G. As the accretion rate decreases, the Alfvén radius moves outward, implying a lower magnetic field at the Alfvén radius, and a lower ion temperature and Thomson depth of the active shell. At the same time, however, Compton cooling becomes less efficient due to the reduced soft photon luminosity of the neutron star surface and due to the larger distance of the active region from the surface. This reduction of the soft photon compactness leads to an increasing equilibrium electron temperature in the accretion shell. Consequently, as the accretion rate decreases, the photon spectra change in the following way: For $`ł_{}=1`$, the hard X-ray spectrum smoothly connects to the peak of the thermal blackbody bump from the neutron star surface, and shows a quasi-exponential cutoff at high energies. For lower $`ł_{}`$, the spectrum turns into a typical low-hard state spectrum of X-ray binaries with the thermal bump clearly distinguished from a hard power-law + exponential cutoff at high energies. The hard X-ray power-law becomes harder with decreasing $`ł_{}`$. Fig. 4 shows the same sequence of decreasing $`l_{}`$ for $`\mu _{30}=10^3`$, corresponding to $`B_{\mathrm{surf}}10^9`$ G. At high accretion rates, a strong reduction of the equilibrium electron temperature with respect to the strong-field case results, leading to a softening of the photon spectrum with decreasing $`\mu _{30}`$. This is due to the increasing magnetic field at the Alfvén radius as the neutron star magnetic moment decreases (because the Alfvén radius decreases, overcompensating for the decreasing magnetic moment, see Eq. 8), resulting in a lower electron temperature caused due to increasing cyclotron/synchrotron cooling. For low accretion rates ($`ł_{}`$ a few %) the electron and photon spectra for the strong-field case and the low-field case are virtually indistinguishable from each other.

Figs. 5 and 6 illustrate the effects of a varying accretion rate and turbulence level on the electron and photon spectra. An increasing turbulence level, obviously, leads to a more pronounced nonthermal tail or bump in the electron spectrum. At the same time, as a result of increased Compton cooling on synchrotron photons produced by the nonthermal electrons, the temperature of the quasi-thermal portion of the electron spectrum decreases as the turbulence level is increasing. Consequently, the photon spectrum at soft to medium-energy X-rays becomes softer with increasing $`\delta ^2`$, while at hard X-ray and $`\gamma `$-ray energies an increasingly hard tail develops. As in the case without turbulence, a decreasing accretion rate leads to a higher electron temperature and the transition from a smoothly connected thermal blackbody + thermal Comptonization power-law spectrum to a typical low/hard state X-ray binary spectrum.

Fig. 7 demonstrates the moderate dependence, in particular of the resulting photon spectra, on the magnetic moment of the neutron star for an intermediate value of the accretion rate, $`ł_{}=0.25`$. Nonthermal tails in the electron spectra become more pronounced with increasing neutron star magnetic moment. For the weak-field case with $`\mu _{30}=10^3`$, we find that the high-energy end of the electron spectra are always truncated with respect to a thermal distribution due to strong Compton losses, rather than developing a non-thermal tail.

At very low accretion rates, resulting in very low proton and electron densities in the accretion shell, we find that even at low turbulence level the heating due to stochastic acceleration strongly dominates over Coulomb heating. Consequently, the resulting electron spectrum becomes strongly non-thermal. We need to point out that in those cases, our treatment is no longer self-consistent since in our code the attenuation of Alfvén waves is calculated under the assumption that it is dominated by Landau damping in a thermal background plasma.

The results illustrated in Figs. 4 and 6 are in excellent agreement with the general trend (Barret & Vedrenne (1994); Tavani & Liang (1996)) that hard tails in LMXBs and soft X-ray transients containing weakly magnetized neutron stars are only observed during episodes of low soft-X-ray luminosity. While at large soft-X-ray luminosity, the X-ray emission out to $`20`$ keV is dominated by the thermal component from the neutron star surface and the hard tail is very soft with a low cut-off energy at $`E_\mathrm{c}100`$ keV, in lower-luminosity states the hard X-ray tail becomes very pronounced with photon indices $`\mathrm{\Gamma }2`$ – 3 and extends out to $`E_\mathrm{c}1`$ MeV. Dedicated deep observations by the upcoming INTEGRAL mission should be able to detect this hard X-ray emission from neutron-star-binary soft X-ray transients and atoll sources in both high and low soft-X-ray states and thus provide a critical test of this type of accretion model for weakly magnetized neutron stars.

Interestingly, a similar anti-correlation of the hard X-ray spectral hardness with the soft X-ray luminosity seems to exist for high surface magnetic fields (see Figs. 3 and 5) only in the case of very low accretion rates ($`ł_{}0.01`$) and rather high magnetic turbulence levels ($`\delta ^20.01`$). The hard X-ray spectral indices resulting from our simulations are in excellent agreement with the values of $`\mathrm{\Gamma }3`$ – 4 generally observed in anomalous X-ray pulsars. This may provide additional support for accretion-powered emission models for AXPs. Our simulations predict cutoff energies of $`E_\mathrm{c}100`$ – 500 keV.

## 5 Summary and Conclusions

We have reported on the development of a new, time-dependent code for radiation transport and particle dynamics. The radiation transport, accounting for Compton scattering, bremsstrahlung emission and absorption, cyclotron and synchrotron emission and absorption, and pair processes, is done using a Monte-Carlo method, while the electron dynamics, including radiative cooling, Compton heating/cooling, and stochastic acceleration by resonant interaction with Alfvén/whistler wave turbulence, are calculated using an implicit Fokker-Planck scheme, coupled to the Monte-Carlo radiation transfer code.

In the second part of this paper, we have applied our code to the static situation of a shell at the Alfvén radius of a magnetized neutron star. This situation is representative for accretion onto weakly magnetized neutron stars (bursting atoll sources or soft X-ray transients containing wekaly magnetized neutron stars) as well as to recently suggested models of accretion-powered emission from anomalous X-ray pulsars.

The main results of our parameter study are:

(1) The lepton thermal temperature increases, and the hard X-ray photon spectrum becomes harder as the accretion rate is decreasing. At the same time, the normalization of the hard X-ray power-law, relative to the thermal blackbody from the neutron star surface, becomes smaller.

(2) The nonthermal tails in the electron and photon spectra become more dominant and harder as the turbulence level is increasing. At the same time, the quasi-thermal electron temperature decreases, leading to a softer hard-X-ray spectrum.

(3) For low accretion rates ($``$ a few %), the photon spectra are only very weakly dependent on the magnetic moment of the neutron star. For higher accretion rates, an increasing neutron star magnetic moment leads to a moderate hardening of the hard X-ray spectrum.

Our results are in good agreement with the non-detection of hard X-ray tails during soft X-ray high states of systems containing weakly magnetized neutron stars. However, we predict that a hard X-ray excess beyond $`20`$ keV should exist even in the high/outburst state. This excess should have a cutoff energy of $`E_\mathrm{c}100`$ – 200 keV. In the low/quiescent state, the predicted spectral indices of $`\mathrm{\Gamma }2`$ – 3 are in excellent agreement witht he observed hard X-ray excesses observed from sources believed to contain weakly magnetized neutron stars. We predict that these spectra should extend out to $`E_\mathrm{c}1`$ MeV.

In the strong-field case, representative of a recently suggested model for anomalous X-ray pulsars, we expect a strong correlation between the hard X-ray spectral hardness and the soft X-ray luminosity only in the case of very low accretion rates and high magnetic turbulence level. The predicted hard X-ray spectral indices are generally in very good agreement with the observed values of $`\mathrm{\Gamma }3`$ – 4. Cutoff energies of $`E_\mathrm{c}100`$ – 500 keV are predicted.

Dedicated, deep observations by the upcoming INTEGRAL mission should be able to establish the existence of the hard power-law tails predicted in the models discussed here, and to constrain the high-energy cutoff of these tails. These measurements will be essential ingredients for a more detailed modeling of the physical conditions governing the accretion onto magnetized neutron stars.

The work of MB is supported by NASA through Chandra Postdoctoral Fellowship Award No. 9-10007, issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS 8-39073. This work was partially supported by NASA grant NAG5-4055. We wish to thank the referee for very helpful comments, and D. Marsden for pointing out the relevance of our work to anomalous X-ray pulsars.

## Appendix A Estimate of equilibrium temperature

The code outlined above can be used for equilibrium situations by simply letting it evolve until both the photon and electron distributions have relaxed to a steady state. In order to provide a consistency test of our numerical methods, we compute a quasi-analytical estimate for the expected equilibrium electron temperature. In our equilibrium simulations, we use this estimated equilibrium temperature as initial condition in order to accelerate the convergence.

Assuming that the electron distribution is roughly thermal, we may estimate the heating and cooling rates due to the various processes as follows. Let

$$W=\frac{1}{n_e}\frac{dE}{dV}=\frac{3}{2}kT_e$$
(A1)
be the average energy per electrons. We approximate the cyclotron emissivity by the high-frequency continuum limit given in Mahadevan et al. (1996):

$$j_\nu (\mathrm{\Theta }_e)=\frac{2^{1/6}}{3^{5/6}}\frac{\pi ^{3/2}e^2n_e\nu }{cK_2\left(\frac{1}{\mathrm{\Theta }_e}\right)v^{1/6}}\mathrm{exp}\left[\left(\frac{9v}{2}\right)^{1/3}\right],$$
(A2)
where $`\mathrm{\Theta }_e=kT_e/(m_ec^2)`$, $`v\nu /(\nu _c\mathrm{\Theta }_e^2)`$, and $`\nu _c=eB/(2\pi m_ec)`$.

Denoting by $`\mathrm{\Delta }R`$ the thickness of the emitting shell or slab, the synchrotron-self-absorption frequency is then determined by solving for

$$1=\mathrm{\Delta }R\kappa _{\mathrm{SSA}}=\mathrm{\Delta }R\frac{j_\nu (\mathrm{\Theta }_e)}{B_\nu (\mathrm{\Theta }_e)}.$$
(A3)
The cyclotron/synchrotron cooling rate may then be estimated as

$$\left(\frac{dW}{dt}\right)_{\mathrm{sy}}=\frac{1}{n_eV}\left(A_{\mathrm{surface}}\underset{0}{\overset{\nu _{\mathrm{SSA}}}{}}𝑑\nu B_\nu (\mathrm{\Theta }_e)+V\underset{\nu _{\mathrm{SSA}}}{\overset{\mathrm{}}{}}𝑑\nu j_\nu (\mathrm{\Theta }_e)\right),$$
(A4)
where $`V`$ and $`A_{\mathrm{surface}}`$ are the volume and the surface area of the emitting shell or slab.

Assuming that the photon field inside the emitting volume is dominated by low-energy photons with $`ϵ\mathrm{\Theta }_e`$, where $`ϵ=h\nu /(m_ec^2)`$, the Compton cooling rate may be approximated as

$$\left(\frac{dW}{dt}\right)_{\mathrm{IC}}=4c\sigma _Tu_{\mathrm{ph}}\mathrm{\Theta }_e\frac{K_3\left(\frac{1}{\mathrm{\Theta }_e}\right)}{K_2\left(\frac{1}{\mathrm{\Theta }_e}\right)}.$$
(A5)
The photon energy density $`u_{\mathrm{ph}}`$ is calculated iteratively from the contributions of the various soft photon fields, and repeated Compton scatterings: Assume that the soft photon energy density $`u_\mathrm{s}=u_{\mathrm{sy}}+u_{\mathrm{br}}+u_{\mathrm{ext}}`$ (sycnhrotron + bremsstrahlung + external radiation field) is centered around a soft photon energy $`\overline{ϵ}_s`$, and that $`\tau _\mathrm{T}1`$. The average energy change of a photon with mean photon energy $`\overline{ϵ}_k`$ (Compton scattered $`k`$ times) changes on average by a factor $`\mathrm{\Delta }\overline{ϵ}_k=\overline{ϵ}_k(4\mathrm{\Theta }\overline{ϵ}_k)`$ in the course of the $`k+1`$. scattering. Then, the total internal photon energy density is

$$u_{\mathrm{ph}}=u_\mathrm{s}\left(1+\underset{n=1}{\overset{\mathrm{}}{}}\tau _\mathrm{T}^n\underset{k=0}{\overset{n1}{}}\left[4\mathrm{\Theta }_e\overline{ϵ}_k\right]\right),$$
(A6)
where $`\overline{ϵ}_0=ϵ_s`$ and $`\overline{ϵ}_{k+1}=\overline{ϵ}_k(1+4\mathrm{\Theta }_e\overline{ϵ}_k)`$.

In the cases we are interested in, the bremsstrahlung photon input will be negligible compared to cyclotron/synchrotron and the external soft photon fields. Thus, we assume $`u_\mathrm{s}=u_{\mathrm{ext}}+u_{\mathrm{sy}}`$ with

$$u_{\mathrm{ext}}=\frac{F_{\mathrm{ext}}}{c}$$
(A7)
and

$$u_{\mathrm{sy}}=\frac{1}{c}\underset{0}{\overset{\nu _{\mathrm{SSA}}}{}}𝑑\nu B_\nu (\mathrm{\Theta }_e)+\frac{\mathrm{\Delta }R}{c}\underset{\nu _{\mathrm{SSA}}}{\overset{\mathrm{}}{}}𝑑\nu j_\nu (\mathrm{\Theta }_e).$$
(A8)
The photon energy $`\overline{ϵ}_s`$ is $`\overline{ϵ}_{\mathrm{ext}}`$ if $`u_{\mathrm{ext}}>u_{\mathrm{sy}}`$, or $`\mathrm{max}\{ϵ_{\mathrm{SSA}},ϵ_0\}`$ else, where $`ϵ_0=(h\nu _c/m_ec^2)`$.

The bremsstrahlung cooling rate is

$$\left(\frac{dW}{dt}\right)_{\mathrm{br}}=\{\begin{array}{cc}\frac{128}{3\sqrt{\pi }}\alpha r_e^2m_ec^3n_e\sqrt{\mathrm{\Theta }_e}\hfill & \text{for }\mathrm{\Theta }_e1\hfill \\ 96\alpha r_e^2m_ec^3n_e\mathrm{\Theta }_e\left(\mathrm{ln}[2\mathrm{\Theta }_e]+0.673\right)\hfill & \text{for }\mathrm{\Theta }_e1\hfill \end{array}$$
(A9)
(Haug (1985)) where $`\alpha =1/137`$ is the fine structure constant and $`r_e`$ is the classical electron radius.

The heating rate due to Alfvén wave turbulence is either specified as a fixed parameter, or can be evaluated using the acceleration rates according to Eq. (2) as

$$\left(\frac{dW}{dt}\right)_\mathrm{A}=m_ec^2\underset{1}{\overset{\mathrm{}}{}}𝑑\gamma \dot{\gamma }_Af_e(\gamma ).$$
(A10)
Finally, the Coulomb heating / cooling rate is

$$\left(\frac{dW}{dt}\right)_{\mathrm{Coul}}=\frac{3}{2}\frac{m_e}{m_p}\mathrm{ln}\mathrm{\Lambda }n_\mathrm{p}c\sigma _T(kT_\mathrm{p}kT_\mathrm{e})h(\mathrm{\Theta }_e,\mathrm{\Theta }_\mathrm{p}),$$
(A11)
where

$$h(\mathrm{\Theta }_e,\mathrm{\Theta }_\mathrm{p})=\frac{\sqrt{\mathrm{\Theta }_e}}{(\mathrm{\Theta }_e+\mathrm{\Theta }_\mathrm{p})^{3/2}}\frac{2(\mathrm{\Theta }_e+\mathrm{\Theta }_\mathrm{p})^2+2(\mathrm{\Theta }_e+\mathrm{\Theta }_\mathrm{p})+1}{K_2\left(\frac{1}{\mathrm{\Theta }_e}\right)}\mathrm{exp}\left(\frac{1}{\mathrm{\Theta }_e}\right)$$
(A12)
(Dermer (1986)), $`\mathrm{ln}\mathrm{\Lambda }`$ is the Coulomb logarithm, $`n_p`$ is the number density of ions, and $`\mathrm{\Theta }_p=kT_p/(m_pc^2)`$, .

The total heating/cooling rate consists the sum of all elementary processes,

$$\left(\frac{dW}{dt}\right)_{\mathrm{total}}=\left(\frac{dW}{dt}\right)_{\mathrm{Coul}}+\left(\frac{dW}{dt}\right)_\mathrm{A}+\left(\frac{dW}{dt}\right)_{\mathrm{sy}}+\left(\frac{dW}{dt}\right)_{\mathrm{IC}}+\left(\frac{dW}{dt}\right)_{\mathrm{br}},$$
(A13)
and we are iteratively solving for the equilibrium solution(s) with $`\left(\frac{dW}{dt}\right)_{\mathrm{total}}=0`$.