text:
# Extensive Air Shower Radio Detection: Recent Results and Outlook 1footnote 11footnote 1Invited talk presented by J. Rosner at RADHEP-2000 Conference, UCLA, Nov. 16–18, 2000, proceedings published by AIP. Enrico Fermi Institute Report No. EFI 2000-57, astro-ph/0101089.

## I Introduction

The observation of the radio-frequency (RF) pulse associated with extensive air showers of cosmic rays has had a long and checkered history. In the present report we describe an attempt to observe such a pulse in conjunction with the Chicago Air Shower Array (CASA) and Michigan Muon Array (MIA) at Dugway, Utah. Only upper limits on a signal have been obtained at present, though we are still processing data and establishing calibrations.

In Section 2 we review the motivation and history of RF pulse detection. Section 3 is devoted to the work at CASA/MIA, while Section 4 deals with some future possibilities, including ones associated with the planned Pierre Auger observatory. We conclude in Section 5. This report is an abbreviated version of a longer one nim .

## II Motivation and history

### II.1 Auxiliary information on shower

Present methods for the detection of an extensive air shower of cosmic rays leave gaps in our information. The height above ground at which showers develop cannot be provided by ground arrays, though stereo detection by air fluorescence detectors is useful. Composition of the primary particles, another unknown, is correlated with shower height, with heavy primaries leading to showers which begin higher in the atmosphere.

Radio detection can help fill such gaps. The electric field associated with a charge $`|e|`$ undergoing an apparent angular acceleration $`\ddot{\theta }`$ is $`||=1.5\times 10^{26}\ddot{\theta }`$ V/m, where time is measured in seconds Allan ; Feyn . For typical showers the charges of radiating particles are expected to be able to act coherently to give a pulse with maximum frequency component $`\nu _{\mathrm{max}}(\mathrm{MHz})10^6/R^2(\mathrm{m})`$, where $`R`$ is the distance of closest approach of the shower axis to the antenna. Showers originating higher in the atmosphere are expected to have higher-frequency components. Thus RF detection may be able to add information on shower height and primary composition, and to provide a low-cost auxiliary system in projects such as the Pierre Auger array Auger .

### II.2 Pulse generation mechanisms

Several possibilities have been discussed for generation of a pulse by air showers. Cosmic rays could induce the atmosphere to act as a giant spark chamber, triggering discharges of the ambient field gradient RRW . Compton scattering and knock-on electrons can give rise to a negative charge excess of some 10 to 25% at shower maximum Ask . Separation of positive and negative charges can occur in the Earth’s magnetic field as a result of a $`q𝐯\times 𝐁`$ force KL . This last mechanism is thought to be the dominant one accounting for atmospheric pulses with frequencies in the 30–100 MHz range Allan , and will be taken as the model for the signal for which the search was undertaken. The charge-excess mechanism is probably the major source of an RF signal in a dense material such as polar ice ZHS , but is expected to be less important in the atmosphere.

### II.3 Early measurements

The first claim for detection of the charge-separation mechanism utilized narrow-band techniques at 44 and 70 MHz Jelley ; Weekes . A Soviet group reported signals at 30 MHz Sov , while a University of Michigan group at the BASJE Cosmic Ray Station on Mt. Chacaltaya, Bolivia Chac  studied pulses in the 40–90 MHz range.

The collaboration of H. R. Allan et al. Allan  at Haverah Park in England studied the dependence of signals on primary energy $`E_p`$, perpendicular distance $`R`$ of closest approach of the shower core, zenith angle $`\theta `$, and angle $`\alpha `$ between the shower axis and the magnetic field vector. Their results indicated that the electric field strength per unit of frequency, $`_\nu `$, could be expressed as

$$_\nu =s\frac{E_p}{10^{17}\mathrm{eV}}\mathrm{sin}\alpha \mathrm{cos}\theta \mathrm{exp}\left(\frac{R}{R_0(\nu ,\theta )}\right)\mu \mathrm{V}m^1\mathrm{MHz}^1,$$
(1)
where $`R_0`$ is an increasing function of $`\theta `$, equal (for example) to $`(110\pm 10)`$ m for $`\nu =55`$ MHz and $`\theta <35^{}`$. The constant $`s`$ was originally claimed to be 20. The Haverah Park observations were recalibrated to yield $`s=1.6`$ (0.6 $`\mu `$V m-1 MHz-1 for a $`10^{17}`$ eV shower at $`R=100`$ m) while observations in the U.S.S.R. gave $`s=9.2`$ (3.4 $`\mu `$V m-1 MHz-1 at $`R=100`$ m) Atrash . To estimate the corresponding signal strength in Jelley ; Weekes , we note that the signal power for showers of average primary energy $`E=5\times 10^{16}`$ eV was measured to be about 4 times that of galactic noise, for which Allan  $`_\nu ^{\mathrm{Gal}}1`$$`2\mu `$V m-1 MHz-1. Thus, for such showers, one expects $`_\nu 2`$$`4\mu `$V m-1 MHz-1. Similarly, the estimate Weekes  of an average pulse power $`V_{\mathrm{peak}}^2/2R=10^{12}`$ W gives $`V_{\mathrm{peak}}=10\mu `$V for $`R=50\mathrm{\Omega }`$. Using the relation between pulse voltage and $`_\nu `$ Allan

$$V=30G^{1/2}(\delta \nu /\nu )_\nu ,$$
(2)
where $`G`$ is the antenna gain, and $`\delta \nu `$ is the bandwidth centered at frequency $`\nu `$, we find for an assumed $`G=5`$ (7 dB) (it is not quoted in Ref. Jelley ) and $`\delta \nu /\nu =2.75/44`$ Jelley , a value of $`_\nu 2.4\mu `$V m-1 MHz-1 at a primary energy of $`5\times 10^{16}`$ eV, or about $`5\mu `$V m-1 MHz-1 at $`10^{17}`$ eV if $`_\nu `$ scales linearly with primary energy Allan . For $`G=5`$ the data of Refs. Jelley ; Weekes  thus would favor the higher field-strength claims of the U.S.S.R. group cited in Ref. Atrash .

More recent claims include pulses with components at or below several MHz Agasa ; Yak ; GS ; Gau , and at VHF frequencies GS ; Gau . The Gauhati University group has reviewed evidence for pulses at a wide range of frequencies Gau .

### II.4 Pulse characteristics

The Haverah Park observations are consistent with a model in which the pulse’s onset is generated by the start of the shower at an elevation of about 10 km above sea level, while its end is associated with the greater total path length (shower $`+`$ signal propagation distance) associated with the shower’s absorption about 5 km above sea level. If a vertical shower is observed at a distance of 100 m from its core, the pulse should rise and fall back to zero within about 10 ns, with a subsequent longer-lasting negative component. High frequencies should be less visible far from the shower axis. Heavy primaries should lead to showers originating higher in the atmosphere, with consequent higher-frequency RF components as a result of the geometric aspect ratio with which they are viewed by the antenna, and possibly a greater $`_\nu `$ for a given primary energy Allan . The polarization of the pulse should be dictated by the mechanism of pulse generation: e.g., perpendicular to the line of sight with component along $`𝐯\times 𝐁`$ for the charge-separation mechanism.

### II.5 RF backgrounds

Discharges of atmospheric electricity will be detected at random intervals at a rate depending on local weather conditions and ionospheric reflections. Man-made RF sources include television and radio stations, police and other communications services, broad-band sources (such as ignition noise), and sources within the experiment itself. The propagation of distant noise sources to the receiver is a strong function of frequency and of solar activity.

Galactic noise can be the dominant signal in exceptionally radio-quiet environments for frequencies in the low VHF (30–100 MHz) range Allan . For higher frequencies in such environments, thermal receiver noise becomes the dominant effect.

## III The installation at CASA

The CASA/MIA detector is located about 100 km southwest of Salt Lake City, Utah, at the Dugway Proving Ground CASAnim . The Chicago Air Shower Array (CASA) is a rectangular grid of $`33\times 33`$ stations on on the desert’s surface. The inter-station spacing is 15 m. A station has four 61 cm $`\times `$ 61 cm $`\times `$ 1.27 cm sheets of plastic scintillator each viewed by its own photomultiplier tube (PMT). When a signal appears on 3 of 4 PMTs in a station, a “trigger request pulse” of 5 mA with 5 $`\mu `$s duration is sent to a central trailer, where a decision is made on whether to interrogate all stations for a possible event. Details of this trigger are described in Ref. CASAnim . When this experiment was begun the CASA array had been reconfigured to remove the 4 westernmost “ribs” of the array. For runs performed in 1998, the easternmost rib had also been removed.

The University of Michigan designed and built a muon detection array (MIA) to operate in conjunction with CASA. It consists of sixteen “patches,” each having 64 muon counters, buried 3 m below ground at various locations in the CASA array. Each counter has lateral dimensions 1.9 m $`\times `$ 1.3 m. Four of the patches (numbered 1 through 4), each about 45 m from the center of the array, lie on the corners of a skewed rectangle; four (numbered 5 through 8), each about 110 m from the center of the array, lie on a rectangle with slightly different skewed orientation, and eight (numbered 9 through 16) lie on the sides and corners of a rectangle with sides $`x\pm 180`$ m and $`y\pm 185`$ m, where $`x`$ and $`y`$ denote East and North coordinates.

### III.1 Expected integral rates

The CASA trigger threshold is a few $`\times 10^{14}`$ eV and corresponds to a trigger rate of 10–20 Hz. The expected rates above ($`10^{15},10^{16},10^{17},10^{18}`$) eV are about (1 Hz, 1 per 2 min, 1 per 4 hr, and 2 per mo), respectively, over the 1/4 km2 area of the array. A primary energy of at least $`10^{17}`$ eV seems to be needed if radio signals are to exceed the galactic noise level of 1–2 $`\mu `$V m-1 MHz-1. Such pulses could be generated by a $`10^{17}`$ eV vertical shower with axis 100 m from the antenna under the most optimistic estimates. Since the whole array should see such showers only every few hours, and most have axes farther from the antenna than 100 m, the possibility of accidental noise pulses during such long time intervals reduces the expected sensitivity considerably.

### III.2 The “radio shack” at CASA

A survey of the CASA/MIA site determined that within the array, broad-band noise associated with computers, switching power supplies, and other electronics was so intense that no RF searches could be undertaken. The same was true at any position within the perimeter of the array. Consequently, an antenna was mounted on top of a mobile searchlight tower at a height of 10 m about 24 m east of the eastern edge of the array, corresponding to $`x=263.8`$ m, $`y=0`$ m. The antenna, a 9-element portable log-periodic antenna manufactured by Dorne and Margolin, was acquired from FairRadio Co. in Lima, Ohio, for about $60. Its nominal bandwidth is 26–76 MHz but it was measured to have usable properties up to 170 MHz.

The signal was fed through 60 ft. of RG-58U cable, filtered by a high-pass filter admitting frequencies above 23 MHz, preamplified using a Minicircuits ZFL-500LN preamplifier with 26 dB of gain, low-pass-filtered to admit frequencies below 250 MHz, and fed to the oscilloscope at a sensitivity of 5 mV/division. This constituted the “wide-band” configuration used for most data acquisition runs. The filters were Minicircuits BNC coaxial models. A “narrow-band” configuration with response between 23 and 37 MHz and two preamplifiers had substantially poorer signal-to-noise ratio in distinguishing transient signals from background.

A trigger based on the coincidence of seven of the eight outer muon “patches” was set to select large showers. Each muon patch was set to produce a trigger pulse when at least 5 of its 64 counters registered a minimum-ionizing pulse within $`5.2`$ $`\mu `$s of one another. The pulses were then combined to produce a summed pulse, fed to a discriminator, whose output was amplified and sent over a cable (with measured delay time 2.15 $`\mu `$s) to the RF trailer. No evidence for pickup of the trigger pulse from the antenna was found. The trigger corresponded to a minimum shower energy somewhat below $`10^{16}`$ eV, based on the integral rate REF  at $`10^{18}`$ eV of 0.17/km2/day/sr.

A Tektronix TDS-540B digitizing oscilloscope registered filtered and preamplified RF data on a rolling basis. These data were then captured upon receipt of a large-event trigger and stored on hard disk using a National Instruments GPIB interface. Data were taken using separate computers (at different times), allowing for analysis both at the University of Washington and at Chicago. The Washington system used a Macintosh Quadra 950 running Labview, while the Chicago system used either a Dell XPS200s Pentium desktop computer or a Dell Latitude LM laptop computer running a C program adapted from those provided by National Instruments. Each trigger caused 50 $`\mu `$s of RF data, centered around the trigger and acquired at 1 GSa/s, to be saved.

The total trigger rate ranged between about 20 and 50 events per hour, depending on intermittent noise sources in the trigger system. Concurrently, the CASA on-line data acquisition system was instructed to write files of events in which at least 7 out of the 8 outermost muon patches produced a patch pulse. These files typically overlapped with the records taken at the RF trailer to a good but not perfect extent as a result of occasional noise on the trigger line.

### III.3 Raw data and RF backgrounds

To remove strong Fourier components associated with narrow-band RF signals which were approximately constant over the duration of each data record, a MATLAB routine performed the fast Fourier transform of the signal and renormalized the large Fourier components to a given maximum intensity. Fig. 1 shows the fast Fourier transform of a typical RF signal before and after this procedure was applied. In each case the data were acquired using the “wide-band” filter configuration mentioned above, whose response cuts off sharply below 23 MHz.

The effect of digital filtering on detectability of a transient is illustrated in Fig. 2. The top panel shows the RF record whose Fourier transform was given in Fig. 1, on which has been superposed a simulated transient of peak amplitude 14.5 digitization units. (The data acquisition scale ranges from $`128`$ to $`+127`$ digitization units; one scale division on the oscilloscope corresponds to 25 units.) The transient is invisible beneath the large amplitude associated with television and FM radio signals. The middle panel shows the result after application of the Fourier coefficient shrinkage algorithm. The bottom panel shows the same record after denoising with a wavelet routine DFS . The records in Figs. 1 and 2 were obtained for 32,768 data points obtained at a 1 ns sampling interval, with the trigger at the 20,000th point.

### III.4 Signal simulation

To quantify signal processing efficiency, we generated the expected signal, fed it through the same preamplifier and filter configurations used for data acquisition, and superposed it on records otherwise free of transients. We successively reduced the amplitude of the superposed test signal until it was no longer detectable, thereby obtaining an estimate of sensitivity.

A Hewlett-Packard Arbitrary Waveform Generator was used to generate signals whose typical characteristics are illustrated in Fig. 3. These signals were taken to have the form $`f(t)=\theta (t)At^2(e^{Bt}Ce^{Dt})`$ with the coefficient $`C`$ chosen so that $`f(t)`$ has no DC component, and $`D`$ corresponding to a long duration of the negative-amplitude component. For all pulses we chose $`D=B/20`$, so that $`C=(8000)^1`$ cancels the DC component. The Fourier components of the test pulse fall off smoothly with frequency. The initial $`t^2`$ behavior was chosen so that both the test pulse and its first derivative vanish at $`t=0`$, as might be expected for a pulse from a developing shower. We chose $`B=0.8,0.4,0.2,0.1`$ corresponding to a time difference between pulse onset and maximum of $`\delta =2.5,5,10,20`$ ns and simulated both narrow-band (23–37 MHz) and broad-band (23–250 MHz) configurations.

The shape of the pulse of Fig. 3 is affected by preamplification and filtration as shown in Fig. 4 for the broad-band example. (The narrow-band configuration leads to a longer ringing time.) The noise is associated with the system used to generate the test pulse, and the fact that the Fourier transform is taken over a much longer time than the duration of the pulse. The sharp feature at 125 MHz is a local artifact.

Systematic studies of signal-to-noise ratios have been performed so far only for the simulated pulses with $`\delta =5`$ ns applied to the broad-band front end. A typical pulse of this type gave a front end output of 21 mV peak-to-peak, acquired at an oscilloscope sensitivity of 5 mV per division. Since each division corresponds to 25 digitization units, the peak-to-peak range is about 104 digitization units, or slightly less than half the dynamic range (255 units, or 8 bits). Positive and negative peaks are thus about 52 digitization units each.

The stored test signal is then multiplied by a scale factor and added algebraically to a collection of RF records in which, in general, randomly occurring transients will be present. One then inspects these records to see if the transient can be distinguished from random noise.

For the broad-band data we estimated that pulses with input voltages corresponding to about 1/5 the original test pulse can be distinguished from average noise (not from noise spikes!). Since the original test pulse had a peak value of 1.3 mV, this corresponds to sensitivity to an antenna output of about $`V_{\mathrm{pk}}260\mu `$V. The ability to detect such a pulse with an effective bandwidth of about 30 MHz corresponds to a threshold sensitivity at the level of order $`3\mu `$V/m/MHz nim .

Preliminary studies of simulated pulses applied to the narrow-band front end suggest a considerably poorer achievable signal-to-noise ratio, despite the expectation that the signal should have a large portion of its energy between 23 and 37 MHz. It appears difficult to detect a pulse from the antenna below about 0.7 mV, which for a bandwidth of 14 MHz corresponds to a threshold sensitivity of $`7\mu `$V/m/MHz nim . Studies of possible improvements of the analysis algorithm for the narrow-band data are continuing.

### III.5 Transients detected under various conditions

Several means were used to characterize transients. One method with good time resolution involved the shrinkage of large Fourier coefficients. One can then search for peaks of each data record, plotting their amplitude against time relative to the trigger. One such plot is shown in Fig. 5 for a data run in which CASA HV was delivered to all boxes. A strong accumulation of transients, mostly with amplitude just above the arbitrarily chosen threshold (mean + 3 $`\sigma `$), is visible at times $`5`$ to $`7\mu `$s relative to the trigger. In a comparable plot for a run in which CASA HV was completely disabled (Fig. 6), only a small accumulation at times $`6`$ to $`7\mu `$s is present. This excess appears due to transients with predominantly high-frequency components (over 100 MHz). Since signal pulses are expected to have more power below 100 MHz (see Fig. 4, bottom) we believe that this accumulation is not due to shower radiation, but most likely arises from the muon patches, one of which is within 75 m of the antenna.

A typical transient occurring in a run with CASA HV on is shown in Fig. 7. The transients are highly suppressed (though not in all runs) when CASA boxes within 100 m of the antenna are disabled.

The time distribution of event maxima above an arbitrary threshold for 880 events taken with CASA HV on (one run from January 1998) is shown in the bottom panel of Fig. 5. The mean arrival time is about $`6\mu `$s before the trigger, with a distribution which is slightly broader for pulses arriving earlier than the mean. This broadening may correspond to some jitter in forming the trigger pulse from the sum of muon patch pulses.

As mentioned earlier, the time for the trigger pulse to propagate from the central station to the RF trailer was measured to be $`2.15\mu `$s. One expects a similar travel time for pulses to arrive from muon patches to the central station. Moreover, the muon patch signals are subjected to delays so that they all arrive at the central station at the same time for a vertically incident shower. Thus, the peak in Fig. 5 is consistent with being associated with the initial detection of a shower by CASA boxes. This circumstance was checked by recording CASA trigger request signals simultaneously with other data; they coincide with transients such as those illustrated in Fig. 7 within better than $`1/2\mu `$s.

The RF signals from the shower are expected to arrive around the same time as, or at most several hundred nanoseconds before, the transients associated with CASA operation. They would propagate directly from the shower to the antenna, whereas transients from CASA stations are associated with a slightly longer total path length from the shower via the CASA station to the antenna. There will also be some small delay at a CASA station in forming the trigger request pulse. Thus, we expect a genuine signal also to show up around 6–7 $`\mu `$s before the trigger. However, for data recorded with CASA boxes disabled, no significant peak with the expected frequency spectrum is visible in this time window. The upper limit on the rate of events giving rise to such a peak can be used to set a limit on RF pulses associated with air showers.

In Table 1 we summarize the triggers taken under optimum conditions, corresponding to broad-band data acquisition under conditions of minimum ambient noise. These triggers represent about 50 hours of data. Data were taken with both East-West and North-South antenna polarizations. Since noise from CASA boxes was found to be a significant source of RF transients, data were taken with some or all CASA boxes disabled by turning off high voltage (HV) supply to the photomultipliers.

### III.6 Sensitivity estimate

The main difficulty associated with pulse detection is that signal pulses are not easily distinguishable from large spurious pulses originating from atmospheric discharges. Both air shower pulses and these background noise pulses can considerably exceed the average noise level. We estimate that all signal pulses should arrive between 7 and 6 $`\mu `$s before the trigger, while the time distribution of noise pulses is assumed to be uniform. Hence, a sufficiently large relative accumulation of pulse maxima in that time bin was adopted as a key criterion in search of signal pulses. We also employed several criteria to increase the fraction of signal pulses: The pulse should be larger than some specified magnitude threshold, bandwidth limited and approximately uniform within its limited bandwidth. To date, no significant accumulation has been detected.

At present we are only able, after accounting for the rate of accidental noise pulses, using Monte-Carlo simulation, and taking an antenna gain of $`G=2.5`$ and a cable attenuation factor of 1.4 dB, to set an upper limit of $`s=54`$ in Eq. (1). This is to be compared with the original Haverah Park result $`s=20`$ Allan , the recalibrated result $`s=1.6`$ Atrash , and the Soviet group’s result $`s=9.2`$ Atrash . The noise level at Dugway is too high and the acquired sample is too limited in statistics and dynamic range to allow us to place upper limits strict enough to check the claims of the two groups. We hope that further improvement of the data processing technique will reduce the noise contribution. The magnetic dip angle $`\gamma `$ is much smaller at the Auger site in Argentina ($`34^{}`$ versus $`68^{}`$ at Utah), leading one to expect bigger electric fields for vertical showers and facilitating the detection of shower radiation. We also expect that the Argentina site will be quieter than Dugway and some clarity as regards the calibrating factor will be established.

## IV Outlook

### IV.1 Further possibilities for processing present data

We are still hoping to improve our sensitivity to the point that we can see a true RF signal from a shower even when CASA HV is not disabled. Such a signal should precede CASA-related transients by at least the delay of formation of a phototube signal. An important calibration will be achieved if we can determine whether we are sensitive to galactic noise, which may be responsible for the continuum between 23 and 88 MHz visible in Fig. 1.

We are still exploring improved methods for removing constant RF signals from our records. In this respect we are limited by the 8-bit dynamic range of the TDS-540B oscilloscope. Data taken with various configurations of boxes near the antenna disabled may help us to better characterize the CASA-related transients.

The triggered data may be useful in a rather different context. It has been proposed that radar methods be used to detect ion trails associated with extensive air showers PG . In our case we may be able to investigate sudden enhancements of the signals of distant television signals (on Channels 3, 6, 8, and 12) correlated with receipt of a large-event trigger.

### IV.2 Considerations for Auger site

A number of questions have been suggested by the present investigation if RF pulse detection is to be considered as an adjunct to the Auger array.

* How far from the shower axis can antennas detect pulses from showers with energies above $`10^{18}`$ or $`10^{19}`$ eV? The answer determines whether a sparse array (e.g., one with the same density as Auger stations) would be sensitive to the RF pulse.
* What dynamic range for a data-acquisition system is needed so that a transient signal survives digital filtering? Apparently 8-bit range is not enough.
* What RF interference exists at each site? Surveys would be desirable. They could pinpoint not only narrow-band sources due to broadcasting stations but potentially dangerous broadband sources from switching power supplies, computers, etc. It would be best to undertake such studies after prototype systems are in place.
* What power budget would an RF detection system need? Each Auger solar-powered station is limited to a total budget of 10 watts. Presumably an RF system would have to use auxiliary power, particularly for its fast-digitization and memory components.
* What is the envisioned minimum cost per RF station? It would presumably be dominated by the data-acquisition system; the antennae and preamplifiers would probably be cheap by comparison. A preliminary estimate is less than $3K per station nim .

The Southern Hemisphere Auger site is progressing well toward an engineering array of 40 stations, as we have heard at this conference Zas . It is hoped that in a couple of years an investigation at that site of the feasibility of RF pulse detection can be undertaken.

## V Conclusions

No “golden signal” has been seen for an RF transient associated with extensive air showers of cosmic rays at the CASA site. With further processing, the data may permit the setting of useful upper limits on signals relevant to at least some of the previous claims. A number of useful lessons have been learned if a similar technique is to be tried in conjunction with the Auger project.

## VI Acknowledgements

It is a pleasure to thank Mike Cassidy, Jim Cronin, Brian Fick, Lucy Fortson, Joe Fowler, Rachel Gall, Kevin Green, Brian Newport, Rene Ong, Scott Oser, Daniel F. Sullivan, Fritz Toevs, Kort Travis, Augustine Urbas, and John Wilkerson for collaboration and support on various aspects of this experiment. Thanks are also due to Bruce Allen, Dave Besson, Maurice Givens, Peter Gorham, Kenny Gross, Dick Gustafson, Gerard Jendraszkiewicz, Larry Jones, Dave Peterson, John Ralston, Leslie Rosenberg, David Saltzberg, Dave Smith, M. Teshima, and Stephan Wegerich for useful discussions. This work was supported in part by the Enrico Fermi Institute, the Louis Block Fund, and the Physics Department of the University of Chicago and in part by the U. S. Department of Energy under Grant No. DE FG02 90ER40560.


format: text

text:
# The Cosmic Star-Formation History: The UV finds most

## 1 Introduction

The assigned subtitle of this talk was “The UV Takes it All!” The absurdity of that subtitle must have been obvious to many conference participants when they first read the program, and now, after four days of talks, I hope it is obvious to everyone here. The UV cannot possibly “take all” of the cosmic star-formation history. Data from low and high redshift alike overwhelmingly show that rapidly star-forming galaxies emit the bulk of their bolometric luminosities in the far-infrared, not the UV. Some extremely luminous high-redshift galaxies have not been detected in the UV at all (see, e.g., Hughes’ and Blain’s contributions to these proceedings); many others dominate relatively short $`850\mu `$m exposures but are barely detected in even the deepest (rest-frame) UV images. I was asked by the organizers to make a case for studying the distant universe in the rest-frame UV, and I will; but I cannot pretend that UV surveys are the best way to study all star-forming galaxies at high redshift. Some fraction of the stars in the universe formed in extremely dusty galaxies that are best studied in the infrared. Nevertheless the available data—even IR data!—suggest that most stars formed in objects that are easiest to detect in the rest-frame ultraviolet, and so I have modified the subtitle of this talk to a statement I am prepared to defend. UV observations may not detect most of the luminosity emitted by any single galaxy at high redshift, and may not detect some high redshift galaxies at all, but the average high redshift star-forming galaxy is much more easily detected through its UV radiation than its dust emission. In the first half of my talk I will lay out the arguments supporting that statement. In the second half I will discuss the attempts to estimate star-formation rates for UV-selected galaxies after they have been detected. It is here (as we shall see) that dust obscuration poses the greatest problems for UV surveys.

## 2 The UV finds most

Fig. 1 helps illustrate the argument I’m going to make. The plot shows the level of dust obscuration observed among different types of star-forming galaxies in the local universe. The abscissa, $`L_{\mathrm{bol},\mathrm{dust}}+L_{\mathrm{UV}}L_{\mathrm{SFR}}`$, is the sum of the star-formation luminosity that is and that is not absorbed by dust; this provides a rough measure of total star-formation rate. The ordinate, $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}`$, is the ratio of obscured to unobscured star-formation luminosity; it provides a rough measure of the level of dust obscuration. Here $`L_{\mathrm{bol},\mathrm{dust}}1.5L_{\mathrm{FIR}}`$ is an estimate of a galaxy’s bolometric dust luminosity and $`L_{\mathrm{UV}}\lambda l_\lambda `$ evaluated at $`\lambda 1600`$Å. The well known trend of increasing dust obscuration with increasing star-formation rate is obvious. The expression $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}L_{\mathrm{SFR}}`$ (dashed line) provides a crude fit to the data. Some galaxies in the local universe are very heavily obscured ($`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}>100`$), but these galaxies are rare and together they host only $`5`$% of the known local star-formationsm . The vast majority of local star-formation occurs among less luminous and less obscured spiral and starburst galaxies.

Now imagine taking the local universe and placing it at redshift $`z3`$ (say) and observing it at $`850\mu `$m with SCUBA. The only galaxies bright enough to detect would be the heavily dust-obscured ultraluminous infrared galaxies (ULIRGs) with $`100<L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}<5000`$. When analyzing this hypothetical $`850\mu `$m data it would be natural to conclude that all high-redshift star formation occurred in galaxies similar to the ones that were detected, in extremely dusty galaxies that would be difficult to detect in the rest-frame UV. But this conclusion would be wrong. It relies on the false assumption that the comparatively faint galaxies which host most of the star formation would be as dusty as the ultraluminous galaxies that were detected; it neglects the strong correlation of dust opacity and star-formation rate.

The real situation has obvious similarities to this hypothetical case. The typical galaxy detected by SCUBA appears to be heavily obscured by dust, orders of magnitude brighter in the (rest-frame) far-IR than the UV. But only $``$25% of the $`850\mu `$m background is produced by objects brighter than SCUBA’s 2mJy detection limitbcs , and only $``$8% by objects bright enough to be included in the $`5`$mJy samples of bcr  and crlwb —the only SCUBA samples with well established optical counterparts111see also siobk  and iet , which discuss sub-mm samples selected in a different way but with the similarly bright typical $`850\mu `$m fluxes.. The 5mJy sub-mm/radio samples have provided the strongest support for the popular belief that most high-redshift star formation occurred in extremely dusty galaxies that are not detected in UV surveys, but the fraction of high-redshift star formation that occurs in the extremely dusty objects of the 5mJy sub-mm/radio samples is as small as the fraction of local star formation that occurs in LIRGs and ULIRGs. Should we believe that objects in the 5mJy samples are any more representative of typical star-forming galaxies at high redshift than LIRGs and ULIRGs are of star-forming galaxies in the local universe?

The data suggest that they are not. In the remainder of this section I will argue that dust opacity and star-formation rate are as strongly correlated at high redshift as in the local universe, and that this implies that the fainter galaxies responsible for most of the $`850\mu `$m background are much less obscured in the UV than typical objects detected by SCUBA, sufficiently unobscured, in fact, that their UV emission is easier to detect than their dust emission.

The correlation of dust opacity and star-formation rate can be established at $`z1`$ by comparing galaxies selected at $`850\mu `$m with SCUBA to those selected at $`15\mu `$m with ISO. Owing to the strength of the 7.7 and $`8.6\mu `$m PAH emission features, $`z1`$ star-forming galaxies are comparatively bright at $`15\mu `$m and deep $`15\mu `$m ISO images detect dusty $`z1`$ galaxies many times fainter than those detectable with SCUBA. Fig. 2a compares the star-formation rates and dust opacities of five $`850\mu `$m sources at $`z1`$ to those of the 17 ISO $`15\mu `$m sources in the HDF with $`0.8<z<1.2`$. The trend of decreasing dust opacity with decreasing star-formation rate is clear; the comparatively faint ISO sources are nowhere near as heavily dust obscured as the brighter $`5`$mJy $`850\mu `$m sources. Since neither $`15\mu `$m nor $`850\mu `$m observations detect a very large fraction of $`z1`$ galaxies’ bolometric dust luminosities, large multiplicative corrections were required when estimating $`L_{\mathrm{bol},\mathrm{dust}}`$ for this plot from the detected $`15\mu `$m or $`850\mu `$m fluxes. I adopted corrections that assume these $`z1`$ galaxies have dust SEDs similar to those of rapidly star-forming galaxies in the local universe. Details can be found in as ; but the main result, the correlation of star-formation rate and dust opacity, is insensitive to the adopted corrections: altering $`L_{\mathrm{bol},\mathrm{dust}}`$ moves galaxies diagonally on the plot, parallel to the correlation.

At higher redshifts, galaxies’ PAH features are redshifted well outside ISO’s $`15\mu `$m window and another approach is required to measure the dust opacity of galaxies significantly less luminous than the 5mJy SCUBA sources. One approach, adopted by cet00  and cet01 , is to point SCUBA at many galaxies which are too faint to be detected individually, and sum the measured fluxes. The 33 Lyman-break galaxies at $`z3`$ targeted by cet00  and cet01  were found to have a mean $`850\mu `$m flux of $`0.72\pm 0.31`$mJy, corresponding to $`L_{\mathrm{bol},\mathrm{dust}}2.7\times 10^{11}h^2L_{}`$ for the $`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ cosmology adopted throughout. The mean value of $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}`$ for the galaxies in this sample was 13, significantly lower than the values observed among brighter SCUBA sources at similar redshifts (Fig 2b). Since 0.7mJy is close to the typical flux of the sources that dominate the $`850\mu `$m backgroundbkis bcs , one might expect Lyman-break galaxies’ moderate ratio of dust to UV luminosity $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}10`$ to be more representative of average high-redshift star-forming galaxies than the extreme ratios $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}\stackrel{>}{}200`$ of the bright $`850\mu `$m samples. In any case, the $`850\mu `$m observations of Lyman-break galaxies provide further evidence that star-formation rate and dust opacity are correlated at high redshift.

If we accept that dust opacity and star-formation rate are correlated at high redshift in a way similar to what Figs. 2a and b suggest, we can proceed to estimate the characteristic level of dust obscuration among the comparatively faint galaxies that dominate the $`850\mu `$m background and see whether most of them are easier to detect in the sub-mm or the UV. Given a rest-frame UV detection limit of $`f_\nu (1600`$Å$`\times (1+z))>0.2\mu `$Jy ($`m_{\mathrm{AB}}=25.5`$), roughly appropriate for the Lyman-break survey of set99 , and a sub-mm flux limit of $`f_\nu (850\mu \mathrm{m})=2`$mJy, the confusion limit of SCUBA on the JCMT, a galaxy will be easier to detect in the rest-frame UV than the sub-mm if it has $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}\stackrel{<}{}500`$ at $`z1`$ or $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}\stackrel{<}{}80`$ at $`z3`$as . Fig. 2a suggests that $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}=500`$ is characteristic of galaxies with $`L_{\mathrm{SFR}}10^{12}h^2L_{}`$ at $`z1`$. Galaxies with lower star-formation luminosities will be less dust obscured on average and easier to detect in the UV than the sub-mm. $`10^{12}h^2L_{}`$ corresponds to $`f_\nu (850\mu \mathrm{m})3.3`$mJy for $`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and the local dust SED shape assumed throughout; roughly 85% of the $`850\mu `$m background is produced by objects fainter than 3.3mJybcs . This suggests that the majority of star-formation at $`z1`$ occurs in galaxies that are easier to detect in the rest frame UV than the sub-mm.

A similar conclusion follows at $`z3`$. Fig. 2b suggests that $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}=80`$ is characteristic of galaxies with $`L_{\mathrm{SFR}}10^{12}h^2L_{}`$ ($`f_\nu (850\mu \mathrm{m})2.5`$mJy) at $`z3`$. 80% of the $`850\mu `$m background is produced by objects fainter and presumably less obscured.

This completes my simple argument: a correlation between star-formation rate and dust opacity is seen to exist at every redshift we can probe ($`0<z\stackrel{<}{}3`$), and its slope is sufficient to imply that most of the $`850\mu `$m background was produced by galaxies that are easiest to detect in the rest-frame UV. It is an empirical argument based solely upon the observed ratios of dust to UV luminosity among galaxies whose dust and UV emissions have both been detected. It largely sidesteps the long and uncertain chain of reasoning – the incompleteness corrections, the luminosity-function extrapolations, the conversions of fluxes to bolometric luminosities to star-formation rates, and so on (see Carilli’s contribution) – that makes attempts to address the same question through “Madau diagrams” so famously contradictory and unreliable.

## 3 Interpreting UV surveys

The reason UV surveys detect a large fraction of high redshift star formation is that the UV luminosities of star-forming galaxies are largely independent of dust opacity. This is a consequence of the correlation of star-formation rate and dust opacity: the dustiest objects tend to be the most luminous, and the two effects – increased obscuration and increased luminosity – mostly cancel out at $`\lambda 1500`$Å. As can be seen in Fig. 3, observed UV luminosities are similar for galaxies with dust obscurations $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}`$ spanning four orders of magnitude. ULIRGs in the local universe are as bright in the UV as starburst dwarfs; high redshift galaxies that have been detected through their dust radiation are typically as bright in the UV as those that have been detected through their UV radiation.

This fact makes UV surveys of star-forming galaxies easy to construct but hard to interpret. Referring to Fig. 3, one can see that a UV luminosity-limited survey of the local universe would net a tremendously broad range of objects, an indiscriminate haul of ULIRGs and blue compact dwarfs and spirals like the Milky Way. The same appears to be true of luminosity-limited surveys at high redshift. But UV luminosities don’t give us much to distinguish between these vastly different objects: if we were told only that a local galaxy had $`L_{\mathrm{UV}}=10^9h^2L_{}`$, for example, we would have no idea whether it was a ULIRG or a spiral or a starburst dwarf. How can we tell what sort of objects we have detected in a high redshift UV selected survey? How can we estimate any single object’s star-formation rate to within even an order of magnitude? Measuring mid-IR, far-IR, sub-mm, or radio fluxes for the objects would be ideal, but most high-redshift star formation occurs in objects too faint to be detected at any of these longer wavelengths, and in any case the small fields of view of current instruments at mid-IR to mm wavelengths ($`2^{}\times 2^{}`$) are poorly matched to the large fields of view of modern optical instruments (up to $`40^{}\times 40^{}`$). Until the next generation of long wavelength instruments becomes available, much of our knowledge of galaxy formation at high redshift will have to be derived from UV data alone.

Meurer and collaborators have developed one technique for estimating dust luminosities (and hence star-formation rates) from UV observations (e.g. mhc , hereafter MHC). Dustier local starbursts are not fainter in the UV on average, but they are redder, and MHC observed that a galaxies’ redness in the UV was a surprisingly good predictor of its dustiness $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}`$. Among the galaxies in MHC’s sample, a measurement of the UV spectral slope $`\beta `$ (in $`l_\lambda \lambda ^\beta `$ at $`1200\stackrel{<}{}\lambda \stackrel{<}{}2000`$Å) was sufficient to predict $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}`$ to within a factor of 2. This $`\beta `$/far-IR correlation is the basis of many attempts to estimate the star-formation rates of UV-selected high-redshift galaxies.

Do high-redshift galaxies obey MHC’s $`\beta `$/far-IR correlation? The evidence is ambiguous. In some cases the $`\beta `$/far-IR correlation appears to hold. One example is the lensed Lyman-break galaxy SMMJ14011+0252 at $`z=2.565`$. Its observed $`850\mu `$m and 20cm fluxes of $`15\pm 2`$mJy, $`115\pm 30\mu `$Jyiet  agree well with the predictions of $`23_9^{+14}`$ mJy, $`150_{67}^{+95}\mu `$Jy that the $`\beta `$/far-IR correlation impliesas . In this case the total star-formation rate can be deduced with reasonable accuracy from UV observations alone. The same appears to be true for the ten brightest Lyman-break galaxies in the HDF. According to MHC, the $`\beta `$/far-IR and far-IR/radio correlations predict a total 20cm flux for these galaxies of $`105\pm 24\mu `$Jy, nicely consistent with the total measured flux of $`100\pm 33\mu `$Jy. But there are counter-examples as well. Baker and van den Werf discussed two in their talks, the galaxies MS1512+36-cB58 and MS1358+62-G1, where the $`\beta `$/far-IR correlation appears to overpredict the observed $`850\mu `$m flux by an order of magnitude or more. This suggests – worryingly – that failures of the $`\beta `$/far-IR correlation at high redshift are not restricted to the dustiest galaxies ($`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}\stackrel{>}{}100`$) as they appear to be in the local universe (see Meurer’s contribution to these proceedings).

Analysis of the largest available sample of UV-selected high-redshift galaxies with independent constraints on $`L_{\mathrm{bol},\mathrm{dust}}`$ suggests that most high-redshift galaxies obey the $`\beta `$/far-IR correlation better than MS1512+36-cB58 and MS1358+62-G1 but perhaps not as well as SMMJ14011+0252. Figure 4 shows the predicted and observed 20cm fluxes of 69 “Balmer-break” galaxies in the Hubble Deep and Flanking Fields with spectroscopic redshifts $`z=1.00\pm 0.11`$as . The objects that are predicted (via the $`\beta `$/far-IR and far-IR/radio correlations) to have larger 20cm fluxes clearly do on average; the data do not support the view that MHC’s $`\beta `$/far-IR correlation generally misestimates dust luminosities by more than an order of magnitude (cf. van den Werf’s contribution to these proceedings). But nor do they obviously support the view that the $`\beta `$/far-IR correlation can predict far-IR fluxes as accurately at high redshift as in the local universe. The measured 20cm fluxes of the galaxies in the brightest bin of predicted flux are roughly a factor of 3 times lower than expected, for example. Assessing the significance of this shortfall is difficult because of the large error bars in the predicted fluxes – a reasonable fraction of the objects in the brightest predicted bin were probably scattered out of the next brightest bin by various errors – and in any case a significant discrepancy between the observed and predicted fluxes might be due to a failure of the far-IR/20cm rather than the $`\beta `$/far-IR correlation. But it is safe to say that the available data do not inspire absolute confidence in the validity of the $`\beta `$/far-IR correlation at high redshift.

An additional cause for concern is the fact that a red spectral slope $`\beta `$ does not appear to be the most distinctive characteristic of the high-redshift galaxies known to have large dust luminosities. For example, these are the differences in mean color between $`z1`$ galaxies in the HDF that are and are not ISO $`15\mu `$m sources: $`\mathrm{\Delta }(U_nG)=0.11`$, $`\mathrm{\Delta }(G)=0.45`$, $`\mathrm{\Delta }(I)=0.17`$. The SEDs of the dusty and luminous ISO sources differ from those of fainter and less obscured galaxies primarily in the range $`2500\stackrel{<}{}\lambda _{\mathrm{rest}}\stackrel{<}{}3500`$Å, not $`\lambda _{\mathrm{rest}}\stackrel{<}{}2000`$Å where $`\beta `$ is measured. This does not show that MHC’s prescription for dust correction is wrong, only that it can be improved, but it adds to concerns that their prescription may not perform as well at high redshift as in the local universe.

Even if the $`\beta `$/far-IR correlation fails on some objects, however, there is reason to hope that it may perform reasonably well on an ensemble of high redshift galaxies. For example, the correlation implies a mean dust obscuration for $`z3`$ Lyman-break galaxies of $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}8`$, a plausible value that is midway between the dust obscuration observed in local spirals ($`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}5`$) and local UV-selected starbursts ($`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}15`$)as  and that is consistent with the mean obscuration $`L_{\mathrm{bol},\mathrm{dust}}/L_{\mathrm{UV}}13\pm 6`$ derived from sub-mm observations of Lyman-break galaxies cet01 . Moreover many of the major results to emerge from recent analyses of $`850\mu `$m data – the brightness of $`850\mu `$m background, the domination of the $`850\mu `$m background by $`1`$mJy sources, the $``$ thousand-fold increase in number density of ULIRGs from $`z0`$ to $`z\stackrel{>}{}3`$ – could have been predicted from UV observations alone by applying the $`\beta `$/far-IR correlation to UV-selected high-redshift galaxy populationsas . Applying MHC’s $`\beta `$/far-IR correlation to high-redshift galaxy populations at least leads to results that are reasonable, and so it seems sensible to continue using it cautiously until further long-wavelength observations can validate it or suggest better alternatives.

## 4 Conclusions

This conference concluded with a debate on the best way to explore high-redshift star formation. The best way? Understanding the history of galaxy formation is a large and difficult task, a shared enterprise that will benefit from observations at many wavelengths. I am not sure how much we stand to gain from arguing whose contributions are the most important.

Nevertheless I would have felt derelict in my assigned role as defender of the UV if I had not attacked one statement that was repeated several times at the meeting. This was the claim that UV and IR surveys detect “orthogonal” populations of galaxies, that a significant fraction of stars form in a “hidden” population of dusty galaxies that is not detected in UV surveys. I am not aware of any evidence supporting this claim. Where is the dusty hidden population in the local universe? A luminosity-limited UV survey deep enough to detect most spirals and blue starbursts would detect most LIRGs and ULIRGs as well (cf. Fig. 3). Where is it at $`z1`$? Existing UV-selected surveys easily detect most if not all of the $`z1`$ $`15\mu `$m sources detected by ISO. Where at higher redshift? A UV-selected survey would need to reach $`26.5`$ to detect the majority of the $`1800`$Å luminosity density at $`z3`$set99 , and at this magnitude limit most if not all of the $`850\mu `$m sources in the samples of bcr , iet 222including correction for magnification by gravitational lensing, and crlwb  would be detected. It is easy to construct a “Madau diagram” in which galaxies detected at wavelength $`x`$ contribute far more to the comoving star-formation density than galaxies detected at wavelength $`y`$. Data will confess to anything if you torture them enough. The proper approach is not to compare Madau diagrams – not to see who benefited most from the grotesquely compounded uncertainties – but to turn to the galaxies themselves and measure their luminosities at both wavelengths. When this is done (as we have seen) the answer is unequivocal: any UV survey at a given redshift $`z\stackrel{<}{}3`$ deep enough to detect majority of the UV luminosity density will detect the majority of IR-selected galaxies as well.

But there is one sense in which most star formation is hidden from UV surveys: only a small fraction of the energy emitted by massive stars emerges from typical galaxies in the UV. The majority is absorbed by dust and reradiated in the far-IR. Unless there is some way to estimate the total luminosity of a galaxy’s massive stars from the (often) trace amount that is detected in the UV, it will be impossible to estimate star-formation rates for the detected galaxies and the usefulness of UV surveys will be diminished. Sect. 3 discussed one method for estimating far-IR fluxes from far-UV fluxes that is known to work in the local universe. As we saw, this method appears to work within a factor of $`3`$ at high redshift, but it is unclear if it works much better. The available data inspire more hope than confidence.

The strength of UV surveys is that they detect large numbers of high-redshift galaxies, even ones that are intrinsically very faint, in large and representative comoving volumes. The weakness is that star-formation rates are difficult to estimate for the detected galaxies. IR surveys complement them perfectly: star-formation rates can be estimated with reasonable confidence, but only small regions of the sky can be surveyed and only the most luminous sources can be detected. Harnessing the strengths of both types of surveys would lead to a rapid advance in our understanding of galaxy formation at high redshift. One obvious strategy would be to use IR observations of some UV-selected high-redshift galaxies to attempt to derive (or validate) a method for estimating dust luminosities from UV observations alone, a method that can then be applied to the large numbers of UV galaxies which have not been (and often cannot be) detected in the IR. This sort of multiwavelength cooperation will take us much farther than arguments about which wavelength is superior – despite the best efforts of wavelength partisans like the author of this screed.

I would like to thank the organizers for support at this enjoyable meeting. Special thanks are due to my collaborators S. Chapman, M. Dickinson, M. Giavalisco, M. Pettini, A. Shapley, and C. Steidel for their contributions to this work.


text:
# Radiative Transfer in Obliquely Illuminated Accretion Disks

## 1 INTRODUCTION

The structure and dynamics of an accretion disk around a neutron star or a black hole is often strongly influenced by its irradiation from the compact object or the disk itself. In galactic sources, for example, the optical emission (see, e.g., de Jong, van Paradijs, & Augusteijn 1996), the stability (see, e.g., van Paradijs 1996; King, Kolb, & Burderi 1996), and the warping of the accretion disk (Pringle 1996; Wijers & Pringle 1999) are probably determined by X-ray irradiation from the central source. In extragalactic sources, the observed broad iron emission lines (see Nandra 1997 for a review), the weak or absent Lyman edges (Sincell & Krolik 1997; Antonucci 1999), as well as the disk warping (Pringle 1996) may also be caused by X-ray irradiation of the accretion disk.

Calculating the structure even of a non-illuminated accretion disk is not trivial, because of the complicated physical processes involved, such as the non-LTE character of the transport of radiation (see, e.g., Hubeny & Hubeny 1998), and assumptions regarding the vertical profile of the viscous heating (see, e.g., Laor & Netzer 1989). The problem becomes even more difficult when illumination is taken into account, both because additional physical process, such as photoionization (see, e.g., Raymond 1993, Ko & Kallman 1994), become dominant, and because of the multidimensional character of the problem.

The transport of radiation through an accretion disk illuminated by a distant source is in general two-dimensional and not axisymmetric. The spectrum of the reflected radiation, when only Compton scattering is taken into account, has often been solved in terms of the distribution of photon escape probabilities from the reflecting medium (e.g., Lightman & Rybicki 1980) or using Green’s functions (e.g., Poutanen, Nagendra, & Svenson 1996). When absorption and emission processes are taken into account, the problem is often solved assuming plane-parallel symmetry (e.g., Sincell & Krolik 1997), or employing variants of the $`\mathrm{\Lambda }`$-iteration (e.g., Basko, Sunyaev, & Titarchuk 1974; Matt, Fabian, & Ross 1993) and Monte-Carlo methods (e.g., George & Fabian 1991).

Solving directly the radiative transfer equation in two spatial dimensions is challenging, especially when the energies of the photons change by Compton scattering and the problem becomes non-local both in space and in photon energy (see, e.g., Mihalas, Auer, & Mihalas 1978; Auer & Paletou 1994; Dykema, Klein, & Castor 1996; Dullemond & Turolla 2000; Busche & Hillier 2000 for methods of solving problems in two spatial dimensions). In general, a problem that is two-dimensional in coordinate space requires the solution of the radiative transfer equation in five dimensions, two in coordinate space and three in photon momentum space. In the case of illumination of a geometrically thin accretion disk, however, the problem can be simplified significantly. Because the photon mean free paths in the radial and azimuthal directions are much smaller than any characteristic length scale in the accretion flow, a geometrically thin accretion disk can be decomposed into a finite number of plane-parallel, obliquely illuminated slabs. Solving the radiative transfer equation for each slab requires only four dimensions, one in coordinate space and three in photon momentum space. When the angular dependence of the interaction cross section between photons and matter can be expanded into a finite number of Legendre polynomials, the problem can be simplified even further (Chandrasekhar 1960). In this case, the radiative transfer equation is equivalent to a finite number of equations over one dimension in coordinate space and two dimensions in photon momentum space, i.e., to a finite number of radiative transfer problems in one spatial dimension. These equations can be solved in general using existing numerical methods.

In this paper I derive the system of equations (reduced to one spatial dimension) that describe, in an obliquely illuminated slab, isotropic absorption and emission, as well as Compton scattering to first order in $`ϵ/m_\mathrm{e}`$ and $`T_\mathrm{e}/m_\mathrm{e}`$, where $`ϵ`$ is the photon energy and $`T_\mathrm{e}`$ and $`m_\mathrm{e}`$ are the electron temperature and rest mass. I then explore a variant of the Feautrier method for solving the resulting transfer equations, in which the scattering kernels do not have forward-backward symmetry (see Milkey, Shine, & Mihalas 1975). I illustrate the use of the derived equations and method of solution by solving simple problems related to the albedo of a cold disk. More detailed calculations of the coupled radiation and gas properties using the approach described in this paper will be reported elsewhere.

## 2 RADIATIVE TRANSFER

I solve the radiative transfer problem in a plane-parallel slab illuminated by some arbitrary, external source of photons. I describe all physical quantities in this system using an orthogonal, Cartesian reference frame, with its $`z`$-axis parallel to the finite dimension of the slab. I also set $`c=k_\mathrm{B}=h=1`$, where $`c`$ is the speed of light, $`k_\mathrm{B}`$ is Boltzman’s constant, and $`h`$ is Planck’s constant.

I assume that the electrons in the slab have density $`n_\mathrm{e}(z)`$ and temperature $`T_\mathrm{e}(z)m_e`$, where $`m_e`$ is the electron rest mass. I also assume that the illuminating radiation is a parallel pencil beam of net flux $`\pi F`$. Neglecting induced processes, the radiative transfer problem is linear and therefore the effect on any arbitrary illumination pattern can be computed by summing the solutions obtained for each plane-parallel pencil beam. The direction of illumination is described by the vector $`\stackrel{}{r}=\mathrm{sin}\theta _\mathrm{i}\mathrm{sin}\varphi _\mathrm{i}\widehat{x}+\mathrm{sin}\theta _\mathrm{i}\mathrm{cos}\varphi _\mathrm{i}\widehat{y}+\mathrm{cos}\theta _\mathrm{i}\widehat{z}`$, where $`\theta _\mathrm{i}`$ and $`\varphi _\mathrm{i}`$ are the directional angles.

I describe the radiation field in terms of the monochromatic specific intensity $`I(z,\mu ,\varphi ,ϵ)`$, where $`z`$ is the distance from the edge of the slab, $`\mu \mathrm{cos}\theta `$ and $`\varphi `$ are the directional angles of the propagation vector $`\widehat{l}`$, and $`ϵ`$ is the photon energy. Because I study non-polarized radiation, I have suppressed the dependence of the specific intensity on polarization mode. I assume that absorption and emission in the slab are isotropic and denote the absorption coefficient by $`\chi (z,ϵ)`$ and the source function by $`S(z,ϵ)`$.

Keeping only terms to first order in $`ϵ/m_e`$ and $`T_e/m_e`$, the radiative transfer equation that describes absorption, emission, and Compton scattering can be written as (Pomraning 1973)

$`\widehat{l}I(z,\widehat{l},ϵ)`$ $`=`$ $`\chi (z,ϵ)[S(z,ϵ)I(z,\widehat{l},ϵ)]n_e\sigma _\mathrm{T}\left(12{\displaystyle \frac{ϵ}{m_e}}\right)I(z,\widehat{l},ϵ)`$ (1)
$`+n_e\sigma _\mathrm{T}{\displaystyle 𝑑\mathrm{\Omega }^{}\underset{n=0}{\overset{3}{}}\left(\frac{2n+1}{4\pi }\right)P_n(\widehat{l}\widehat{l}^{})S_nI(z,\widehat{l}^{},ϵ)},`$
where $`\sigma _\mathrm{T}`$ is the angle-integrated cross section for Thomson scattering, $`d\mathrm{\Omega }^{}`$ is the solid-angle element around $`\widehat{l}^{}`$, $`P_n(\widehat{l}\widehat{l}^{})`$ is the Legendre polynomial of order $`n`$, and

$`S_0`$ $`=`$ $`1{\displaystyle \frac{ϵ}{m_e}}\left(1ϵ{\displaystyle \frac{}{ϵ}}\right){\displaystyle \frac{T_e}{m_e}}\left(2ϵ{\displaystyle \frac{}{ϵ}}ϵ^2{\displaystyle \frac{^2}{ϵ^2}}\right),`$
$`S_1`$ $`=`$ $`{\displaystyle \frac{2}{5}}\left[{\displaystyle \frac{ϵ}{m_e}}\left(1ϵ{\displaystyle \frac{}{ϵ}}\right){\displaystyle \frac{T_e}{m_e}}\left(12ϵ{\displaystyle \frac{}{ϵ}}+ϵ^2{\displaystyle \frac{^2}{ϵ^2}}\right)\right],`$
$`S_2`$ $`=`$ $`{\displaystyle \frac{1}{10}}\left[1{\displaystyle \frac{ϵ}{m_e}}\left(1ϵ{\displaystyle \frac{}{ϵ}}\right){\displaystyle \frac{T_e}{m_e}}\left(6+2ϵ{\displaystyle \frac{}{ϵ}}ϵ^2{\displaystyle \frac{^2}{ϵ^2}}\right)\right],`$
$`S_3`$ $`=`$ $`{\displaystyle \frac{3}{70}}\left[{\displaystyle \frac{ϵ}{m_e}}\left(1ϵ{\displaystyle \frac{}{ϵ}}\right)+{\displaystyle \frac{T_e}{m_e}}\left(4+2ϵ{\displaystyle \frac{}{ϵ}}ϵ^2{\displaystyle \frac{^2}{ϵ^2}}\right)\right].`$ (2)
Defining the energy-dependent optical depth as

$$d\tau (ϵ)\left[\left(1\frac{2ϵ}{m_e}\right)n_e\sigma _\mathrm{T}+\chi \right]dz,$$
(3)
the relative importance of absorption and scattering as

$$\rho _\mathrm{a}=\frac{\chi }{(12ϵ/m_e)n_e\sigma _\mathrm{T}+\chi }$$
(4)
and

$$\rho _\mathrm{s}=\frac{n_e\sigma _\mathrm{T}}{(12ϵ/m_e)n_e\sigma _\mathrm{T}+\chi },$$
(5)
and the redistribution function in the scattering Kernel as

$$p(\widehat{l}\widehat{l}^{})\underset{n=0}{\overset{3}{}}\varpi _nP_n(\widehat{l}\widehat{l}^{})$$
(6)
with $`\varpi _n=(2n+1)\rho _\mathrm{s}S_n`$, the radiative transfer equation (1) becomes

$`\mu {\displaystyle \frac{}{\tau }}I(\tau ,\mu ,\varphi ,ϵ)`$ $`=`$ $`I(\tau ,\mu ,\varphi ,ϵ)\rho _\mathrm{a}S(\tau ,ϵ){\displaystyle \frac{1}{4\pi }}{\displaystyle 𝑑\mathrm{\Omega }^{}p(\widehat{l}\widehat{l}^{})I(\tau ,\mu ^{},\varphi ^{},ϵ)}.`$ (7)
Equation (7) is a second-order, intergrodifferential equation in a four-dimensional phase space and is equivalent to a system of four equations in a three-dimensional phase space (see Chandrasekhar 1960, §48.1). This is true because the expansion in Legendre polynomials of the redistribution function (6) terminates after the first four terms.

Following Chandrasekhar (1960), I will write the radiative transfer equation (7) in terms of the specific intensity $`I_\mathrm{d}(\tau ,\mu ,\varphi ,ϵ)`$ of the diffuse radiation field, i.e., of the photons that have interacted with the gas at least once. Expanding the specific intensity of the diffuse radiation field as

$$I_\mathrm{d}(\tau ,\mu ,\varphi ,ϵ)=\underset{m=0}{\overset{3}{}}I_\mathrm{d}^m(\tau ,\mu ,ϵ)\mathrm{cos}\left[m(\varphi _\mathrm{i}\varphi )\right],$$
(8)
where $`\varphi _\mathrm{i}`$ defines the plane of illumination, the transfer equation for the diffuse field becomes equivalent to the system of equations (cf. Chandasekhar 1960, §48.1)

$`\mu {\displaystyle \frac{}{\tau }}I_\mathrm{d}^m(\tau ,\mu ,ϵ)`$ $`=`$ $`I_\mathrm{d}^m(\tau ,\mu ,ϵ)\delta _0^m\rho _\mathrm{a}S(\tau ,ϵ){\displaystyle \frac{1}{2}}{\displaystyle \underset{n=m}{\overset{3}{}}}\varpi _n^mP_n^m(\mu ){\displaystyle _1^1}P_n^m(\mu ^{})I_\mathrm{d}^m(\tau ,\mu ^{},ϵ)𝑑\mu ^{}`$ (9)
$`\left[{\displaystyle \frac{1}{4}}(2\delta _0^m){\displaystyle \underset{n=m}{\overset{3}{}}}\varpi _n^m(1)^{m+n}P_n^m(\mu )P_n^m(\mu _\mathrm{i})\right]Fe^{\tau /\mu _\mathrm{i}},`$
where $`\delta _0^m`$ is Kronecker’s delta,

$$\varpi _n^m\varpi _n\frac{(nm)!}{(n+m)!},$$
(10)
and $`P_n^m(\mu )`$ are the associated Legendre’s functions of the first kind defined by

$$P_n^m(\mu )=(1)^m(1\mu ^2)^{m/2}\frac{d^mP_n(\mu )}{d\mu ^m}.$$
(11)
Written in an explicit form, the zeroth-order equation is

$`\mu {\displaystyle \frac{I_\mathrm{d}^0(\mu )}{\tau }}`$ $`=`$ $`I_\mathrm{d}^0(\mu )\left[{\displaystyle \frac{1}{2}}S_0{\displaystyle \frac{5}{8}}S_2(3\mu ^21)\right]\rho _\mathrm{s}{\displaystyle _1^1}I_\mathrm{d}^0(\mu ^{})𝑑\mu ^{}`$ (12)
$`\left[{\displaystyle \frac{3}{2}}S_1\mu {\displaystyle \frac{21}{8}}S_3(5\mu ^33\mu )\right]\rho _\mathrm{s}{\displaystyle _1^1}\mu ^{}I_\mathrm{d}^0(\mu ^{})𝑑\mu ^{}`$
$`{\displaystyle \frac{15}{8}}S_2(3\mu ^21)\rho _\mathrm{s}{\displaystyle _1^1}\mu ^2I_\mathrm{d}^0(\mu ^{})𝑑\mu ^{}{\displaystyle \frac{35}{8}}S_3(5\mu ^33\mu )\rho _\mathrm{s}{\displaystyle _1^1}\mu ^3I_\mathrm{d}^0(\mu ^{})𝑑\mu ^{}`$
$`{\displaystyle \frac{1}{4}}\left[S_03S_1\mu \mu _\mathrm{i}+{\displaystyle \frac{5}{4}}S_2(3\mu ^21)(3\mu _\mathrm{i}^21){\displaystyle \frac{7}{4}}S_3(5\mu ^33\mu )(5\mu _\mathrm{i}^33\mu _\mathrm{i})\right]\rho _\mathrm{s}Fe^{\tau /\mu _\mathrm{i}}.`$
The first-order equation is

$`{\displaystyle \frac{\mu }{1\mu ^2}}{\displaystyle \frac{_\mathrm{d}^1(\mu )}{\tau }}`$ $`=`$ $`{\displaystyle \frac{_\mathrm{d}^1(\mu )}{1\mu ^2}}{\displaystyle \frac{3}{4}}\left[S_1{\displaystyle \frac{21}{8}}S_3(5\mu ^23)\right]\rho _\mathrm{s}{\displaystyle _1^1}_\mathrm{d}^1(\mu ^{})𝑑\mu ^{}`$ (13)
$`{\displaystyle \frac{15}{4}}S_2\mu \rho _\mathrm{s}{\displaystyle _1^1}\mu ^{}_\mathrm{d}^1(\mu ^{})𝑑\mu ^{}{\displaystyle \frac{105}{32}}S_3(5\mu ^23)\rho _\mathrm{s}{\displaystyle _1^1}\mu ^2_\mathrm{d}^1(\mu ^{})𝑑\mu ^{}`$
$`{\displaystyle \frac{3}{4}}[S_1(1\mu _\mathrm{i})^{1/2}5S_2\mu \mu _\mathrm{i}(1\mu _\mathrm{i}^2)^{1/2}`$
$`+{\displaystyle \frac{7}{8}}S_3(5\mu ^23)(5\mu _\mathrm{i}^23)(1\mu _\mathrm{i}^2)^{1/2}]\rho _\mathrm{s}Fe^{\tau /\mu _\mathrm{i}},`$
where $`_\mathrm{d}^1(\mu )(1\mu ^2)^{1/2}I_\mathrm{d}^1(\mu )`$. The second-order equation is

$`\mu {\displaystyle \frac{I_\mathrm{d}^2(\mu )}{\tau }}`$ $`=`$ $`I_\mathrm{d}^2(\mu ){\displaystyle \frac{15}{16}}S_2(1\mu ^2)\rho _\mathrm{s}{\displaystyle _1^1}I_\mathrm{d}^2(\mu ^{})𝑑\mu ^{}`$ (14)
$`{\displaystyle \frac{105}{16}}S_3\mu (1\mu ^2)\rho _\mathrm{s}{\displaystyle _1^1}\mu ^{}I_\mathrm{d}^2(\mu ^{})𝑑\mu ^{}+{\displaystyle \frac{15}{16}}S_2(1\mu ^2)\rho _\mathrm{s}{\displaystyle _1^1}\mu ^2I_\mathrm{d}^2(\mu ^{})𝑑\mu ^{}`$
$`+{\displaystyle \frac{105}{16}}S_3\mu (1\mu ^2)\rho _\mathrm{s}{\displaystyle _1^1}\mu ^3I_\mathrm{d}^2(\mu ^{})𝑑\mu ^{}`$
$`{\displaystyle \frac{15}{16}}\left[S_2(1\mu ^2)(1\mu _\mathrm{i}^2)7S_3\mu \mu _\mathrm{i}(1\mu ^2)(1\mu _\mathrm{i}^2)\right]\rho _\mathrm{s}Fe^{\tau /\mu _\mathrm{i}}.`$
Finally, the third-order equation is

$`\mu {\displaystyle \frac{I_\mathrm{d}^3(\mu )}{\tau }}`$ $`=`$ $`I_\mathrm{d}^3(\mu ){\displaystyle \frac{35}{32}}S_3(1\mu ^2)^{3/2}\rho _\mathrm{s}{\displaystyle _1^1}(1\mu ^2)^{3/2}I_\mathrm{d}^3(\mu ^{})𝑑\mu ^{}`$ (15)
$`{\displaystyle \frac{35}{32}}S_3(1\mu ^2)^{3/2}(1\mu _\mathrm{i}^2)^{3/2}\rho _\mathrm{s}Fe^{\tau /\mu _\mathrm{i}}.`$
Equations (12)–(15) are four, second-order, partial differential equations. Because they describe the evolution of the diffuse radiation field, they can be solved with the following boundary conditions,

$$I_\mathrm{d}^\mathrm{m}(\tau =0,\mu <0,ϵ)=I_\mathrm{d}^\mathrm{m}(\tau _{\mathrm{max}},\mu >0,ϵ)=0,$$
(16)
where $`\tau _{\mathrm{max}}`$ is the total vertical optical depth of the slab, and

$$I_\mathrm{d}^\mathrm{m}(\tau ,\mu ,ϵ=0)=I_\mathrm{d}^\mathrm{m}(\tau ,\mu ,ϵ\mathrm{})=0.$$
(17)
In practice, the vertical optical depth of an accretion disk can be very large and, therefore, the second of boundary conditions (16) can be exchanged with another condition that is easier to handle numerically.

The first two moments of the specific intensity of the diffuse radiation field, which are necessary for calculating the energy and momentum exchange between photons and matter, can be calculated as

$$J_\mathrm{d}(\tau ,ϵ)=\frac{1}{2}_1^1I_\mathrm{d}^0(\tau ,ϵ,\mu )𝑑\mu $$
(18)
and

$$\stackrel{}{H_\mathrm{d}}(\tau ,ϵ)=H_1(\tau ,ϵ)\mathrm{cos}(\varphi _\mathrm{i})\widehat{x}+H_1(\tau ,ϵ)\mathrm{sin}(\varphi _\mathrm{i})\widehat{y}+H_2(\tau ,ϵ)\widehat{z},$$
(19)
where

$$H_1(\tau ,ϵ)=\frac{1}{4}_1^1(1\mu ^2)^{1/2}I_\mathrm{d}^1(\tau ,ϵ,\mu )𝑑\mu $$
(20)
and

$$H_2(\tau ,ϵ)=\frac{1}{4}\underset{m=0}{\overset{3}{}}_1^1\mu I_\mathrm{d}^m(\tau ,ϵ,\mu )𝑑\mu .$$
(21)
## 3 NUMERICAL METHOD

In §2, I showed that the specific radiative transfer problem in two spatial dimensions described by equation (1) has been reduced to four problems in one spatial dimension each, which are easier to solve. However, solving the latter problems still requires special care, for a number of reasons. First, the interaction of the illuminating radiation with the disk material takes place in the outermost layers of the slab, which are optically thin. As a result, the method of solution of the transfer equation must be accurate in the limit of low optical depth. Second, in a typical disk-illumination problem, the interaction of the illuminating radiation with the disk material is dominated by true-absorption at low photon energies but is scattering-dominated at high photon-energies. For this reason, simple $`\mathrm{\Lambda }`$-iteration procedures are not adequate and either accelerated iterative procedures (e.g., the accelerated $`\mathrm{\Lambda }`$-iteration or the method of variable Eddington factors) or other non-iterative procedures (e.g., the Feautrier method) must be employed. However, even the latter methods are not directly applicable to the problem studied here because the redistribution integrals in the right-hand sides of the four transfer equations (12)–(15) are not forward-backward symmetric. In this section I describe a variant of the Feautrier method that has the desired properties for solving the four one-dimensional radiative transfer equations (12)–(15). I follow in general the procedure outlined by Milkey et al. (1975), pointing out the differences that arise from the particular properties of the problem studied here.

I choose as independent variables the Thomson scattering optical depth

$$\tau _{\mathrm{es}}(z)=_0^zn_\mathrm{e}(z)\sigma _\mathrm{T}𝑑z,$$
(22)
which is independent of photon energy, as well as the photon energy $`ϵ`$, and the direction of propagation $`\mu `$. I then write the four one-dimensional radiative transfer equations in the general form

$$\rho _\mathrm{s}\mu \frac{I^m}{\tau _{\mathrm{es}}}=I^m𝒮^{\mathrm{m},\mathrm{s}}𝒮^{\mathrm{m},\mathrm{a}}m=0,..,3,$$
(23)
where $`𝒮^{\mathrm{m},\mathrm{s}}`$ and $`𝒮^{\mathrm{m},\mathrm{a}}`$ are the symmetric and antisymmetric parts of the source functions and scattering integrals in equations (12)–(15), which depend implicitly on the radiation field. Defining the Feautrier variables

$$u^m(\tau _{\mathrm{es}},ϵ,\mu )\frac{1}{2}[I^m(\tau _{\mathrm{es}},ϵ,\mu )+I^m(\tau _{\mathrm{es}},ϵ,\mu )]$$
(24)
and

$$v^m(\tau _{\mathrm{es}},ϵ,\mu )\frac{1}{2}[I^m(\tau _{\mathrm{es}},ϵ,\mu )I^m(\tau _{\mathrm{es}},ϵ,\mu )],$$
(25)
the transfer equation takes the form of the system of equations (Milkey et al. 1975)

$`\rho _\mathrm{s}\mu {\displaystyle \frac{u^m}{\tau _{\mathrm{es}}}}`$ $`=`$ $`v^m𝒮^{\mathrm{m},\mathrm{s}}`$ (26)
$`\rho _\mathrm{s}\mu {\displaystyle \frac{v^m}{\tau _{\mathrm{es}}}}`$ $`=`$ $`u^m𝒮^{\mathrm{m},\mathrm{a}}.`$ (27)
Note, that because of the lack of forward-backward symmetry in the scattering kernels, equations (26)-(27) cannot be combined into a single second-order equation, as in the usual Feautrier method (see, e.g., Mihalas 1978).

I then discretize all quantities over $`N_\mathrm{d}`$ grid points in the variable $`\tau _{\mathrm{es}}`$, $`N_ϵ`$ grid points in photon energy, and $`N_\mu `$ grid points in the direction of propagation $`\mu `$. For simplicity, I use, e.g., $`u_{ij}^m`$ to denote the first Feautrier variable of order $`m`$ of the diffuse radiation field, evaluated on the $`i`$th grid point in the quantity $`\tau _{\mathrm{es}}`$, on the $`k`$th grid point in photon energy, and on the $`l`$th grid point in the direction of propagation $`\mu `$, such that $`j=k+(l1)N_\mu `$. In order to recover the diffusion of photons in energy space because of Compton scattering by thermal electrons, I use a second-order differencing scheme in photon energy. For example, I denote the first and second derivatives of the first Feautrier variable with respect to energy by

$$\frac{u_{ij}^m}{ϵ}|_j=\underset{k=1}{\overset{1}{}}X_{jk}u_{i,j+k}^m$$
(28)
and

$$\frac{^2u_{ij}^m}{ϵ^2}|_j=\underset{k=1}{\overset{1}{}}Y_{jk}u_{i,j+k}^m,$$
(29)
where I have defined the operators

$`X_{j,1}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{ϵ_jϵ_{j1}}}\right)`$ (30)
$`X_{j,0}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{ϵ_jϵ_{j1}}}{\displaystyle \frac{1}{ϵ_{j+1}ϵ_j}}\right)`$ (31)
$`X_{j,+1}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{ϵ_{j+1}ϵ_j}}\right)`$ (32)
and

$`Y_{j,1}`$ $`=`$ $`{\displaystyle \frac{2}{(ϵ_{j+1}ϵ_{j1})(ϵ_jϵ_{j1})}}`$ (33)
$`Y_{j,0}`$ $`=`$ $`2\left({\displaystyle \frac{1}{ϵ_{j+1}ϵ_j}}{\displaystyle \frac{1}{ϵ_jϵ_{j1}}}\right){\displaystyle \frac{1}{ϵ_{j+1}ϵ_{j1}}}`$ (34)
$`Y_{j,+1}`$ $`=`$ $`{\displaystyle \frac{2}{(ϵ_{j+1}ϵ_{j1})(ϵ_{j+1}ϵ_j)}}.`$ (35)
In differencing with respect to the variable $`\tau _{\mathrm{es}}`$, I note that the quantity $`u_{ij}^m`$ is density-like and I, therefore, use a center-differencing scheme for all interior grid points, i.e.,

$$\frac{u_{ij}^m}{\tau _{\mathrm{es}}}|_{i+1/2}=\frac{u_{i+1,j}^mu_{i,j}^m}{\tau _{\mathrm{es}}^{i+1}\tau _{\mathrm{es}}^i}.$$
(36)
On the other hand, the quantity $`v_{ij}^m`$ is flux-like and I therefore use for all interior grid points

$$\frac{v_{ij}^m}{\tau _{\mathrm{es}}}|_i=2\frac{v_{i+1/2,j}^mv_{i1/2,j}^m}{\tau _{\mathrm{es}}^{i+1}\tau _{\mathrm{es}}^{i1}}.$$
(37)
Using the above differencing schemes, the difference equations in all interior grid points become

$`2(\rho _\mathrm{s})_{ij}\mu _j{\displaystyle \frac{1}{\tau _{\mathrm{es}}^{i+1}\tau _{\mathrm{es}}^{i1}}}\left(v_{i+1/2,j}^mv_{i1/2,j}^m\right)`$ $`=`$ $`u_{ij}^m{\displaystyle \underset{k}{}}R_{ijk}^{m,s}u_{ik}^mG_{ij}^{m,s}`$ (38)
$`(\rho _\mathrm{s})_{i+1/2,j}\mu _j{\displaystyle \frac{1}{\tau _{\mathrm{es}}^{i+1}\tau _{\mathrm{es}}^i}}\left(u_{i+1,j}^mu_{ij}^m\right)`$ $`=`$ $`v_{i+1/2,j}^m{\displaystyle \underset{k}{}}R_{i+1/2,jk}^{m,a}v_{i+1/2,k}^mG_{i+1/2,j}^{m,a},`$ (39)
where I have written explicitly the dependence of the source functions on the radiation field as

$`𝒮_{ij}^{m,s}`$ $`=`$ $`{\displaystyle \underset{k}{}}R_{ijk}^{m,s}u_{ik}^m+G_{ij}^{m,s}`$ (40)
$`𝒮_{i+1/2,j}^{m,a}`$ $`=`$ $`{\displaystyle \underset{k}{}}R_{i+1/2,jk}^{m,a}v_{i+1/2,k}^m+G_{i+1/2,j}^{m,a}.`$ (41)
Note here that the above source functions have a different structure than those of Milkey et al. (1975), because of the presence of the antisymmetric term $`G_{i+1/2,j}^a`$ that does not depend on the diffuse radiation field.

At the illuminated surface of the slab, which I denote by $`i=1`$, the boundary condition (16) translates into $`u_{1j}^m=v_{1j}^m`$ and applying equation (38) on the first half of the first grid cell, denoted by $`(1,3/2)`$, I obtain

$$v_{3/2,j}^m=u_{1j}^m+\frac{1}{2\mu _j(\rho _\mathrm{s})_{1j}}(\tau _{\mathrm{es}}^2\tau _{\mathrm{es}}^1)\left(u_{1j}^m\underset{k}{}R_{1jk}^{m,s}u_{1k}^mG_{1j}^{m,s}\right).$$
(42)
At a very large optical depth, which I denote by $`i=N_\mathrm{d}`$, I set the flux equal to zero, i.e., $`v_{N_dj}^\mathrm{m}=0`$, and applying equation (38) on the last half of the last grid cell, denoted by $`(N_\mathrm{d}1/2,N_\mathrm{d})`$, I obtain

$$v_{N_\mathrm{d}1/2,j}^m=\frac{\tau _{\mathrm{es}}^{N_\mathrm{d}1}\tau _{\mathrm{es}}^{N_\mathrm{d}}}{2(\rho _\mathrm{s})_{N_\mathrm{d}1/2,j}\mu _j}\left(u_{N_\mathrm{d},j}^m\underset{k}{}R_{N_\mathrm{d}jk}^{m,s}u_{N_\mathrm{d}k}^mG_{N_\mathrm{d}j}^{m,s}\right).$$
(43)
Note that this boundary condition is different than equation (16). Finally, at the first and last energy grid-points, boundary conditions (17) become simply

$$u_{i,1+(l1)N_\mu }^m=u_{i,N_ϵ+(l1)N_\mu }^m=v_{i,1+(l1)N_\mu }^m=v_{i,N_ϵ+(l1)N_\mu }^m=0.$$
(44)
The system of algebraic equations (39)–(44) can be written in the general matrix form

$`A_i^mu_i^m+B_i^mu_{i+1}^m`$ $`=`$ $`C_i^mv_{i+1/2}^mG_i^{m,a},i=1,\mathrm{},N_d1`$
$`D_i^mv_{i1/2}^m+E_i^mv_{i+1/2}^m`$ $`=`$ $`F_i^mu_i^mG_i^{m,s},i=1,\mathrm{},N_d,`$ (45)
where $`A_i^m`$, $`B_i^m`$, $`D_i^m`$, and $`E_i^m`$ are $`N_ϵN_\mu \times N_ϵN_\mu `$ diagonal matrices, $`C_i^m`$ and $`F_i^m`$ are full $`N_ϵN_\mu \times N_ϵN_\mu `$ matrices, and $`G_i^{m,a}`$ and $`G_i^{m,s}`$ are $`N_ϵN_\mu \times 1`$ vectors. Equations (45) can then be solved recursively from $`i=N_\mathrm{d}`$ to $`i=1`$, using the relations

$`u_i^m`$ $`=`$ $`V_i^mv_{i+1/2}^m+W_i^m`$ (46)
$`v_{i+1/2}^m`$ $`=`$ $`U_{i+1/2}^mu_{i+1}^m+T_{i+1/2}^m,`$ (47)
where

$`V_i^m`$ $`=`$ $`(F_i^mD_i^mU_{i1/2}^m)^1E_i^m`$ (48)
$`W_i^m`$ $`=`$ $`(F_i^mD_i^mU_{i1/2}^m)^1(D_i^mT_{i1/2}^m+G_i^{m,s})`$ (49)
$`U_{i+1/2}^m`$ $`=`$ $`(C_i^mA_i^mV_i^m)^1B_i^m`$ (50)
$`T_{i+1/2}^m`$ $`=`$ $`(C_i^mA_i^mV_i^m)^1(A_i^mW_{i1/2}^m+G_{i1/2}^{m,a}).`$ (51)
Note here a misprint in the elimination scheme of Milkey et al. (1975) as well as a difference with the above scheme that arises because of the presence of the antisymmetric term $`G_{i1/2}^{m,a}`$.

It is important to point out here that the method presented in this section is very efficient since it requires the solution of only four equations of simplicity equal to the minimum required for calculating the interaction of an external radiation field with an accretion disk. Therefore, the computational cost of this method is exactly equal to four times the cost required for solving the simplest problem of normal illumination of a plane-parallel slab. As an example, the solution of a problem with 20 grid points in optical depth, 10 grid points in angle, and 60 grid points in energy (i.e., similar to the one shown in Fig. 1) requires only 2 CPU minutes on a 500 MHz Alpha processor.

## 4 APPLICATION TO ILLUMINATED, COLD ACCRETION DISKS

In this section, I use the method described in §3 for calculating the albedo of an illuminated, cold accretion disk. I consider concentric annuli of the disk, which I approximate by plane-parallel slabs. I take into account Compton scattering as well as bound-free absorption from a cold gas (Morrison & McCammon 1983) and neglect all other radiation processes. When the radiation field is weak, the ratio of the bound-free absorption coefficient to the scattering cross section is independent of the electron density and temperature. As a result, I can solve the radiative transfer problem in terms of the electron-scattering optical depth $`\tau _{\mathrm{es}}`$ without the need to consider the vertical disk structure. When the radiation field is strong, the illumination of the disk will affect both its ionization balance and the absorption coefficients, and hence the overall solution will depend explicitly on the vertical structure of the disk (see, e.g., Nayakshin & Kallman 2001).

In general, the radiative transfer equation (7) is linear in the specific intensity. I can, therefore, write

$$I(\tau ,\mu ,\varphi ,ϵ)=I_\mathrm{v}(\tau ,\mu ,\varphi ,ϵ)+I_\mathrm{i}(\tau ,\mu ,\varphi ,ϵ),$$
(52)
where $`I_\mathrm{v}`$ and $`I_\mathrm{i}`$ describe the radiation field due to viscous heating and illumination respectively, and solve for the two radiation fields independently. Note here that decomposition (52) is only formally valid, since the absorption and emission coefficients, as well as the scattering kernel may depend on the electron density and temperature that are determined by the total radiation field. The model parameters for the calculations include the angle of illumination, $`\theta _\mathrm{i}`$ with respect to the normal and the spectrum of the illuminating radiation, which I set to a power-law of photon index $`\mathrm{\Gamma }`$, i.e., $`F(ϵ)ϵ^{\mathrm{\Gamma }+1}`$.

Figure 1 shows the contribution to the vertical flux (Eq. ) of the various orders of the decomposition (8), for an illumination angle of $`\theta _\mathrm{i}=80^{}`$ and a power-law spectrum of photon index $`\mathrm{\Gamma }=2.5`$. Even though the term of zeroth-order provides the dominant contribution to the vertical flux, the contribution of higher-order terms is not negligible. This is shown in Figure 2, where the contribution of the high-order terms to the vertical flux is plotted for different angles of illumination. The total correction, caused by the obliqueness of the illumination, ranges between $`520`$%. Note here that the relative contributions of the different orders of the specific intensity plotted in Figure 2 correspond to the case of zero electron temperature and increase for finite electron temperatures because of the asymmetry of the Compton terms in the scattering kernel.

In the context of accretion onto galactic compact objects, the albedo $`\eta `$ of an accretion disk is usually defined in terms of the fraction of the illuminating flux that does not heat the disk gas (see, e.g., de Jonk et al. 1996). In this section, I first give the general expression for the albedo of an accretion disk that is illuminated obliquely and then evaluate it for the case of geometrically thin, optically thick, cold accretion disks.

Starting from equation (52) and integrating the transfer equation for $`I_\mathrm{i}`$ over photon energy and directional angle I obtain

$$\stackrel{}{H}_\mathrm{i}=\chi J_\mathrm{i}n_e\sigma _\mathrm{T}\left(4\frac{T_\mathrm{e}}{m_\mathrm{e}}\frac{ϵ}{m_\mathrm{e}}\right)J_\mathrm{i},$$
(53)
where $`J_\mathrm{i}`$ and $`\stackrel{}{H}_\mathrm{i}`$ are the energy-integrated, zeroth and first angular moments of $`I_\mathrm{i}`$ and

$$ϵ\frac{ϵJ_\mathrm{i}(ϵ)𝑑ϵ}{J_\mathrm{i}(ϵ)𝑑ϵ}.$$
(54)
The rate $`Q_\mathrm{i}`$, at which the illuminating radiation heats the gas is therefore

$$Q=4\pi \left[\chi +n_e\sigma _\mathrm{T}\left(\frac{ϵ}{m_\mathrm{e}}4\frac{T_\mathrm{e}}{m_\mathrm{e}}\right)\right]J_\mathrm{i},$$
(55)
and the albedo can be written in terms of the volume integral of $`Q`$ as

$`\eta `$ $`=`$ $`1{\displaystyle \frac{4}{\mu _\mathrm{i}}}{\displaystyle _0^{\tau _{\mathrm{es},\mathrm{max}}}}\left[\chi +n_e\sigma _\mathrm{T}\left({\displaystyle \frac{ϵ}{m_\mathrm{e}}}4{\displaystyle \frac{T_\mathrm{e}}{m_\mathrm{e}}}\right)\right]{\displaystyle \frac{J_\mathrm{i}}{F}}𝑑\tau _{\mathrm{es}}`$ (56)
$`=`$ $`1{\displaystyle \frac{2}{\mu _\mathrm{i}F}}{\displaystyle _0^{\tau _{\mathrm{es},\mathrm{max}}}}𝑑\tau _{\mathrm{es}}\left[\chi +n_e\sigma _\mathrm{T}\left({\displaystyle \frac{ϵ}{m_\mathrm{e}}}4{\displaystyle \frac{T_\mathrm{e}}{m_\mathrm{e}}}\right)\right]{\displaystyle _ϵ}𝑑ϵ{\displaystyle _1^1}𝑑\mu I_\mathrm{d}^0(\tau _{\mathrm{es}},ϵ,\mu ),`$
where $`F`$ is the energy-integrated flux of the illuminating radiation and $`\tau _{\mathrm{es},\mathrm{max}}`$ is the vertical height of the accretion disk. Expression (56) allows the study of the vertical profile of the heating of the accretion disk by the irradiating atmosphere as well as the calculation of its albedo, by solving only the zeroth-order transfer equation for $`I_\mathrm{d}^0(\tau _{\mathrm{es}},ϵ,\mu )`$.

Figure 3 shows the albedo of a cold accretion disk, calculated in the $`220`$ keV energy range, for different power-law spectra and angles of illumination. As expected, for larger photon indices of the illuminating radiation, the fraction of low-energy photons, which are efficiently absorbed, is higher and hence the albedo of the disk is lower. At the same time, as the angle of illumination increases, the photons interact with the electrons in a shallower layer of the accretion disk and therefore have a higher chance of escaping after one interaction, increasing the value of the albedo. The overall effect of our treatment of the obliqueness of irradiation is this systematic increase of the disk albedo with illumination angle, which can be up to a factor of $`2`$ larger that in the case of normal illumination.

## 5 CONCLUSIONS

In this paper, I studied the transfer of radiation in an accretion disk that is obliquely illuminated by an external source of radiation. I showed that the resulting transport problem can be decomposed exactly into four one-dimensional problems, which I solved using a variant of the Feautrier method. I then applied this method in calculating the albedos of cold accretion disks.

The calculated values for the albedos are $`0.30`$, even for the softer spectra and larger illumination angles considered here. These values are small and cannot account, for example, for the observed optical magnitudes of galactic low-mass X-ray binaries, which require albedos in excess of $`0.90`$ (see de Jong et al. 1996). Figure 3 shows that relying on a near-grazing illumination of the accretion disk is not enough to account for the observations. High ionization fractions at the surface layers of the disk, which would reduce the absorption of photons, or even the existence of a highly ionized scattering wind above the accretion disk is probably required for the calculated albedos to reach the high values inferred from observations.

In this study, I assumed for simplicity that all heavy elements in the accretion disk are neutral and, therefore their interaction with the illuminating photons is described by the bound-free opacities of Morrison & McCammon (1983). In reality, however, the heated skin of the accretion disk will be collisionally- and photo-ionized and its vertical ionization and thermal balance will need to be calculated self-consistently with the radiation field (e.g., Raymond 1993; Ko & Kallman 1994). Note, however, that the calculation of both the ionization balance and the radiative equilibrium depend only on the zeroth moment of the specific intensity (see eq.  and ) and, therefore, require the solution of only the zeroth-order transfer equation. As a result, the properties of the disk gas can be calculated exactly in a simple, one-dimensional configuration and the full angular dependence of the radiation field can then be calculated with prescribed gas properties. The results of a self-consistent caclulation of the radiation and gas properties will be reported elsewhere.

I am grateful to G. Rybicki for bringing to my attention the possibility of decomposing a multi-dimensional transfer equation into a small number of one-dimensional equations and for carefully reading the manuscript. I also thank Feryal Özel for many useful discussions, especially on the implementation of the Feautrier method in problems with no forward-backward symmetry. This work was supported by a postdoctoral fellowship of the Smithsonian Institution and also, in part, by NASA.