text:
# Spin–dependent recombination – an electronic readout mechanism for solid state quantum computers

## 1 Introduction

The state of the art of classical computer concepts is rapidly approaching the physical limits as fabrication technology of metal oxide semiconductor logic is minimized to scales where quantum effects determine device properties. While these natural limitations of classical electronics are the dead end for the development of conventional electronics, they open up possibilities for new alternative concepts such as spintronics and quantum computing (QC) . A concept for a silicon based solid state spin–quantum computer as outlined by Kane , combines the advantages of conventional semiconductor technology with regard to the high degree to which this technology has been developed and the fundamental concepts of QC, the massive parallel processing of information by coherent quantum states. Kane’s concept takes advantage of the two nuclear–spin energy eigenstates of 31P-donor nuclei which can be used as well–isolated (long relaxation times) quantum bits (qbits) if they are embedded in a nuclear spin–free crystalline 28Si matrix. Interaction between these nuclear spin qbits can be controlled by electric fields from charged metal gates above and between the donor atoms which can selectively increase the hyperfine interaction between the localized electron donor states as well as the exchange interaction between electron donor states of different 31P-atoms . Before an implementation of the silicon based spin–QC is possible, many technological challenges have to be overcome among which the problem of a single spin readout is particularly difficult. In the original proposal , the readout of nuclear spin states is done by charge measurements of the qbits’ electronic shell which can contain one or two donor electrons from adjacent 31P-atoms. Recently, other proposals for the measurement of a single nuclear spin state have been made utilizing single electron transistors  or spin–transport in combination with spin refrigeration/spin–readout devices .

In this study, a concept is presented, which utilizes spin–dependent charge carrier recombination in silicon for nuclear spin measurements. Electronic transitions between localized, singly occupied

paramagnetic band gap states in silicon have known to be spin–dependent since Lepine  discovered that electron spin resonant (ESR) changes of such defect states can change recombination rates and therefore photoconductivity. The spin selection rules on these recombination transitions arise from the weakness of spin–orbit coupling in silicon which imposes spin–conservation. Thus, before a transition into a singly occupied state can occur, the two electrons participating have to be in a spin state with singlet content. A model of spin–dependent recombination has been developed by Kaplan et al. , who suggested, that spin–dependent recombination processes are always preceeded by a selective pair formation of the recombining electrons and holes. Figure 1 illustrates an example of such a process, which was observed by Kanschat et al.  in hydrogenated microcrystalline silicon ($`\mu `$c-Si:H). Due to the disorder of this material, shallow trap states exist close to the conduction band which can be singly occupied with excess conduction electrons (CE) at low temperatures. Beside these CE centers, additional dangling bond defect states (db) exist which are excellent recombination centers. A formation of a spin–pair as described by Kaplan et al.  takes place, when a CE state in proximity of a db is occupied. Similar to this example, spin–dependent recombination mechanisms exist in other silicon based materials such as dangling–bond recombination in amorphous silicon , various donor–acceptor and donor radiation–defect recombination pathes in crystalline silicon (c-Si) , in silicon devices  and even at c-Si/silicon dioxide interfaces  where spin–dependent recombination occurs between stress induced traps and deep interface states.

With regard to silicon QC, the selectiveness of these spin–dependent recombination mechanisms raises the question, whether they could be utilized for a readout of electronic and hence also nuclear spin states, especially since 31P donors are known to show nuclear spin resonance influence on recombination . Precondition for the information readout of coherent spin states with recombination processes is the ability of a given transition to reflect the coherence of the spin states involved, which means the spin–pair states that determine whether a recombination transition takes place or not must not fall into one of the four eigenstates, before a second transition, the actual electronic transfer takes place. This however has never been proven in the past since the experiment used for the detection of all the processes mentioned above (often refered to as electrically detected magnetic resonance, EDMR) is a pure incoherent steady–state magnetic field–sweep experiment, where a constant current of light or electrically injected excess charge carriers is measured while a constant microwave radiation is imposed on the sample. When the magnetic field reaches the resonance of one of the involved centers, a change of the steady state recombination and thus a current change takes place. While EDMR shows that certain magnetic centers have an influence on recombination, it fails completely to reveal any information about coherence times of the systems involved which are not only the spin–relaxation times of the respective spin centers but also electronic transition times such as recombination or dissociation of the spin pair.

In the following an experiment is presented that shows how coherent spin precession of electrons in localized band gap states can govern recombination rates. Based on this observation we propose an architecture for a recombination based readout of 31P-qbits.

## 2 Coherent spin motion effects on recombination

Time–domain measurements of spin–dependent recombination (TSR) were carried out on the CE-db mechanism in $`\mu `$c-Si:H, that is mentioned above. TSR is the electrical detection of pulsed ESR and therefore the time resolved equivalent of

EDMR . It combines the advantages of EDMR with regard to the selective detection of distinct recombination processes and the advantages of pulsed ESR with regard to its ability to detect coherent spin-motion effects.

The idea of the experiment carried out is to manipulate the steady state of the CE-db pair–ensemble with pulsed ESR on a nanosecond time scale and to observe the transient photocurrent response that is determined by the recombination rate through these centers. Excess charge carriers are generated by a continuous light source. Since singlet states have only short lifetimes, a high density of triplet states is present in the steady state when resonant microwave radiation is absent .

The orientation of the magnetic moments of $`|T_{}`$–pairs is illustrated in fig. 2 in the rotating frame Bloch–sphere representation  of the two respective spins contained in the charge carrier pairs. If a microwave that is in resonance with the dangling bond and polarized along the $`x`$–axis is switched on, Rabi–oscillations of the db spins take place. The influence of the microwave pulses on the CE centers is small due to their different resonance frequencies. The Rabi–oscillation of the db spin would lead to an oscillating recombination rate which could be detected as an oscillation of the photocurrent, if the inhomogeneities present in $`\mu `$c-Si:H  did not cause a rapid dephasing of the db ensemble within less than one oscillation period (illustrated in fig. 2). This dephasing is still a coherent process, however, its effect on the recombination rate (strong attenuation of the Rabi induced oscillation) makes it indistinguishable from incoherence. Thus, a rephasing of the ensemble has to be induced, so that incoherence and dephasing becomes distinguishable. As illustrated in fig. 2, a reversal of the direction of spin precession caused by a reversed $`B_1`$-microwave field causes such a rephasing. The $`B_1`$–field reversal, done by a microwave phase change that is introduced at an arbitrary time $`\tau _{180}`$, causes a maximum recovery of the $`|T_{}`$–density at time $`2\tau _{180}`$ and hence, a minimum of the recombination takes place. This minimum, which we call a recombination echo  will only be detectable, if incoherent processes are not faster than the time scale of the pulse sequence. The detection is a proof of the coherence of the charge carrier pair’s spin motion.

Experimental details with regard to the $`\mu `$c-Si:H sample and the TSR setup are identical to those of ref.  execept that the laser intensity was 300mW, the sample temperature $`T`$ = 30K and the intensity of the coherent microwave radiation 500W.

The experimental results are displayed in figs. 4 and 3. Note that the detection of fast photocurrent changes in the pA-range with constant offsets in $`\mu `$A-range is extremely difficult due to noise limitations. The inital decrease of the real–time photocurrent displayed in fig.4 is therefore due to the rise time of the detection setup. As explained elsewhere , the amplitude of the exponential photocurrent relaxation back to the steady state is determined by the coherent spin state prepared during the ns-pulse sequence. Hence, a measurement of the real time signal at a certain time after the rise time versus the pulse length reveals the spin motion during the pulse sequence on a ns-scale. This pulse length dependence, measured at the photocurrent minimum ($`t=19.5\mu `$s) after the pulses, is displayed in fig.3. It clearly displays the photocurrent increase after a sequence of two 100ns long microwave pulses with equal $`B_1`$-field strength but opposite phase orientations. The width of the echo displayed in fig. 3 turns out to be antiproportional to the $`B_1`$-field strength which proves the involvement of Rabi–oscillation in the recombination echo effect (not shown here). The two real time transients displayed in fig.4 are recorded after resonant excitation of the db center.

In both cases, the total length and the power of the pulses were identical and only a phase change of 180o was introduced in one case. The smaller amplitude of one transient is solely due to this phase change after 100ns which led to the rephasing of triplet states. An identification of singlet and triplet states with the two binary digits (labeled “1” and “0” for instance) used for information processing allows a direct readout with current measurements.

## 3 Spin–dependent recombination and nuclear spin measurements

Based on the observation above, a setup for a recombination readout of Kanes’ silicon based quantum computer is proposed as illustrated in fig. 5. The device proposed consists of a point defect that introduces a deep level paramagnetic state in proximity to the 31P-donor which is to be read out, a readout–gate (R-gate), and an additional A-gate above the defect. The latter is only necessary if the defect is induced by a deep level donor with nuclear spin ($`I0`$). When the nature and the implementation of the paramagnetic deep level state will be discussed in the following section, it will be refered to as db in this section in consistency with the dangling bond used for the experimental demonstration discussed above.

The idea behind this setup is to take advantage of Pauli’s principle, similar to the charge–readout proposed by Kane . As long as the R-gate is charged negatively, the wave function overlap between the donor electron of the 31P and the db is small which minimizes a transition probability between the two states and keeps the spin–exchange negligibly small as well. When the R-gate is charged positively, such that the wave function overlap between the 31P and the db is sufficiently large, exchange interaction increases. If this increase is introduced slowly, an adiabatic change of the spin-pairs’ energy eigenstates can take place from an uncoupled product base into a set of singlet and triplet states. At low temperatures ($`T100`$mK), the uncoupled pair is polarized in a $`|T`$-state as long as the coupling is absent. When the exchange is increased slowly, this $`|T`$-state remains either unchanged or shifts into a $`|S`$-state, depending on the orientation of the 31P-nucleus (illustrated in fig. 6).

Note that in presence of a second nuclear spin at the db-site, the A-gate above this site would have to be charge positively in order to minimize hyperfine coupling. After the R-gate has been “opened” and the electronic pair is in singlet or triplet state, an electronic transition can take place which charges the db; this transition however is possible only when the electronic pair has singlet content in advance and thus, the charging of the db state depends on the nuclear state of the 31P qbit.

The timing and a scetch of the transitions during a readout sequence are displayed in fig. 7. As long as quantum operations on the qbit take place, no light is imposed on the sample and the R-gate is charge negatively.

When the readout sequence begins at a time $`t_0`$, the bias of the R-gate is slowly inverted towards a positive voltage. Once exchange coupling is established after a time $`\tau _{\mathrm{s}lope}`$, it remains unchanged until the actual recombination transition has taken place. The time necessary for this transition is of the order of the electronic pair states’s lifetime $`\tau _{\mathrm{l}ife}`$, after which the exchange interaction is switched off again. The nuclear spin state of the 31P donor is then coded in the donor and db charge state. A short ($`\tau _{\mathrm{f}lash}`$, ns-Range) and weak (nW-range) laser pulse imposed on the sample will then increase conductivity by the generation of a few pairs of excess electrons and holes. If the donor and the db state are not charged (no transition), a slow decrease of the photoconductivity will follow, which is determined by slow band–band recombination in the ultra pure 28Si bulk. If the two states are charged, a fast decay will take place since charge carriers will be trapped in the charged states. Thus, the level of the photoconductivity a time $`\tau _{\mathrm{d}ecay}`$ after the end of the laser pulse will reveal the result of the readout process. Note, that the readout process itself automatically neutralizes the two states such that a new series of operations can take place after its completion.

## 4 Challenges of the implementation

The experiments presented above are the motivation for the concept proposed in this study and many other questions will have to be answered experimentally before an actual proof of this concept can be given. The implementation of the singly occupied deep level state will mostly depend on whether it is possible to control its location and whether its charge carrier capture cross sections and the given transition times of the system are fast enough. The dangling bond present in amorphous and microcrystalline silicon would be an ideal system with regard to the latter properties. However, the high disorder in the microcrystalline morphology of silicon makes it a bad choice for QC. Since no process has been established which allows the creation of a single db with $`\AA `$-accuracy at an arbitrary site, a different way of deep level implementation must be chosen.

Various impurities and defects provide singly occupied, paramagnetic deep level donor states in c-Si, some of which are gold, potassium, strontium or chromium . An additional possibility could be interface recombination. Recombination at c-Si/SiO2–interfaces has shown to have spin-dependent paths . Hence, the very same processes which limited the original pure c-Si QC-concept and led to the proposal of Si/SiGe-heterostructures  could actually be benificial for a recombination readout mechanism.

Beside the questions for a proper implementation of the deep level center, several other issues remain to be investigated: Due to the necessity of error correction, one of the preconditions for QC is that readout does not destroy the state of the qbit. Therefore, only P-deep level–transitions are possible which leave the 31P–relaxation unchanged. Important for the feasibility of the recombination readout will also be the question of the overall readout time which should not exceed a lower microsecond range. When the P-deep level transition takes place, the slow relaxation of non–equilibrium polarization at low temperatures could pose a problem due to charged P-deep level–pairs which are neutralized by excess charge carriers produced in excited spin states. The time needed for the relaxation of these states could be very long. However, this problem may be solved by injection of excess carriers by spin–polarized electronic injection instead of light injection. Finally, the question for the one–qbit sensitivity has to be asked: Time–domain measurements of spin–dependent recombination have shown to be sensitive enough to detect the influence of as little as 50 charge carrier pairs, even though the experiment was carried out in the presence of a strong constant photocurrent offset. For the readout as proposed, a measurement without offset could be made and the detection would be even more sensitive. A 2eV laser pulse of 1ns length and 3.2nW intensity produces 10 electron–hole pairs if the internal quantum efficiency is assumed to be 100%. The detection of these 10 charge carrier pairs would require a pA–current measurement on a microsecond scale which does not pose a problem – however, whether the 10% difference between the two readout results caused by the one recombination center are detectable can only be proven when an implementation of the proposed setup is available.

## 5 Summary

The demonstration that coherent spin–states of electrons and holes trapped in CE and db states can govern their recombination rate has motivated the idea of a readout mechanism for 31P qbits in c-Si. A concept, based on a sequence of exchange gating, light pulse and photocurrent measurements was proposed, various implementations suggested and possible limitations and drawbacks discussed.


format: text

text:
# Studies of the Relativistic Binary Pulsar PSR B1534+12: I. Timing Analysis

## 1 Introduction

The discovery (Hulse & Taylor, 1975) of the first pulsar in a binary system, PSR B1913+16, introduced the possibility of powerful new experimental tests of gravity in the strong field and radiative regimes. Over a quarter century, measurements of this pulsar were shown to be in excellent agreement with the predictions of general relativity (Taylor & Weisberg, 1989; Damour & Taylor, 1991; Taylor, 1994). A dozen years ago the discovery (Wolszczan, 1991) of PSR B1534+12—a bright, relatively nearby pulsar, also with a neutron-star companion—promised the independent opportunity to confirm the tests of gravity done earlier with PSR B1913+16. Furthermore, B1534+12’s favorable, nearly edge-on orbital geometry allowed new experimental tests that were complementary to those done with PSR B1913+16 (Taylor et al., 1992).

In earlier work (Stairs et al., 1998, hereafter Paper I) we analyzed data taken over a six year period with the Arecibo 305 m telescope, the 43-m telescope at Green Bank, West Virginia, and the 76-m Lovell Telescope at Jodrell Bank Observatory, U.K. We found that those data were in agreement with the predictions of general relativity at the limit of the measurement uncertainties: about 1%. We showed that the orbital period was decreasing, as expected if the system is losing energy in the form of gravitational radiation. Because the observed orbital decay rate is contaminated by kinematic effects which depend on the unknown pulsar distance (Damour & Taylor, 1991; Camilo, Thorsett & Kulkarni, 1994), the radiation predictions of general relativity cannot be tested with this pulsar to better than $`30\%`$ without an independent distance measurement. However, with the assumption that general relativity is the correct theory of gravity, we were able to invert the kinematic model to measure the pulsar distance to 20%.

Over the past four years, we have continued our studies of PSR B1534+12, using the recently upgraded Arecibo telescope together with a new data acquisition system (Stairs et al., 2000a) that fully removes the dispersive effects of the interstellar medium from the signal. In addition to extending the length of the data set by more than half, these new data are of substantially higher quality than those available for Paper I. In this paper, we report on timing studies based on the complete 11.5-year Arecibo data set. We present the improved tests of general relativity that are now possible, give a pulsar distance accurate to 5%, and report on dispersion measure variations. In subsequent papers, we will describe the analysis of secular variations in the pulsar radiation pattern that have been interpreted as evidence for geodetic precession of the pulsar, and we will discuss single-pulse studies.

## 2 Observations

In this paper we discuss results of observations made with two different observing systems at the Arecibo telescope. Many of the details have already been discussed in Paper I, and are only summarized here and in Table 1.

Observations were made with the 305-m Arecibo telescope both before its upgrade (1990–1994) and after its upgrade (1997–2002). For our purposes, the most important post-upgrade change was a significant reduction in the system temperature at 1400 MHz, leading to an improvement in the signal-to-noise ratio achieved for a given integration time and observing bandwidth. During the period that the Arecibo telescope was unavailable, observations were made with the National Radio Astronomy Observatory (NRAO) 43-m telescope at Green Bank, and with the 76-m Lovell telescope at Jodrell Bank Observatory, U.K. The latter two datasets were used in Paper I, but, with the abundance of high-quality post-upgrade Arecibo data, they contribute nothing to the present solution and are therefore disregarded in this work.

Two observing systems were used. The Princeton Mark III system (Stinebring et al., 1992) is an “incoherent” system, in which the incoming signal is decomposed into spectral channels with a filterbank. The effective time resolution is dominated by dispersion smearing of the signal during its propagation through the ionized interstellar medium, with some additional instrumental smoothing after detection of the signal. The signals in each frequency channel were folded at the instantaneous topocentric pulsar period for three to five minutes, then the channels were shifted to account for interstellar dispersion and summed to produce a single total-intensity (summed-polarization) pulse profile for each integration.

For the new post-upgrade data that we introduce here, we observed with a “coherent” system: the Princeton Mark IV instrument (Stairs et al., 2000a). For each sense of circular polarization, a 5 MHz bandpass was mixed to baseband using local oscillators in phase quadrature. The four resulting signals were low-pass filtered at 2.35 MHz, sampled at 5 MHz with 4-bit quantization, and written to a large disk array. Some of the 1400 MHz observations were taken with 10 MHz of bandwidth and 2-bit quantization. In the offline processing, the undetected signal voltages were dedispersed using the phase-coherent technique described by Hankins & Rickett (1975). After amplitude calibrations were applied, self- and cross-products of the right- and left-handed complex voltages yielded the Stokes parameters of the incoming signal. These products were folded at the topocentric pulsar period using 1024 phase bins, and pulse TOAs were determined from the resulting total-intensity profiles. A summary of the more important parameters and statistics of all combinations of observing systems and telescopes is presented in Table 1. Many of the post-upgrade data come from observations on one or two days at quasi-regular intervals; in addition we conducted 12-day intensive “campaign-style” observations in each of the summers of 1998, 1999, 2000 and 2001, with full orbital coverage at each epoch. This mixture of data permits a high-quality timing solution and allows us to monitor, at each campaign epoch, the pulse shape changes due to geodetic precession, which we will discuss in a subsequent paper.

As in Paper I, we used the same TOA-fitting procedure for all data sets. Each observed profile was fitted to a standard template, using a least-squares method in the Fourier transform domain (Taylor, 1992) to measure its time offset. The offset was added to the time of the first sample of a period near the middle of the integration, thereby yielding an effective pulse arrival time. A different standard template was used for each observing system and frequency; they were made by averaging the available profiles over several hours or more. Uncertainties in the TOAs were estimated from the least squares procedure, and a minimum error, dependent on the observing system and frequency, was required for each TOA. In addition, a cut was established such that all Mark IV TOAs with least-squares uncertainties greater than 10 $`\mu `$s (roughly 20% of the data) were discarded; tests show that they do not significantly affect the best-fit solution or the parameter uncertainties, and discarding these data points allows us to better assess the data quality (§3.3). The observatory’s local time standard was corrected retroactively to the UTC timescale, using data from the Global Positioning System (GPS) satellites.

## 3 Data Analysis

### 3.1 The Timing Model

A pulse received on Earth at topocentric time $`t`$ is emitted at a time in the comoving pulsar frame given by

$$T=tt_0+\mathrm{\Delta }_CD/f^2+\mathrm{\Delta }_R+\mathrm{\Delta }_E\mathrm{\Delta }_S\mathrm{\Delta }_R\mathrm{\Delta }_E\mathrm{\Delta }_S.$$
(1)
Here $`t_0`$ is a reference epoch and $`\mathrm{\Delta }_C`$ is the offset between the observatory master clock and the reference standard of terrestrial time. The dispersive delay is $`D/f^2`$, where $`D=\mathrm{DM}/2.41\times 10^4`$, with the dispersion measure DM in cm-3pc, the radio frequency $`f`$ in MHz, and the delay in seconds. The $`\mathrm{\Delta }_R`$, $`\mathrm{\Delta }_E`$, and $`\mathrm{\Delta }_S`$ terms are propagation delays and relativistic time adjustments for effects within the solar system, and $`\mathrm{\Delta }_R`$, $`\mathrm{\Delta }_E`$ and $`\mathrm{\Delta }_S`$ are similar terms accounting for phenomena within the pulsar’s orbit.

In Paper I, we set out the equations that define the three orbital terms $`\mathrm{\Delta }_R`$, $`\mathrm{\Delta }_E`$, and $`\mathrm{\Delta }_S`$ in terms of ten parameters. Five are the standard Keplerian parameters: the period $`P_b`$, projected semi-major axis $`xa_1\mathrm{sin}i/c`$, eccentricity $`e`$, longitude of periastron $`\omega `$, and time of periastron $`T_0`$. The other five are Post Keplerian (PK) parameters: the rate of advance of the periastron $`\dot{\omega }`$, orbital period derivative $`\dot{P}_b`$, time dilation and gravitational redshift $`\gamma `$, and range $`r`$ and shape $`s`$ of the Shapiro time delay. These quantities, in conjunction with the time-variable dispersion measure, a simple time polynomial to model the spin of the pulsar, and astrometric parameters to model the propagation of the signal across the solar system, constitute the free parameters to be fit in the theory-independent timing model.

Within the framework of a particular theory of gravity, the five PK parameters can be written as functions of the pulsar and companion star masses, $`m_1`$ and $`m_2`$, and the well-determined Keplerian parameters. In general relativity the equations are as follows (see Damour & Deruelle, 1986; Taylor & Weisberg, 1989; Damour & Taylor, 1992):

$`\dot{\omega }`$ $`=`$ $`3\left({\displaystyle \frac{P_b}{2\pi }}\right)^{5/3}(T_{}M)^{2/3}(1e^2)^1,`$ (2)
$`\gamma `$ $`=`$ $`e\left({\displaystyle \frac{P_b}{2\pi }}\right)^{1/3}T_{}^{2/3}M^{4/3}m_2(m_1+2m_2),`$ (3)
$`\dot{P}_b`$ $`=`$ $`{\displaystyle \frac{192\pi }{5}}\left({\displaystyle \frac{P_b}{2\pi }}\right)^{5/3}\left(1+{\displaystyle \frac{73}{24}}e^2+{\displaystyle \frac{37}{96}}e^4\right)(1e^2)^{7/2}T_{}^{5/3}m_1m_2M^{1/3},`$ (4)
$`r`$ $`=`$ $`T_{}m_2,`$ (5)
$`s`$ $`=`$ $`x\left({\displaystyle \frac{P_b}{2\pi }}\right)^{2/3}T_{}^{1/3}M^{2/3}m_2^1.`$ (6)
Here the masses $`m_1`$, $`m_2`$, and $`Mm_1+m_2`$ are expressed in solar units, $`T_{}GM_{}/c^3=4.925490947\mu `$s, $`G`$ the Newtonian constant of gravity, and $`c`$ the speed of light.

### 3.2 Arrival Time Analysis

We used the standard tempo analysis software (Taylor & Weisberg, 1989, see also http://pulsar.princeton.edu/tempo) together with the JPL DE200 solar-system ephemeris (Standish, 1990) to perform a least-squares fit of the measured pulse arrival times to the timing model. Results for the astrometric, spin, and dispersion parameters of PSR B1534+12 are presented in Table 2. The DM was fit in two blocks, pre-and post-upgrade, with different base values and time derivatives in each block. Because of systematic errors (see §3.3 below), the Mark III 430 MHz data were only used to fit the pre-upgrade DM and DM derivative, and were not incorporated into the rest of the timing solution. This meant that an iterative approach was required to find the best-fit DMs: first the spin and orbital solution was held constant while the DM parameters were fit using the whole dataset, then the Mark III 430 MHz data were dropped as the best overall timing solution was obtained from the other data, holding the pre-upgrade DM parameters fixed. This new timing solution was then used to refine both sets of DM parameters, until no significant changes resulted. The DM and DM-derivative uncertainties reported in Table 2 are the errors derived while holding the rest of the timing solution fixed. Fitting for post-upgrade DM and DM-derivative simultaneously with the rest of the timing solution yielded similar values, but somewhat larger uncertainties due to covariances with other parameters. Figure 1 shows the time variation of the best-fit DM, with the DM and first derivative from Paper I shown for comparison; the pre-upgrade values are consistent within the uncertainties. We note that the offset in DM base values between the pre- and post-upgrade fits is due mostly to the use of different standard profiles with necessarily different frequency alignments; thus this offset should not be interpreted as physically meaningful. The change in DM derivative, however, is real, and is discussed in §4.2 below. An arbitrary time offset was also allowed between the pre- and post-upgrade data sets to accommodate the different standard profiles and changes in the instrumental signal path length.

Two models of the pulsar orbit were used to fit the data. The theory-independent “DD” model (Damour & Deruelle, 1986) treats all five PK parameters defined in §3.1 as free parameters in the fit. The alternate “DDGR” model (Taylor, 1987; Taylor & Weisberg, 1989) assumes that general relativity is correct and uses equations 2 through 6 to tie the PK parameters to $`Mm_1+m_2`$ and $`m_2`$; consequently it requires only two post-Keplerian free parameters.

Table 3 presents our best-fit orbital parameters. Uncertainties given in the table are approximately twice the formal “$`1\sigma `$” errors from the fit; we believe them to be conservative estimates of the true 68%-confidence uncertainties, including both random and systematic effects. For each model, we list the Keplerian parameters and the relevant post-Keplerian parameters. We also include, in italic numbers, the computed PK parameters derived from the measured masses in the DDGR fit; these computed parameters are in excellent agreement with the DD fit values. As in Paper I, the DDGR solution includes the “excess $`\dot{P}_b`$” parameter, which accounts for an otherwise unmodeled acceleration resulting from galactic kinematics. We note that the measured orbital period derivative, $`0.137\pm 0.003\times 10^{12}`$, is very close to the sum of the DDGR model’s intrinsic value, $`0.1924\times 10^{12}`$, and the fitted excess, $`+0.055\pm 0.003\times 10^{12}`$.

Figure 2 shows the post-fit residuals for the pre-upgrade 1400 MHz data, and the post-upgrade 1400 and 430 MHz data, plotted as functions of date. Interstellar scintillation causes the TOA uncertainties to vary somewhat, but in the interest of clarity we have not drawn error bars on these residual plots. Figure 3 illustrates the averaged post-fit residuals for the same datasets, plotted as functions of orbital phase.

### 3.3 Comparison of Coherent and Incoherent Observing Systems

As discussed in Paper I, solutions for the fitted parameters of PSR B1534+12, as determined from the pre-upgrade 430 MHz data, are biased by small systematic errors in the TOAs, likely caused by imperfect post-detection dispersion removal combined with variable interstellar scintillation patterns. These measurements are unreliable at the several microsecond level required for high-precision timing of this pulsar.

One of our major goals in building the Princeton Mark IV instrument (Stairs et al., 2000a) was to improve the timing accuracy at low frequencies for fast systems such as PSR B1534+12, by minimizing any dispersion-related systematic errors in TOAs. We have determined that the new instrument does indeed appear to have circumvented the residual dispersion problems that plague the Mark III data at low frequencies. In Figure 4 we follow Taylor & Weisberg (1989) and Paper I in examining the statistical properties of the post-fit residuals (relative to the solution of Tables 2 and 3) to see if they “integrate down” as $`n^{1/2}`$ when $`n`$ consecutive values are averaged. The pre-upgrade 1400 MHz data and both post-upgrade data sets are fairly consistent with the expected slope, while the pre-upgrade 430 MHz data deviate significantly. This confirms that we have been justified in leaving the pre-upgrade 430 MHz data out of our preferred timing solution, that these Mark III 430 MHz data indeed contain systematic errors related to the incoherent dispersion removal, and that the new Mark IV instrument greatly reduces these systematics, as we had hoped. The minimum TOA uncertainty requirement mentioned above was introduced to compensate for the small remaining deviations from the $`n^{1/2}`$ law visible in the three good datasets.

## 4 Discussion

### 4.1 Observed Change in Orbital Period: Distance to PSR B1534+12

The observed $`\dot{P}_b`$ is measured in the reference frame of the solar system barycenter, and must be corrected to the center-of-mass frame of the binary pulsar system before it can be compared to the value predicted by GR. The largest correction is kinematic, involving the relative acceleration of the two reference frames. It can be broken down into the vertical acceleration in the Galactic potential, acceleration in the plane of the Galaxy, and an apparent centripetal acceleration due to the transverse velocity of the pulsar binary (Damour & Taylor, 1991):

$$\left(\frac{\dot{P_b}}{P_b}\right)^{\mathrm{gal}}=\frac{a_z\mathrm{sin}b}{c}\frac{v_0^2}{cR_0}\left[\mathrm{cos}l+\frac{\beta }{\mathrm{sin}^2l+\beta ^2}\right]+\mu ^2\frac{d}{c}.$$
(7)
Here $`a_z`$ is the vertical component of galactic acceleration, $`l`$ and $`b`$ the pulsar’s galactic coordinates, $`R_0`$ and $`v_0`$ the Sun’s galactocentric distance and galactic circular velocity, $`\mu `$ and $`d`$ the pulsar’s proper motion and distance, and $`\beta =d/R_0\mathrm{cos}l`$. As we do not have a precise distance to the pulsar through either a timing or an interferometric parallax measurement, the best independent distance estimate still comes from the Taylor & Cordes (1993) model of the free electron content of the Galaxy; this model puts the pulsar at a distance of $`0.7\pm 0.2`$ kpc. At this distance the Kuijken & Gilmore (1989) model of the Galactic potential yields an estimate of $`a_z/c=(1.60\pm 0.13)\times 10^{19}\text{s}^1`$. We assume $`v_0=222\pm 20`$km s-1 and $`R_0=7.7\pm 0.7`$ kpc, as in Damour & Taylor (1991). Then, summing the terms in equation (7) and multiplying by $`P_b`$, we find the total kinematic bias to be

$$\left(\dot{P}_b\right)^{\mathrm{gal}}=(0.037\pm 0.011)\times 10^{12}.$$
(8)
The uncertainty in this correction is dominated by the uncertainty in distance, which is only roughly estimated by the Taylor and Cordes model. The slight decrease in the uncertainty from that given in Paper I results from an increase by a factor of ten in the precision of the proper motion measurement.

Our measurement of the intrinsic rate of orbital period decay is therefore

$$\left(\dot{P}_b\right)^{\mathrm{obs}}\left(\dot{P}_b\right)^{\mathrm{gal}}=(0.174\pm 0.011)\times 10^{12}.$$
(9)
The uncertainty is completely dominated by the uncertainty on the kinematic correction, which is nearly a factor of 4 larger than the measurement uncertainty on $`\dot{P}_b^{\mathrm{obs}}`$. In GR, the orbital period decay due to gravitational radiation damping, $`(\dot{P_b})^{\mathrm{GR}}`$, can be predicted from the masses $`m_1`$ and $`m_2`$ (eq. 4), which in turn can be deduced from the high precision measurements of $`\dot{\omega }`$ and $`\gamma `$. As listed in Table 3, the expected value is

$$\left(\dot{P}_b\right)^{\mathrm{GR}}=0.192\times 10^{12}.$$
(10)
As in Paper I, our measured value differs from this prediction, now by some 1.7 standard deviations, and again, assuming that GR is the correct theory of gravity, we can derive the true distance to the pulsar from the “excess $`\dot{P}_b`$” parameter and equation 7 above (Bell & Bailes, 1996). Our improved distance is $`d=1.02\pm 0.05`$ kpc (68% confidence limit). The uncertainty is still dominated by the measurement uncertainty of $`(\dot{P}_b)^{\mathrm{obs}}`$, rather than uncertainties in the galactic rotation parameters or the acceleration $`a_z`$, though the Galactic model uncertainties will ultimately limit the distance measurement to a few percent accuracy. The timing parallax for this system is still not significantly measured, but is constrained to be less than 1.5 mas, in good agreement with the GR-derived result.

Our new distance result is consistent with the $`1.1\pm 0.2`$ kpc determined in Paper I, but we have improved the uncertainty by a factor of 4. Our previous distance led to a downward revision of the estimated double-neutron-star inspiral rate visible to gravitational-wave observatories such as LIGO. The new measurement leads to a small (15%) increase in the number density of similar systems in the local universe, with the uncertainty now dominated by the uncertain scale height of such binaries (e.g., Kalogera et al., 2001).

### 4.2 Dispersion Measurements and the Local Interstellar Medium

As noted in §3.2, we find a significantly different DM and DM derivative in the post-upgrade era than before the upgrade (Figure 1). We argue that while the change in reference DM value is largely due to the use of new standard profiles with a slightly different frequency alignment, the difference in the rate of change (more than a factor of 4 slower since the upgrade) is physical. It is true that the profile at both frequencies is undergoing a secular evolution due to geodetic precession of the pulsar’s spin axis (Arzoumanian, 1995, Stairs et al., in prep.) and thus one might suspect a contribution to the DM derivative from different temporal evolution at the two observing frequencies. We explore this possibility by calculating the DM difference between our 1999 May and 2001 June observing campaigns, for the best-fit post-upgrade DM derivative versus the best-fit pre-upgrade DM derivative, finding that if the pre-upgrade DM trend had continued, the 2001 June DM would have been smaller by about 0.0005 cm-3 pc. This large a difference would imply a difference in offset between the 430 MHz and 1400 MHz profiles of 10 $`\mu `$s, or 0.27 bins in a 1024-bin profile. By comparing the actual cumulative profiles at these epochs using the same cross-correlation algorithm used to calculate TOAs, we estimate that the 430–1400 MHz offset has shifted by no more than about 0.05 bins. As the TOAs are determined by the strong main peak of the pulse and the most noticeable evolution is taking place in the low-level emission near the base of the pulse (Stairs et al., 2000b) this is perhaps not surprising. We therefore conclude that the observed changes in DM slope do represent changes in the structure of the interstellar medium between the Earth and the pulsar.

With our current data set, it is not possible to construct a satisfactory structure function with which to characterize the length scales of turbulence in the interstellar medium (e.g., Phillips & Wolszczan, 1991). Because of strong refractive scintillation at 1400 MHz, the high-frequency data set remains quite sparse (Figure 2). Even though we may use a wider-bandwidth instrument such as the recently commissioned 100-MHz “WAPP” spectrometer for future 1400 MHz observations, the broadband nature of the refractive interstellar scintillation may make it difficult to accurately sample the dispersion measure variations with resolution finer than roughly 6 months.

### 4.3 The Neutron Star Masses: Testing Evolutionary Theory

The best mass estimates for the two neutron stars in the PSR B1534+12 system come from the DDGR timing model. The measured values are $`m_1=1.3332\pm 0.0010M_{}`$ and $`m_2=1.3452\pm 0.0010M_{}`$. The masses derived from the DD model and equations 5 and 6 are in good agreement with these values, though the measurement uncertainties are much larger (Figure 5).

Although the masses of the two neutron stars are very similar, it is now clear that the pulsar is significantly less massive than its companion. This is contrary to initial expectations from binary evolution. The pulsar is rapidly spinning, and has the low magnetic field characteristic of a “recycled pulsar” that has been spun-up by mass transfer from a companion star. Assuming a monotonic relation between progenitor mass and neutron star mass, we therefore might predict that the pulsar would be more massive at birth, and that the mass difference would only increase during the later evolution of the system.

This perhaps naive prediction has already been challenged by the binary pulsar system PSR B2303+46. Recently, the companion in this close, eccentric binary was shown to be a white dwarf rather than a second neutron star (van Kerkwijk & Kulkarni, 1999). The eccentric orbit implies that the white dwarf formed first, from the star that was originally more massive. Presumably, mass transfer to the initially less massive star pushed it above the minimum mass needed to form a supernova and a neutron star.

A similar process may have been important in the formation of the PSR B1534+12 system. Although the pulsar was formed from the initially more massive star, mass transfer to its companion star probably resulted in a mass inversion before either neutron star was formed.

To move beyond the simplest (or most naive) predictions about the relative masses of the pulsar and companion in the B1534+12 system requires significant improvements in our theoretical understanding of several areas of stellar and binary evolution. First, we must understand the relationship between progenitor and neutron star masses, which may not be purely monotonic. Then we must understand how the mass lost by the pulsar’s progenitor or gained by the companion’s progenitor affected the evolution of the stellar cores and the amount of fallback during the supernova events. Finally, we must understand in some detail the mass transfer that spun up the pulsar. Although these are all challenging problems, there is some hope for future guidance from the growing number of precisely measured neutron star masses.

### 4.4 Test of Relativity

PSR B1534+12 permits the second test of general relativity based on the $`\dot{\omega }`$, $`\gamma `$, and $`\dot{P}_b`$ parameters of a binary pulsar system, and adds significant measurements of the Shapiro-delay parameters $`r`$ and $`s`$. The left-hand sides of equations (26) represent quantities measured in a theory-independent fashion and listed in the DD column of Table 3. If GR is consistent with these measurements and there are no significant unmodeled effects, the five curves corresponding to equations (26) should intersect at a single point in the $`m_1`$-$`m_2`$ plane. These curves are presented in Figure 6, in which a pair of lines delimit the 68% confidence limit for each PK parameter (a single line is drawn for $`\dot{\omega }`$, whose uncertainty is too small to show). It is clear that the $`\dot{\omega }`$, $`\gamma `$, $`r`$, and $`s`$ curves intersect, though $`r`$ is still poorly measured. The curve obtained from the observed DD value of $`\dot{P}_b`$ can be made to intersect the others, as discussed above, by setting the pulsar distance to $`1.02\pm 0.05`$ kpc rather than the $`0.7\pm 0.2`$ kpc estimated from the dispersion measure. A filled circle at $`m_1=1.3332M_{},m_2=1.3452M_{}`$ marks the DDGR solution of Table 3, and its location on the $`\dot{\omega }`$ line agrees to within 0.05% with the measured DD values of $`\gamma `$ and $`s`$. This provides a highly precise test of the validity of general relativity using only non-radiative timing parameters, an important complement to the mixed $`\dot{\omega }`$$`\gamma `$$`\dot{P}_b`$ test provided by PSR B1913+16.

The Arecibo Observatory, a facility of the National Astronomy and Ionosphere Center, is operated by Cornell University under a cooperative agreement with the National Science Foundation. We thank Zaven Arzoumanian, Fernando Camilo, Andrew Lyne and David Nice for their earlier contributions to this long-term project, Kiriaki Xilouris, Duncan Lorimer, David Nice, Eric Splaver, Andrea Lommen, Paulo Freire, and Ian Hoffman for assistance with observations, and an anonymous referee for comments on the manuscript.. I. H. S. is a Jansky Fellow. S. E. T. is supported by the NSF under grant AST-0098343. J. H. T. is supported by the NSF under grant AST96-18357. A. W. is supported by the NSF under grant AST-9988217.


text:
# A New Kinetic Equation for Compton Scattering

## 1 Introduction

The most fundamental approach to treating Compton scattering of photons from thermal electrons is to use a Boltzmann-like kinetic equation for the photon distribution function. The basic physics of the Compton process is incorporated into certain scattering functions, or kernel functions, that specify the probability of an initial photon scattering into various frequency and angular ranges. This Boltzmann equation gives highly accurate results, but can also involve a heavy computational burden, especially when the number of scatterings is large.

The Boltzmann equation can be substantially simplified when certain conditions are met. Consider the case where the photons and electrons are both non-relativistic, so that $`h\nu /mc^21`$ and $`kT/mc^21`$. Here $`h`$ is Planck’s constant, $`\nu `$ is the photon frequency, $`m`$ is the electron mass, $`c`$ is the speed of light, $`k`$ is Boltzmann’s constant, and $`T`$ is the electron temperature. In this case the scattering kernel is relatively narrow in frequency. In fact the spreading is due to the thermal Doppler width is of order $`\mathrm{\Delta }\nu \nu \alpha ^{1/2}`$, where,

$$\alpha =\frac{kT}{mc^2}.$$
(1)
When the width of the scattering kernel is small compared to the scale over which the radiation field changes substantially, it is possible to approximate the scattering terms in the Boltzmann equation by a second-order differential operator acting on frequency. Such approximate equations are generally called Fokker-Planck equations, but for Compton scattering the term Kompaneets equation is used, in honor of its originator (Kompaneets, 1957).

In its original form the Kompaneets equation applies strictly only to homogeneous, isotropic radiation fields in which there is time dependence but no spatial transport. However, since the Compton scattering process itself is not very anisotropic, the Kompaneets scattering term is often used for non-isotropic radiation fields, introducing some additional degree of approximation. In this paper, for simplicity, the transport equation is also written without the spatial transport terms, but it should be understood that the scattering term, within the isotropic approximation, is directly applicable to cases involving spatial transport as well.

The Kompaneets equation has been a very useful tool for describing the process of Compton scattering from thermal electrons in astrophysics, when the conditions mentioned above are met. Second-order partial differential equations are much easier to handle numerically than full Boltzmann equations. Within its limitations the simple Kompaneets equation manages to incorporate a number of important physical properties and effects, namely,

1. Conservation of photon number.
2. Detailed balance (correct equilibrium solution).
3. The frequency spreading due to the thermal Doppler effect.
4. The frequency shift due to thermal Doppler effect (inverse Compton effect).
5. The frequency shift due to electron recoil (Compton effect).
6. Stimulated scattering.

The Kompaneets equation incorporates the first two properties exactly. The remaining effects are accurately represented only to lowest order in the parameters $`h\nu /mc^2`$ and $`\alpha =kT/mc^2`$.

For all its many positive features, however, the Kompaneets equation suffers from one major shortcoming: it is unable to treat cases where the radiation fields varies significantly on the scale of Compton frequency shifts, such as can occur in the neighborhood of spectral lines and continuum jumps. The usual approach in such cases has been to revert to the full Boltzmann equation, which contains the detailed scattering kernels, or to use a Monte Carlo method. The simplicity of the second order differential operator is thereby lost.

An entirely different approach to the treatment of Compton scattering was presented by Rybicki & Hummer (1994, hereafter RH), designed to apply to conditions typical in stellar atmospheres. Under these conditions the radiation fields change rapidly over the scale of the electron thermal frequency shift $`\mathrm{\Delta }\nu `$ in the neighborhood of spectral lines and continuum edges. The method of RH handles such rapidly varying radiation fields in a numerically efficient way. An essential feature of the RH method is that it expresses the radiation field as a differential operator acting on the emissivity, not vice versa. This departure from the traditional Fokker-Planck type of operator leads us to call this an inverse operator method.

As originally formulated, the RH method incorporated only numbers 1 and 3 of the above physical effects, namely, photon conservation and the spreading due to the thermal Doppler effect. This was not a serious limitation when applied to normal stellar atmospheres, where these are the predominant effects. However, it seemed desirable for have a method that combined the advantages of the Kompaneets and RH equations, which then would be applicable to a much wider class of problems.

The purpose of this paper is to present such a new kinetic equation for Compton scattering, which (like the Kompaneets equation) incorporates all the physical effects 1–6 above, but which also (like the RH method) has the ability to treat rapidly varying radiation fields. While the idea for this new inverse operator method was motivated by the RH method, the derivation presented here is based on a simple alteration of the Kompaneets equation, which is considerably simpler than the approach used in RH. The new method is only accurate to first order in $`h\nu /mc^2`$ and $`\alpha =kT/mc^2`$, and it assumes that the scattering process can be well approximated by isotropic emission. However, since the Kompaneets equation itself has these same limitations, the new method cannot be judged inferior because of them.

From the standpoint of applications to stellar atmospheres, there are several advantages of the new kinetic equation. Like the Kompaneets equation, it now incorporates the processes of Comptonization of the radiation field and the Compton heating/cooling of the gas, which may be of importance in certain cases involving X-ray irradiation, for example. Another feature is that it satisfies detailed balance, and thus correctly describes the approach to thermal equilibrium. At the same time, like the RH method, it also can handle the distortion of the line profiles due to Doppler broadening due to scattering on the electrons. The numerical implementation of the method involves only a minor modification of the RH method, which essentially maintains the latter’s favorable timing requirements, namely proportional to the first power of the number of frequency points.

For problems other than stellar atmostpheres, there are also advantages of the new method. For example, in studying Comptonization of X-rays, Lightman & Rybicki (1979a, b, 1980) wrote the solution to the transfer problem as a sum of products of probability of scattering $`k`$ times and the $`k`$-th iterated kernel function for Comptonization. As we shall see, the equivalent kernel function for the Kompaneets equation is not wholly positive and is highly singular, being a linear combination derivatives of delta functions up to order two. The “iterated” kernels would involve even higher derivatives, and any solution involving sums over such functions would be entirely unworkable. On the contrary, the kernel functions of the inverse operator method are quite normal functions, and their iterates can be found stably by numerical means, making the method of Lightman & Rybicki more practical.

In §2 the basic properties of the Kompaneets equation will be reviewed and the equations of the new inverse operator method will be derived. It is demonstrated that the inverse operator mehtod, like the Kompaneets equation, satisfies all six of the properties listed above. A comparison with the RH method will be given. In §3 we discuss properties of a rather general Boltzmann equation and show how one may define the equivalent kernels of any approximate version of it. This is then done for the Kompaneets and inverse operator methods. Numerical results are presented for the inverse operator kernel, and its properties, including its accuracy, are discussed. In §4 iterates of the kernel function are defined and discussed, and numerical examples are given. In §5 a useful formula is derived that reduces certain summations occuring in the formalism of Lightman & Rybicki (1979a, b, 1980) to numerically tractable forms. In §6 it is shown how initial value problems can be treated using the inverse operator method. Numerical results show the advantages of the inverse operator method over the Kompaneets method for short times, of the order of the mean scattering time. In §7 a short review is given, and some possibilities for future work are discussed.

## 2 The Inverse Operator Method

For isotropic homogeneous radiation, the kinetic equation for Compton scattering can be written,

$$\frac{1}{c}\frac{n}{t}=n_\mathrm{e}\sigma _\mathrm{T}(n+e).$$
(2)
Here $`t`$ is the time and the photon occupation number $`n=n(x,t)`$ is a function of time and frequency, here given in terms of the scaled frequency $`x=h\nu /kT`$. On the right hand side $`n_\mathrm{e}`$ is the electron density and $`\sigma _\mathrm{T}`$ is the Thomson cross section, which we assume is valid for the energies considered. Within the parenthesis the term $`(n)`$ accounts for extinction and the term $`e=e(x,t)`$ for the scattered emission.

The equation is completed by specifying a form for the scattering term $`e`$. In the well-known Kompaneets equation (Kompaneets, 1957; Rybicki & Lightman, 1979, §7.6) this is given explicitly by,

$$e=n+\alpha \frac{1}{x^2}\frac{}{x}\left[x^4\left(\frac{n}{x}+n+n^2\right)\right].$$
(3)
This expression requires that the photon energy is small compared to the electron rest energy, $`h\nu /mc^21`$, and that the thermal electrons are nonrelativistic, $`\alpha 1`$. In addition, the radiation field $`n(x,t)`$ must vary slowly with frequency $`x`$ over the typical width of the scattering kernel, $`\mathrm{\Delta }xx\alpha ^{1/2}`$.

The equations of the new method are found as a simple alteration of the Kompaneets equation. The reversion of the functional relationship (3) between $`n`$ and $`e`$ is easily done to first order in $`\alpha `$. Since $`e=n`$ to lowest order, this may be used to replace $`n`$ in the term multiplied by $`\alpha `$, giving

$$n=e\alpha \frac{1}{x^2}\frac{}{x}\left[x^4\left(\frac{e}{x}+e+e^2\right)\right].$$
(4)
This implicit equation for $`e`$, combined with equation (2), represents our new formulation of Compton scattering.

The Kompaneets method and the new method differ significantly in how $`e`$ is determined from $`n`$. In the ordinary Kompaneets equation this is done by applying a differential operator to $`n`$, as in equation (3). However, in the new method one solves a differential equation for $`e`$, using equation (4). It might be characterized as applying to the radiation field the inverse of a certain operator related to the Kompaneets operator. For this reason we call this the inverse operator method.

The Kompaneets equation has the advantage that by substitution of equation (3) into equation (2) one obtains a single partial differential equation, which makes it much more amenable to analytic treatment. In the inverse operator method one must deal with two separate equations, (2) and (4). However, we shall show below that for many problems this pair of equations is no more difficult to treat numerically than the Kompaneets equation.

The inverse operator method based on equation (4) has very different mathematical properties than the Kompaneets method. For example, if the radiation field $`n`$ varies rapidly enough, application of the Kompaneets operator can lead to negative emission terms $`e`$. This can be seen easily from the above equations. When the radiation field varies significantly on the scale of the frequency shift $`\mathrm{\Delta }xx\alpha ^{1/2}`$, the second derivative term in the Kompaneets operator become of order $`\alpha ^1n`$, so that the second term in equation (3) is no longer of order $`\alpha `$, but is of order unity. The derivative terms, which can take either sign, are able to dominate the first, strictly positive, term. As we shall see below, the inverse operator method does not lead to such negative values of $`e`$.

The manner of derivation of equation (4) automatically guarantees that it has the same order of accuracy in $`\alpha `$ as the Kompaneets equation, and thus satisfies all the properties 3–6 above to order $`\alpha `$. It can easily be shown that it exactly satisfies the first two properties, namely, photon conservation and detailed balance.

To prove photon conservation, we multiply equation (4) by $`x^2`$ and integrate over all $`x`$. The integrated derivative term vanishes at the endpoints, which gives,

$$_0^{\mathrm{}}x^2n𝑑x=_0^{\mathrm{}}x^2e𝑑x.$$
(5)
The photon density is proportional to

$$𝒩=_0^{\mathrm{}}x^2n𝑑x.$$
(6)
Multiplying equation (2) by $`x^2`$ and integrating over all $`x`$, and using equation (5), then yields the photon conservation result, $`d𝒩/dt=0`$.

Property 2, having the exact equilibrium solution, can be proved as follows. In equilibrium, $`n/t=0`$, so by equation (3) we have $`n(x)=e(x)`$. Substituting this into equation (4) and integrating once, using the fact that $`e(x)0`$ strongly as $`x\mathrm{}`$, we find,

$$\frac{n}{x}+n+n^2=0.$$
(7)
The general solution to this is the general Bose-Einstein distribution,

$$n(x)=\frac{1}{e^{\gamma +x}1},$$
(8)
where the constant $`\gamma `$ is related to the chemical potential. This is the correct photon distribution for thermal equilibrium.

Thus we have shown that the inverse operator equation (4) has the same general properties 1–6 as the Kompaneets equation. Since the two equations differ only at a order of approximation beyond the validity of the Kompaneets equation itself, the inverse operator equation cannot be regarded a priori as inferior to the Kompaneets equation. As we shall show, the inverse operator equation is superior in some important respects.

We conclude this section by making contact with RH. In that work the mean intensity $`J(\nu )`$ and emission integral $`E(\nu )`$ based on the usual specific intensity were used. These are related to the quantities $`n`$ and $`e`$ by,

$$J(\nu )=\frac{2h\nu ^3}{c^2}n(x),E(\nu )=\frac{2h\nu ^3}{c^2}e(x).$$
(9)
In terms of these variables, (4) becomes,

$$\alpha \nu \frac{}{\nu }\left[\nu \frac{E}{\nu }+\left(\frac{h\nu }{kT}3\right)E+\frac{c^2E^2}{2kT\nu ^2}\right]+E(\nu )=J(\nu ).$$
(10)
It is instructive to compare this to equation (25) of RH for $`N=1`$, which in the present notation is

$$\alpha \nu \frac{}{\nu }\left[\nu \frac{E}{\nu }\right]+E(\nu )=J(\nu ).$$
(11)
The new inverse operator method of equation (10) is seen to be a generalization of the RH method for $`N=1`$ that includes extra terms involving first order derivatives of the radiation field. These extra terms will merely involve redefinitions of the coefficients of the second-order difference equations, so, from the numerical point of view, the new method is easily incoporated into the iterative solution scheme given in RH.

However, now, in addition to the treatment of Doppler broadening by scattering off thermal electrons, all the above Compton processes 1–6, including Compton heating/cooling, will be properly taken into account. The new method will now be applicable to a wider range of stellar atmosphere problems, including Comptonization of X-rays. In cases where the radiation should approach thermal equilibrium, this will occur correctly.

## 3 Associated Kernel Functions

Some understanding of the relationship between the Kompaneets method and the inverse operator method can be gained by comparison of the equivalent kernel functions associated with each method. In this section we shall define these kernel functions and give some of their general properties. Then we shall derive explicit forms for them corresponding to the two methods.

The scattering process can be conveniently characterized by means of its kernel function, also called scattering function or redistribution function. This function determines the probability that an initial photon of given frequency and direction is scattered into a final photon of some other frequency and direction. For the simplest case of a homogeneous, isotropic radiation field (the case treated in this paper), the kinetic equation can be written,

$$\frac{1}{n_\mathrm{e}\sigma _\mathrm{T}c}\frac{n}{t}=\{\left[1+n(x)\right]K(x,x^{})n(x^{})\left[1+n(x^{})\right]\frac{x^2}{x^2}K(x^{},x)n(x)\}𝑑x^{}.$$
(12)
The kernel function is denoted by $`K(x,x^{})`$. The ratio $`x^2/x^2`$, which ensures the conservation law for photons, has its origin in certain phases space factors in the derivation (see, e.g., Rybicki & Lightman 1978). Stimulated scattering is taken into account by the factors $`1+n`$, which give rise to quadratic, as well as linear, terms in the radiation field.

The integral on the right side of equation (12) must vanish when the the radiation field is in thermal equilibrium with the electrons, that is, when it is a Bose-Einstein distribution of the same temperature. The principle of detailed balance makes the stronger statement that, in thermal equilibrium, not only must the integral vanish, but the integrand itself must vanish frequency by frequency. That is,

$$[1+n(x)]K(x,x^{})n(x^{})=[1+n(x^{})]\frac{x^2}{x^2}K(x^{},x)n(x),$$
(13)
for all $`x`$ and $`x^{}`$. Substituting the Bose-Einstein form (8) into this equation leads to the condition for detailed balance,

$$x^2e^xK(x,x^{})=x^2e^x^{}K(x^{},x).$$
(14)
The kernel function does not depend on the radiation field $`n`$, so many of its general properties can be derived by considering radiation fields that are restricted in certain ways. In particular, it is very interesting to consider cases where $`n1`$, so that the quadratic terms in the equation can be ignored. The kinetic equation then becomes,

$$\frac{1}{n_\mathrm{e}\sigma _\mathrm{T}c}\frac{n}{t}=\left[K(x,x^{})n(x^{})\frac{x^2}{x^2}K(x^{},x)n(x)\right]𝑑x^{},$$
(15)
which is linear in the radiation field. Comparing this to equation (2) we have the general normalization property,

$$1=\frac{x^2}{x^2}K(x^{},x)𝑑x^{},$$
(16)
and the identification,

$$e(x)=K(x,x^{})n(x^{})𝑑x^{}.$$
(17)
The kinetic equation can also be written,

$$\frac{1}{n_\mathrm{e}\sigma _\mathrm{T}c}\frac{n}{t}=n+K(x,x^{})n(x^{})𝑑x^{}.$$
(18)
The most significant distinction between the Kompaneets method and the new method presented here is found in the kernel functions associated with each. The kernel functions associated with each method are not immediately obvious, but there is a simple trick for determining them based on equation (17), which applies to the case of no stimulated emission. After first changing the dummy variable of integration to $`x^{\prime \prime }`$ (say), we note that by substituting the special radiation field $`n(x)=\delta (xx^{})`$ into the integral, we obtain the kernel function $`e(x)=K(x,x^{})`$, where $`x^{}`$ is considered as a fixed parameter.

We shall now discuss these kernel functions, first for the Kompaneets equation, then for the new method.

### 3.1 The Kompaneets kernel

Using the trick given at the end of the previous section, the kernel function $`K_\mathrm{K}(x,x^{})`$ corresponding to the Kompaneets equation can be found by substituting $`n(x)=\delta (xx^{})`$ into equation (3) (without the stimulated term), which gives,

$$K_\mathrm{K}(x,x^{})=\delta (xx^{})+\frac{\alpha }{x^2}\frac{}{x}\left[x^4\{\delta ^{}(xx^{})+\delta (xx^{})\}\right].$$
(19)
Performing the differentiation, we can also write the Kompaneets kernel in the form,

$$K_\mathrm{K}(x,x^{})=(1+4\alpha x)\delta (xx^{})+\alpha (x^2+4x)\delta ^{}(xx^{})+\alpha x^2\delta ^{\prime \prime }(xx^{}).$$
(20)
This shows that the Kompaneets kernel is a linear combination of generalized functions, namely delta functions and derivatives of delta functions up to order 2, which keep it concentrated near the point $`x=x^{}`$. In terms of a limiting process used to define the delta function, such derivatives of delta functions take both positive and negative signs, and so in a real practical sense the Kompaneets kernel is not a completely positive function. This can also be verified by evaluating equation (17) for the Kompaneets kernel with a positive, but rapidly varying, function $`n`$, say a sufficiently narrow Gaussian, which will produce a result $`e`$ that can take negative values.

An important property of the kernel is its set of frequency moments, which characterize the frequency of the scattered photon relative to the frequency $`x`$ of the incident photon. These are defined as,

$$(xx^{})^s=(xx^{})^s\frac{x^2}{x^2}K(x,x^{})𝑑x,$$
(21)
where $`s=0`$, $`1`$, $`2`$, $`\mathrm{}`$. The factor $`x^2/x^2`$ is included to convert the kernel $`K`$ into the appropriate kernel for describing photon numbers, rather than occupation number. Substitution of the Kompaneets form for the kernel (19) and integration by parts, one obtains the general result,

$$(xx^{})^s=\delta _{s0}+\alpha (4x^{}x^2)\delta _{s1}+2\alpha x^2\delta _{s2}.$$
(22)
In particular, we have,

$`(xx^{})^0`$ $`=`$ $`1,`$ (23)
$`(xx^{})^1`$ $`=`$ $`\alpha (4x^{}x^2),`$ (24)
$`(xx^{})^2`$ $`=`$ $`2\alpha x^2,`$ (25)
$`(xx^{})^s`$ $`=`$ $`0,\text{for }s3.`$ (26)
The first of these is a simple restatement of the photon conservation law for scattering. The second and third give the results for the mean shift and variance of the emitted frequency relative to the incident frequency, respectively. These values for the moments are, of course, well known, and were used in the derviation of the Kompaneets equation. The above derivation simply verifies that the Kompaneets kernel we have obtained is consistent with these values.

The final result (26) states that all moments of order 3 or greater vanish. This is a consequence of the derivation of the Kompaneets equation, in which a Taylor series in the emitted frequency is only carried out up to second order, implicitly assuming that higher order moments would vanish.

### 3.2 The Kernel for the Inverse Operator Method

We have seen that the Kompaneeets equation has an equivalent kernel that is a generalized function, containing up to second derivatives of delta functions. We now want to find the equivalent kernel for the inverse operator method for treating Comptonization. First of all the equation for the method will be written without the stimulated scattering terms, since, as we have seen, the kernel function can more easily be found for this purely linear equation. Dropping the $`e^2`$ term in equation (4) gives,

$$n=e\alpha \frac{1}{x^2}\frac{}{x}\left[x^4\left(\frac{e}{x}+e\right)\right].$$
(27)
The kernel function determined by this equation will be denoted by $`\stackrel{~}{K}(x,x^{})`$. Applying the same trick as before, we set $`n(x)=\delta (xx^{})`$ and $`e(x)=\stackrel{~}{K}(x,x^{})`$ in equation (27), which yields,

$$\delta (xx^{})=\stackrel{~}{K}\alpha \frac{1}{x^2}\frac{}{x}\left[x^4\left(\frac{\stackrel{~}{K}}{x}+\stackrel{~}{K}\right)\right].$$
(28)
That is, $`\stackrel{~}{K}=\stackrel{~}{K}(x,x^{})`$ is the Greeen’s function of the differential operator on the right hand side of equation (27).

It is actually possible to solve this equation analytically in terms of Bessel functions, but we omit details of this derivation. We only quote one analytic result for very small values of $`\alpha `$. Defining,

$$\mu =\left(\frac{9}{4}+\frac{1}{\alpha }\right)^{1/2},$$
(29)
the kernel is asymptotically,

$$\stackrel{~}{K}(x,x^{})\frac{1}{2\alpha \mu }\left(\frac{x^{}}{x^3}\right)^{1/2}e^{(x^{}x)/2}\left(\frac{x_<}{x_>}\right)^\mu .$$
(30)
Here $`x_>`$ and $`x_<`$ denote the larger and smaller of $`x`$ and $`x^{}`$, respectively. Because $`\mu \alpha ^{1/2}1`$, the last factor dominates the behavior of the kernel, making it very sharply concentrated about $`x=x^{}`$. By virtue of the slowly varying factors $`(x^{}/x^3)^{1/2}\mathrm{exp}[(x^{}x)/2]`$, this approximate form of the kernel still obeys the detailed balance relation (14).

When the kernel is sufficiently sharp, one may set $`x=x^{}`$ in the slowly varying factors. Taking into account the differing notations, the kernel (29) then reduces to the form of the kernel given in equation (17) of RH for $`N=1`$, showing again the relation of the inverse operator method to RH.

It is of great practical interest here, as it was in RH, that the scattering function $`e`$, or the kernel function $`\stackrel{~}{K}`$, is the solution of a second-order differential equation, (27) or (28), respectively, since these are easily solved numerically. The numerical solution begins with the introduction of a frequency grid of $`N_x`$ points. Using an obvious matrix notation over the frequency space, the Kompaneets relation (3) becomes,

$$𝐞=(\mathrm{𝟏}+\alpha 𝐓)𝐧,$$
(31)
where $`\mathrm{𝟏}`$ is the unit matrix and $`𝐓`$ is a tridiagonal matrix formed from coefficients of the discretization of the second-order differential operator.

Similarly, the inverse operator relation (4) becomes,

$$𝐧=(\mathrm{𝟏}\alpha 𝐓)𝐞,$$
(32)
which can be solve for $`𝐞`$ in terms of $`𝐧`$,

$$𝐞=(\mathrm{𝟏}\alpha 𝐓)^1𝐧.$$
(33)
Comparison of equations (31) and (33) should again make clear why we call this an “inverse operator” method. The numerical solution of a tridiagonal system of equations is quite rapid, with timing of order $`N_x`$. This is the same order as the multiplication of a tridiagonal matrix times a vector, so that the numerical burden associated with the new operator for many applications is of the same order as that of the Kompaneets operator.

The inverse operator kernel is similarly found by discretization of equation (28), giving,

$$𝐃=(\mathrm{𝟏}\alpha 𝐓)\stackrel{~}{𝐊}.$$
(34)
Here $`𝐃`$ is a discrete matrix version of the delta function $`\delta (xx_0)`$. This is a diagonal matrix with diagonal elements that depend on the quadrature scheme to evaluate integrations over $`x`$ (for the crudest scheme these are equal to the inverse of the local frequency differences). Equation (34) can be solved column by column in order $`N_x^2`$ operations to give the complete inverse,

$$\stackrel{~}{𝐊}=(\mathrm{𝟏}\alpha 𝐓)^1𝐃.$$
(35)
If the kernel is only wanted for a few values of $`x_0`$, the number of operations is correspondingly smaller.

In figure 1 results for the inverse operator kernel function $`\stackrel{~}{K}(x,x_0)`$, obtained numerically in the above way, are plotted as the solid curves for $`\alpha =10^3`$ and for values of $`x_0=0.1`$, $`1`$, and $`10`$. The normalizing factor $`(x^3/x_0^2)`$ is chosen to give equal areas under all the curves. Results for the asymptotic results (29) are also plotted as solid curves, but these are so close as to be indistinguishable. As expected, the kernel is skewed towards larger (smaller) values of $`x`$ for $`x_01`$ ($`x_01`$), expressing the tendency for the radiation to approach thermal equilibrium. This can be seen most clearly for the $`x_0=10`$, less so for $`x_0=0.1`$.

In order to get some idea of the absolute accuracy of the inverse operator results, we compare them to an expansion of the true kernel to third order in $`\alpha `$ given by equation (19) of Sazonov & Sunyaev (2000). These results are plotted as dashed curves in figure 1. These show that the relative errors in the inverse operator are less than 20% in a central region of the kernel about $`x_0`$, which extends out to values of frequency $`x`$ where the kernel itself falls to a few percent of the central value. Beyond this central region the relative errors increase rapidly, since the true kernel decays faster than any power law. However, we remark that the absolute errors outside the central region are actually quite small. More extensive comparisons show that the above errors are characteristic of all the cases for which $`\alpha 10^3`$ and $`x/x_010`$.

One can imagine cases where the large relative errors outside the central region could be a problem, for instance, for accurately describing the wings of a Compton-broadened line far away from line center. However, it is not common for the kernel to appear in such pure form, and more typically the small absolute errors in the kernel will cause lesser problems when integrations over frequency are done; recall that the kernel does have the proper frequency moments up to second order.

The errors in the inverse operator kernel are larger when $`\alpha >10^3`$, where additional relativistic corrections, omitted in the Kompaneets equation, begin to be important. The errors are also larger when $`x_0\alpha ^{1/2}`$, since the dispersion associated with Compton recoil is not treated adequately by the Kompaneets equation, as first discussed by Ross et al. (1978). However, the Kompaneets equation suffers from the same problems in these regimes.

The inverse operator kernel shares some important properties with the exact kernel: It is positive and has roughly the right shape in the central regions, including the characteristic discontinuity of slope for $`x=x_0`$. Given that no physics beyond the Kompaneets equation has been assumed, even the modest accuracy of the inverse operator kernels seems most fortuitous, especially as the Kompaneets kernel itself is singular and not wholly positive.

## 4 Iterated kernels

A concept that often proves useful in practice is that of iterated kernels. The physical interpretation of equation (17) is that the kernel function $`K(x,x^{})`$ describes how photons of initial frequency $`x^{}`$ are redistributed to other frequencies $`x`$ after one scattering. The distribution after two scatterings is described by a function $`K_2(x,x^{})`$, clearly given by,

$$K_2(x,x^{})=K(x,x^{\prime \prime })K(x^{\prime \prime },x^{})𝑑x^{\prime \prime },$$
(36)
and, more generally, after $`k>2`$ scatterings, the distribution is described by a function $`K_k(x,x^{})`$, which is defined by the recurrence,

$$K_k(x,x^{})=K(x,x^{\prime \prime })K_{k1}(x^{\prime \prime },x^{})𝑑x^{\prime \prime }.$$
(37)
For convenience, the definition of iterated kernels may be extended to $`k=1`$ and $`k=0`$ by the formulas,

$$K_1(x,x^{})=K(x,x^{}),K_0(x,x^{})=\delta (xx^{}).$$
(38)
It follows that for any nonnegative $`k`$ and $`l`$,

$$K_{k+l}(x,x^{})=K_k(x,x^{\prime \prime })K_l(x^{\prime \prime },x^{})𝑑x^{\prime \prime }.$$
(39)
These formulas may, in principle, be applied to the Kompaneets kernel (20). However, without working out the details, it is obvious that the $`k`$-th iterated kernel is highly singular, involving derivatives up to order $`2k`$ of delta functions, and is not wholly positive. Consequently, even for moderate values of $`k`$ and relatively smooth functions, finding iterated Kompaneets kernels by the above recurrence relations will be impractical numerically.

For the inverse operator method, the iterated kernels, like the kernel itself, will be ordinary positive functions, which can be easily computed numerically. Using the discretization introduced above, we note that the interated kernel for the inverse operator method can be expressed as the recurrence relation

$$\stackrel{~}{𝐊}_k=(1\alpha 𝐓)^1\stackrel{~}{𝐊}_{k1}.$$
(40)
The application of the inverse of a tridiagonal matrix involves of order $`N_x`$ opperations, so the recurrence relation can be applied using of order $`N_x^2`$ operations. This is an order than the naive $`N_x^3`$ operations if full matrices were involved. If one only wants the $`k`$-th scattered emission,

$$e_k(x)=\stackrel{~}{K}_k(x,x^{})n(x^{})𝑑x^{},$$
(41)
then the iterative step,

$$e_k(x)=\stackrel{~}{K}(x,x^{})e_{k1}(x^{})𝑑x^{},$$
(42)
can be formulated, in matrix language, as,

$$𝐞_k=\stackrel{~}{𝐊}𝐞_{k1}=(1\alpha 𝐓)^1𝐞_{k1},$$
(43)
which can be computed in order $`N_x`$ operations.

The iterates of the kernels in figure 1 up to $`k=5`$ have been computed in this fashion and are shown in figure 2. The discontinuity in slope of the kernel itself disappears after one iteration, and the higher iterates become broader and quite smooth. The bias towards approach to thermal equilibrium can be seen in the tendency for the higher iterates to move toward the central Wien frequency of $`x3`$.

Repeated application of equations (40) and (43) lead to the explicit forms,

$$\stackrel{~}{𝐊}_k=(\mathrm{𝟏}\alpha 𝐓)^k𝐃,$$
(44)
and

$$𝐞_k=(\mathrm{𝟏}\alpha 𝐓)^k𝐞_0,$$
(45)
in terms of powers of the matrix $`(\mathrm{𝟏}\alpha 𝐓)^1`$.

## 5 A Summation Formula

Equation (45) is the basis of a formula we have found very useful for the computation of Comptonized spectra using the formalism of Lightman & Rybicki (1979a, b, 1980). Putting that formalism into the present notation, a spectrum $`s(x)`$ is expressed as the sum,

$$s(x)=\underset{k=0}{\overset{\mathrm{}}{}}p_ke_k(x),$$
(46)
where $`p_k`$ is the probability of a photon in the spectrum having scattered $`k`$ times before escaping. Again, in matrix form,

$$𝐬=\underset{k=0}{\overset{\mathrm{}}{}}p_k𝐞_k.$$
(47)
It turns out that $`p_k`$ can often be accurately expressed as the linear combinations of exponential terms of the form $`Az^k`$, where $`A`$ and $`z>1`$ are constants. We shall derive a formula for the spectrum associated with just one term, since the large-$`k`$ behavior of $`p_k`$ is often dominated by one term; examples of this can be found in Nishimura et al. (1986). In any case, by linear superposition, we can use the formula to find the spectrum due to a sum of such terms.

Thus we consider the exponential expression for $`p_k`$,

$$p_k=Az^k=(1z^1)z^k,$$
(48)
where for convenience we take $`A=(1z^1)`$ so that the $`p_k`$ are normalized,

$$\underset{k=0}{\overset{\mathrm{}}{}}p_k=1.$$
(49)
The summation (47) for this choice of $`p_k`$ and with equation (45) is,

$$𝐬=A\underset{k=0}{\overset{\mathrm{}}{}}\left[z^1(\mathrm{𝟏}\alpha 𝐓)^1\right]^k𝐞_k=A\left[\mathrm{𝟏}z^1(\mathrm{𝟏}\alpha 𝐓)^1\right]^1𝐞_k,$$
(50)
summing the geometric series. A slight rearrangement yields the desired formula,

$$𝐬=A𝐞_0+Az^1(\mathrm{𝟏}\alpha 𝐓z^1)^1𝐞_0.$$
(51)
The second term may be found as the solution of a tridiagonal system in order $`N_x`$ operations, so the computation is very fast.

Examples of spectra computed using our summation formula are given in figure 3 for a source $`e_0(x)`$ of soft photons, emitted in a Wien law of temperature $`T_0=10^4`$ K, which are Comptonized by scattering in a hot medium of temperature $`T=10^7`$ K, making $`\alpha =1.68\times 10^3`$. The choice of normalizing factors here for $`s(x)`$ makes the plot analogous to a $`\mathrm{log}\nu F_\nu `$ vs. $`\mathrm{log}\nu `$ plot, where $`F_\nu `$ is the photon flux. The values of $`z`$ are conveniently characterized by the parameter,

$$y_{}=\alpha /\mathrm{ln}z,$$
(52)
which is a kind of overall “$`y`$-parameter” for the scattering process. To see this, we examine the limit where $`z`$ is close to unity, so that $`\mathrm{ln}zz1`$. But the mean number of scatterings is $`N=_kkp_k=(z1)^1`$ from equation (48), so that in the limit we have,

$$y_{}\alpha N,$$
(53)
which is a typical definition for a $`y`$-parameter; see, e.g., Rybicki & Lightman (1979, Eq. 7.41a). In figure 3, as the values of $`y_{}`$ range from 0.03 to 10, the spectrum changes from being nearly a Wien law at the initial $`10^4`$ K temperature to being nearly a Wien law at the temperature $`10^7`$ K of the Comptonizing medium. For values of frequency between these two Wien laws, one sees an approximate power law, which is originially steep, but which eventually becomes flat.

## 6 Initial value problems

In an initial value problem the radiation field is specified at an initial time, say $`t=0`$, and the problem is to find its subsequent behavior as a function of frequency and time, that is, given $`n(x,0)`$, find $`n(x,t)`$. In this section we shall describe the results of numerical solutions of initial value problems for the Kompaneets and the inverse operator methods for the spatially homogeneous case. This will help clarify some of the distinctions between the two methods.

For simplicity the stimulated scattering process will be neglected, so that the equations are linear. The prototypical initial value problem is to start with an initial radiation field concentrated at a single frequency $`x_0`$, so that $`n(x,0)=\delta (xx_0)`$. The discretized vector version of this initial function will be denoted $`𝐧_0`$. Any initial value problem can then be solved by linear superposition of such solutions.

It is convenient to replace $`t`$ by the Compton $`y`$-parameter,

$$y=\alpha \tau =\left(\frac{kT}{mc^2}\right)n_\mathrm{e}\sigma _\mathrm{T}ct,$$
(54)
which measures time in units of the time for significant Comptonization to take place. Combining equations (2) and (3) and writing the result in matrix form, we then have,

$$\frac{𝐧}{y}=\mathrm{𝐓𝐧}.$$
(55)
This equation can be solved by the well-known Crank-Nicholson method. Introducing a discretization in $`y`$ (time), indicated by a subscript $`j`$, we write,

$$\frac{𝐧_{j+1}𝐧_j}{\mathrm{\Delta }y}=\frac{1}{2}\mathrm{𝐓𝐧}_j+\frac{1}{2}\mathrm{𝐓𝐧}_{j+1}.$$
(56)
That is, we use a differencing scheme that is semi-implicit in time. This can be rearranged into the form,

$$(\mathrm{𝟏}\mathrm{\Delta }_0𝐓)𝐧_{j+1}=(\mathrm{𝟏}+\mathrm{\Delta }_0𝐓)𝐧_j,$$
(57)
where,

$$\mathrm{\Delta }_0=\mathrm{\Delta }y/2.$$
(58)
For the initial value problem we are given the radiation field $`𝐧_0`$ at $`y=0`$, and we construct the $`𝐧`$’s at succeeding times by using (58) as a recurrence relation. Given the solution $`𝐧_j`$ we evaluate the right hand side (this involves multiplying by a tridiagonal matrix) then solving the resulting equation for $`𝐧_{j+1}`$ (this involves solving a linear system with tridiagonal matrix). Both operations involving tridiagonal matrices requires only $`N_x`$ operations, so they are very rapid.

Now let us turn to the initial value problem using the inverse operator method. Using the same discretization, equations (2) and (4) become,

$$\frac{𝐧}{y}=\alpha ^1\left[𝐧+(1\alpha 𝐓)^1𝐧\right]=(\mathrm{𝟏}\alpha 𝐓)^1\mathrm{𝐓𝐧}.$$
(59)
This can be differenced using the Crank-Nicholson method, as before, so that,

$$\frac{𝐧_{j+1}𝐧_j}{\mathrm{\Delta }y}=\frac{1}{2}(\mathrm{𝟏}\alpha 𝐓)^1\mathrm{𝐓𝐧}_j+\frac{1}{2}(\mathrm{𝟏}\alpha 𝐓)^1\mathrm{𝐓𝐧}_{j+1}.$$
(60)
After multiplication by $`(\mathrm{𝟏}\alpha 𝐓)`$, this can be rearranged into the form,

$$(\mathrm{𝟏}\mathrm{\Delta }_+𝐓)𝐧_{j+1}=(\mathrm{𝟏}+\mathrm{\Delta }_{}𝐓)𝐧_j,$$
(61)
where,

$$\mathrm{\Delta }_+=\mathrm{\Delta }y/2+\alpha ,\mathrm{\Delta }_{}=\mathrm{\Delta }y/2\alpha .$$
(62)
Equation (61) is a trivial modification of (57), the only change being the replacement of the two appearances of $`\mathrm{\Delta }_0`$ with $`\mathrm{\Delta }_+`$ and $`\mathrm{\Delta }_{}`$ on the opposite sides of the equation. Thus the numerical method developed for the Kompaneets equation can be trivially adapted to solve the initial value for the inverse operator method, and it will be equally fast.

One distinction between the two methods is now apparent. For the Kompaneets equation, the time variable $`t`$ and $`\alpha `$ can be combined into a single variable, the $`y`$-parameter. This implies that all initial value problems starting from the same initial radiation field will have exactly the same evolution for different values of $`\alpha `$, except for a simple rescaling of the time variable. For the inverse operator method, even using the $`y`$-parameter does not completely eliminate the dependence on $`\alpha `$, as can be seen from equation (59), or from the fact that $`\alpha `$ appears in the definitions of $`\mathrm{\Delta }_+`$ and $`\mathrm{\Delta }_{}`$. This implies that the inverse operator method describes the solution not only on the time for significant Comptonization, but also on the typically much shorter mean free time for scattering.

An example of the numerical solution of initial value problems is shown in Figure 4. Here $`\alpha =10^3`$, and the initial radiation field is a delta-function at the frequency $`x_0=10^2`$. The values of the $`y`$-parameter range from $`0`$ to $`4`$ in steps of $`0.1`$. For the larger values of $`y`$, the solution is seen to be converging to the Wien law, $`(1/2)x^3\mathrm{exp}(x)`$, shown as the dotted curve. Both the Kompaneets and the inverse operator results are plotted in figure 4, but for this set of $`y`$ values the two solutions cannot be distinguished to plotting accuracy.

The discussion above suggests that the distinctions between the two methods will become manifest for times comparable to the mean free scattering time, $`t(n_\mathrm{e}\sigma _\mathrm{T}c)^1`$, that is, when $`\tau 1`$ or $`y\alpha `$.

Figure 5 demonstrates that the behaviors of the two methods at relatively short times are indeed very different. The upper set of curves show the inverse operator results and the lower set the Kompaneets results. To avoid confusing overlapping curves, the times have been divided into two sets, $`\tau =0.25`$ to $`1.5`$ in steps of $`0.25`$ for the left panels (a), and $`\tau =1.5`$ to $`15`$ in steps of $`1.5`$ for the right panels (b). The initial radiation field $`n(x,0)=\delta (xx_0)`$ with $`x_0=10^2`$ has not been plotted here, but would have appeared as a vertical line in all plots.

For the inverse operator method, one component of the solution is clearly a remnant of the localized initial delta function. Immediately outside this, there is an extended component that increases with time until $`\tau 1.5`$, then falls. The extended component originally has a noticeable discontinuity of slope at $`x=x_0`$, but eventually becomes smoother as it broadens.

This behavior of the inverse operator results is easily understood analytically. One shows by direct substitution that the solution to the time dependent equation (18) with $`n(x,0)=\delta (xx_0)`$ can be written in terms of iterated kernels,

$$n(x,\tau )=\underset{k=0}{\overset{\mathrm{}}{}}\frac{\tau ^ke^\tau }{k!}K_k(x,x_0).$$
(63)
For moderate times $`\tau 1`$, we expect only a few terms of this series to contribute, so that

$$n(x,\tau )=e^\tau \delta (xx_0)+\tau e^\tau K(x,x_0)+\mathrm{}.$$
(64)
Thus, for the inverse operator method, we expect an early time behavior of a decaying delta function plus a rising contribution of the noniterated kernel, peaking at $`\tau =1`$, then decaying away. The presence of the delta function expresses the fact that after a time $`\tau `$, some fraction (in fact, $`e^\tau `$) of the photons have not yet scattered; in principle, this delta function never goes away, but in the numerical solution it eventually becomes negligible. The presence of the noniterated kernel can be recognized by the discontinuity in slope at $`x_0`$, but at later times this feature decays and disappears as further terms involving iterated kernels begin to dominate. This is precisely what is seen in the numerical solutions for the inverse operator method in the upper panels of figure 5.

Equation (63) is not a useful representation of the solution of the Kompaneets equation, due to the singular nature of the iterated kernels, so the preceding discussion is not applicable to it. The numerical solution of the Kompaneets shown in the lower panels of figure 5 never exhibits a trace of a delta function for $`\tau >0`$, but rather shows quasi-Gaussian forms, which continually broaden and shift with time. For the largest time ($`\tau =15`$) the results are virtually indistinguishable from those of the inverse operator method.

The coefficient of the iterated kernel in the representation (63), that is, $`\tau ^ke^\tau /k!`$, has the same form as the Poisson probability distribution, and for large $`\tau `$ (or large $`k`$) is sharply peaked about $`k=\tau `$. This suggests the large-$`k`$ asymptotic approximation to the iterated kernels,

$$K_k(x,x_0)n(x,\tau =k),$$
(65)
where $`n(x,\tau )`$ is the initial value problem with $`n(x,0)=\delta (xx_0)`$. Indeed, this is essentially the basis of the method of Sunyaev & Titarchuk (1980), who used it to determine Compton spectra.

For small $`\tau `$ the mixing of the different values of $`k`$ in equation (63) is significant, and the approximation (65) is not so good. Nonetheless, Sunyaev (1980) suggested using equation (65) for $`k=1`$ to define the approximate kernel,

$$K_{\mathrm{SS}}(x,x_0)=n(x,\tau =1),$$
(66)
which he called the “Kompaneets kernel.” This differs from our “Kompaneets kernel” $`K_\mathrm{K}`$, which is the true kernel of the Kompaneets equation, as given in equation (20). The kernel $`K_{\mathrm{SS}}(x,x_0)`$ does have the desired property that the $`k`$-th iterate of it is exactly equal to $`n(x,\tau =k)`$, so it does give asymptotically accurate iterates through equation (65). However, we find that $`K_{\mathrm{SS}}`$ is not as accurate as the inverse operator kernel in the central regions, and it does not have the characteristic discontinuity of slope at $`x=x_0`$. In addition, $`K_{\mathrm{SS}}`$ is actually harder to compute than the inverse operator kernel, since it is defined by an integration of the time-dependent Kompaneets equation over a finite interval of time.

## 7 Summary and Discussion

We have demonstrated a simple modification of the Kompaneets method that maintains most of its desirable properties, but which has much better behavior at very short times, of order of the mean free scattering time. This inverse operator method has an equivalent kernel function that is positive and nonsingular, unlike that of the Kompaneets equation. The kernel function is roughly of the right shape in the central regions with an accuracy of about 20%. We have shown that for many applications the inverse operator method requires no more numerical effort than the Kompaneets equation.

We have shown how the iterated kernels of the inverse operator method can be efficiently computed numerically. A summation formula involving iterated kernels is given that reduces the effort in computing Comptonized spectra using the method of Lightman & Rybicki (1979a, b, 1980). We have shown that initial value problems can be solved using the inverse operator method with no more effort than the Kompaneets equation, but with noticeable improvement in the solution at times comparable to the mean scattering time.

In this paper there has been no attempt to advance beyond the physics of the Kompaneets equation. The main goal has been simply to point out the special advantages that occur when one reverses the role of the second order differential operator that connects the radiation field with the emission coefficient. However, in future work it would be desirable to overcome some of the limitations of the inverse operator method, as others have done for the Kompaneets equation. There does not seem to be any reason to expect any difficulties in extending the present theory to anisotropic scattering or polarization.

It would be desirable to find ways to increase the accuracy of the inverse operator method. Obvious improvements would come with using higher order expansions of the coefficients in the parameters $`\alpha `$ and $`h\nu /mc^2`$, but this will not solve all of the accuracy issues.

More improvement might be gained by expressing the emission coefficient as a sum of $`N>1`$ terms, as in RH. Alternatively, and more in keeping with the present approach, one could include higher order moments and derivatives in a Kramers-Moyal expansion (see, e.g., Riskin, 1984) of the Boltzmann equation, and then use this to form new, higher order, inverse operators.

The present work has concentrated on one particular Fokker-Planck equation, the Kompaneets equation. However, the idea of using inverse operators clearly could be generalized and applied to other Fokker-Planck equations. Since each physical situation has its own special properties, it is not possible to predict whether an inverse operator approach would provide significant advantages in any particular case, but it might be worth investigating. An intriguing possibility is that higher order Kramers-Moyal expansions might have better properties when modified by the use of inverse operators, which perhaps might avoid some of the difficulties known to exist for the traditional expansions beyond the Fokker-Planck approximation (Pawula, 1967).

The author wishes to thank A. Loeb and R. Narayan for helpful discussions and comments. The paper was also improved through the valuable suggestions of S.Y. Sazonov, D.G. Hummer, and especially the referee, I. Hubeny.