text:
# EXPERIMENTAL EVIDENCE OF GIANT ELECTRON – GAMMA BURSTS GENERATED BY EXTENSIVE ATMOSPHERIC SHOWERS IN THUNDERCLOUDS

## 1 Introduction

A new physical concept of an avalanche type increase of a number of relativistic electrons in gas under the action of the electric field was proposed by Gurevich, Milikh and Roussel - Dupre . The avalanche can grow in electric field $`EE_c`$. The field $`E_c`$ is almost an order of magnitude less than the threshold electric field of conventional breakdown $`E_{th}`$. The growth of number of electrons with energies $`\epsilon >\epsilon _c0.11`$MeV is determined by the fact that under the action of electric field $`E>E_c`$ fast electrons could become runaways , what means that they are accelerated by electric field $`E`$ as suggested by Wilson . Due to collisions with gas molecules they can generate not only large number of slow thermal electrons, but the new fast electrons having energies $`\epsilon >\epsilon _c`$ as well. Directly this process – acceleration and collisions lead to the avalanche type growth of the number of runaway and thermal electrons, which was called in  ”runaway breakdown” (RB). The detailed kinetic theory of RB was developed in .

In atmosphere, the critical field $`E_c`$ is $`100÷150`$ kV/m and exactly these values of electric field are often observed during thunderstorms ,  . When the electric field in thunderstorm cloud reaches the critical value $`EE_c`$ every cosmic ray secondary electron (its energy $`\epsilon >1`$MeV) initiates a micro runaway breakdown. It serves as a source of intensive ionization of air and manifests itself in a strong amplification of $`X`$ and $`\gamma `$ -rays emission in thundercloud ,,.

This emission was observed at airplanes , baloons , and in the mountain experiments .

Extensive atmospheric showers (EAS) are generated by a high energy $`\epsilon `$ cosmic ray particles. EAS are accompanied by a very strong local growth of cosmic ray secondaries number and $`\gamma `$emission in a high energy range $`\epsilon 10÷100MeV`$. This flux interacts effectively with thundercloud electric field if runaway breakdown conditions are fulfilled (RB-EAS interaction). The energy obtained by fast electrons from thundercloud electric field due to the RB-EAS combined effect can serve for the explosive generation of a very large number of a newborn electrons and gamma quanta (see accompanying theoretical paper ).

The experimental study at the Tien Shan Mountain Scientific Station of this phenomena is our goal.

At the station a special system of gamma-spectrometers exist for ”thunderstorm” measurements. It is based on the Geiger-Muller counters SIG5 and NaJ detectors. The system permit to measure intensity of electrons, X and gamma quanta in wide energy range. A detailed description of ”thunderstorm” detector complex can be found in ,.

For the present work a special air shower trigger array was created. The array system has to fix the EAS coming on an extensive region with high accuracy. The system consists of a set of shower electrons and $`\gamma `$quanta detectors (based on proportional Geiger-Muller counters) widely spread over Tien Shan Station territory.

The RB-EAS combined action lead to a strong radio pulses emitted from thundercloud . A specially constructed radio installation was used to study the short electromagnetic pulses of this radio emission . The installation allows to determine pulse form, its maximal intensity and arrival direction (inclination and azimuth angles). The moments of EAS coming is fixed by piping the trigger signal from shower array to the radio installation. That allows to study the EAS and radio emission simultaneously.

The results of measurements are presented. More then 150 simultaneous EAS-radio events are observed. The analysis of the obtained experimental data allows to establish the existence of a new phenomena – giant electron – gamma bursts generated in thunderclouds.

## 2 The cosmic showers trigger array

The cosmic showers trigger array has to fix coming EAS and inform about it radio installation piping there trigger signal. The array is fixed around Tien - Shan Mountain Scientific Station situated 3340 m high above see level ($`43^o02^{}`$N, $`76^o56^{}`$E)

The basic design principles.

The basis of the trigger array topology is the principle of three-fold signal coincidences from a set of detectors placed in the vertexes of a triangle. The whole system was build as an aggregate of such triangles. The common trigger signal is the logical *OR* sum of signals, generated independently by elementary triangle parts.

For the EAS particles detection in trigger array are used the SI5G type self-quenching ionization counters which permit to have a high amplitude signal at relativistic charged particles and high energy gamma quanta passage. The counters are tolerant to sharp changes of weather conditions. The registration efficiency of SI5G counters for relativistic electrons is high ($`>`$ 80%). For gamma quanta it is much lower and depend on quanta energy. The detailed calibration measurements of SI5G counter efficiency for registration of gamma quanta in our conditions were carried out in . The characteristic efficiency is of the order of 1%. Note, that it could be noticeably increased due to photoelectrons emitted from the aluminum walls of the counter box.

The counters are placed inside the boxes made of a 2 mm thick aluminum, which assures a proper shielding from electromagnetic interference. Each detector box contains 18–20 counters. The sensitive area of each SI5G counter being 0.033 m2, the total area of a detector box is around 0.6 $`m^2`$. The detector is placed in a special little housing construction .

Pulse signals from the all distant detector points are transmitted to the data collection centers through the screened electrical cables. To avoid the damage of semiconductor electronics in the moment of lightning discharge, the cables are conjugated both with the detector and with registering apparatus throughout the vacuum tube based circuits (see ,).

Layout of the shower trigger array.

Current status of the shower trigger array is shown in Fig. 1. The array consists of two independent subsystems marked with Roman numerals on the figure; each subsystem, in its turn, — of a set of elementary triangles.

As it is seen from Fig. 1, the total area being covered with trigger array’s perimeter exceeds 0.1 square kilometer.

The size of the elementary detector triangles was determined by such a way to achieve the rate of triple coincidences about one per 10 seconds. Preliminary measurements have shown that our first triangle ABC having the edge sizes 65 x 75 x 65 m3 satisfy this requirement. This triangle was selected as a basic for the trigger array constructing. Fig. 1(down) represents on a large scale the mutual location of detector points in the trigger subsystems. The numbers shown near the edges of elementary triangles mean their length in meters.

Trigger data requisition.

The shower trigger array (so as the whole ”thunderstorm” detector system) has several data registration centers connected with each other and with the trigger detector subsystems by the means of cable lines. Topology of the data transfer cables is shown in Fig. 1.

One of the registration centers — the main center — has a cable connection also with the system of radio emission registration; it is the same place where the common shower trigger signal is elaborated from coincidence signals of separate elementary triangles. The common trigger signal is transmitted to the radio registration system through the special cable; in the moment of its appearance information about the coincidence type (namely, identifiers of the all detector points whose pulses were coinciding with the trigger) is written by a registration computer. Later these data are used to estimate roughly the size of a registered shower.

All the detector points of subsystem *I* (it is also denoted as *CENTER* in Fig. 1) have a direct connection with the main registration point and their pulses are transmitted to that point immediately. This circumstance permits to analyze the trigger coincidence type in particularly and to define precisely the intensity of *ABC*, *BCE*, *CDE* and *DEF* triggers coming from the *CENTER* trigger subsystem.

The signals of the subsystem *II* are analyzed in the peripheral registration point and just the ready coincidence pulses from its basic triangles *G* and *H* are transmitted to the main center (see Fig. 1).

In the main registration center any combination of the *ABC*, *BCE*, *CDE*, *DEF*, *G*, *H* signals is used to generate a common trigger pulse for radio emission registration system. The trigger signal is formed as a coincidence of signals within $`5\mu `$s duration from every three modulus connected in a triangle (ABC for example). That is why the trigger time accuracy of the array is $`5\mu `$s. Its stability was tested with a precision scintillation counter.

Calibration of the array during quiet time

During quiet time (not thunderstorm days) trigger array was carefully calibrated. It was established that in average the full array was registering a shower every 2.45 second (Fig. 2). Calibration was proved by measurements at triangle ABC done simultaneously by SIG5G and scintillation counter counters. A detailed numerical simulations was fulfilled, which demonstrated a reasonable agreement with observations.

It will be shown below that the work of a trigger system in thunderstorm period (which is of our interest) and in quiet days is significantly different. That is why the detailed description of a quiet time array calibration is omitted here. It will be presented in a separate paper.

## 3 Radio registration system

Receiving antennas

Radio measurements were carried out using specially designed installation (analogous to described in ) for short electromagnetic pulse observations in the frequency range from 0.1 to 30 MHz. This installation contains three spaced receiving antennas connected to the receiving apparatus (central unit). It allows signal arrival angle determination using correlation technique as well as waveform recording with high temporal resolution (16 ns). Its block diagram is shown in the Fig. 3.

Each receiving antenna of the installation (antenna assembly) actually consists of three individual antennas—two mutually perpendicular loop antennas for horizontal magnetic field measurements, and an End-Fed antenna to measure vertical electric field. All three antennas are active and include transistor amplifiers. This allowed to reduce the size of the assembly essentially and to reach uniform gain-frequency characteristics for all antennas in the wide frequency range.

Antenna for magnetic field measurements is a screened vertically sited rectangular loop. Loop and prime amplifier screening is needed to provide antenna sensitivity only to magnetic field, in other case its directional pattern becomes non-symmetric and frequency dependent. To equalize frequency-response of the loop it should be loaded by small resistance. Thus it is attached to the amplifier input through a step-up transformer made on long lines working as current transformer.

Electric antenna is an electrically short rod attached to the matched transistor amplifier. Its frequency-response equalization is achieved using principally capacitate character of the electrically short pole impedance as well as of transistor input impedance in the frequency range concerned. So, antenna–amplifier connection is a frequency independent divider with a ratio of antenna capacity to amplifier capacity as a transfer coefficient. Antenna capacity is increased by a plate tip at the end of the rod to increase the transfer coefficient. The increase of antenna capacity extends the antenna frequency range to lower frequencies also. Low boundary of the antenna frequency range is determined by the response time that is equal to the product of the antenna capacity and an input resistance of the amplifier. High input resistance of the amplifier leads to the growth of low frequency noise and assigns a risk of input transistor damage due to breakdown in strong electric fields. Therefore the input resistance was chosen in accordance with low boundary of the antenna frequency range (100 kHz approximately) to be about 1 MOhm. Breakdown protection of the amplifier is provided in three stages by high-voltage feed-through capacitor, gas-filled excess voltage preventer, and two parallel-opposition semiconductor diodes successively mounted at its input.

Data registration and storage

All three antenna amplifiers outputs from each antenna assembly are connected to the central unit by coaxial cables of equal lengths. Separate cables are used to provide antenna amplifiers power supply through appropriate bypasses. The central unit provides switching and selection up to 4 receiving elements (antennas) most suitable under specific conditions of the experiment. It serves also for preliminary analog data reduction, namely, reduction of the dynamic range of received signals to analog-to-digital converter (ADC) input range to prevent loss of information, and both high-pass and low-pass filtering to decrease low-frequency noise and to prevent spectral aliasing during data sampling.

After the central unit chosen signals are digitized by 4-channel radio receiver based on the personal computer (PC) using two synchronized two-channel ADC boards with sampling rate up to 60 MHz and 256 MB per board internal memory. High sampling rate leads to high data traffic that exceeds the ability of PC bus. Therefore pulse operation mode is needed with fast data storage into internal buffer area following by slower data transfer to the PC initiated by trigger pulse. ADC internal buffer memory allows to record up to 1 second sessions with maximum sampling rate of 60 MHz.

Radio receiver is operated in the pretrigger mode allowing signal recording before and after the trigger pulse. It can be triggered by external TTL signal as well as internally. In the last case triggering occurred when an input signal of one of the ADC channels exceeds a preset threshold.

The whole installation is supplied by UPS (Uninterruptible Power Supply) with enhanced up to several hours off-line operation capability (using external accumulator of higher capacity).

Antenna units (assemblies) are mounted at the Tien Shan Scientific Station at the altitude of about 3340 m above the sea level. They form direct angle triangle with perpendicular bases of about 78 m in the $`350^{}`$ azimuth direction with elevation angle of $`7^{}30^{^{}}`$ (below the horizon) and about 74 m in the $`80^{}`$ azimuth direction with elevation angle of $`2^{}25^{^{}}`$ (above the horizon).

An additional receiver operated at about 250 MHz is used as a separate facility to observe VHF radiation from the lightning. This receiver has filter bandwidth of 6 MHz, its output is detected with time constant of 100 $`\mu `$s or 10 $`\mu `$s providing correspondent time resolution. Its data acquisition system can also record EAS trigger pulses and static electric field using additional channels.

## 4 Results of measurements

Quiet time

Observations were carried out from July 19 till October 09, 2003 almost continuously. The radio receiving system was triggered by the EAS facility. Mean time interval between triggers was about 10 s in the beginning of observations decreasing later to 2.45 s due to upgrade of EAS facility and enhancement of its effective area (Fig. 2). Recording time was slightly above 100 $`\mu `$s with 83.2 $`\mu `$s pre-history. Data sampling rate of 60 MHz was used. Short radio pulses with few hundred nanosecond length correspondent to EAS’s were searched. Wide band radio interferometer in frequency range 0.1 – 30 MHz was used in our experiment. The electromagnetic noise in this frequency range is always rather strong. The radio noise spectrum has well defined maximum of the order or less than 0.2 MHz. This frequency is small compared with the frequency of registered signals ($`3`$ MHz). That’s why only noise intensity is significant. It is shown in Fig. 4. Strong changes of noise with time is seen

We note, that the mean electric field amplitude value (in mv/m) approximately 1.5 times larger than the potential (in mv) at antenna rode, shown at the figure.

The general result of the observations is that the short radio pulses under quiet conditions are absent. An example of this type regular record in quiet time is shown at Fig. 5. One can see, that the trigger electric circuit signals lead to the noise amplification around trigger point 82.2$`\mu s`$.

Far thunderstorm

Few records were done for far thunderstorms in the autonomous mode with triggering from the output of one of the receivers as it was described in . The record length was 750 ms with 500 ms pre-history. The same behavior of the radio signal during lightning initiation was observed as it was reported in  for lightning in the middle Russia far from mountains. Every lightning began with a short pulse and there were no other pulses during 0.5 sec before the first one.

Close thunderstorm

Unfortunately, during the period of the observations there were no strong thunderstorms directly above the station. Two thunderstorms on July 23 and September 11, 2003 occurred in its close vicinity with lightning at the line-of-site range. An example of records in thunderstorm conditions is given in Fig. 6. Short radio pulses are seen.

The noise amplification around trigger point in thunderstorm time become stronger

During thunderstorms a number of short electromagnetic pulses were recorded. An example of the record with such a pulse is presented in the Fig. 7. Four shown data records correspond to North, Central, Central, and East electric field antennas. Central antenna signal was recorded twice by different channels of different ADC boards to diminish possible errors in the determination of signal delays between antennas for more accurate arrival angle calculations.

Some records show two short pulses in the 100 $`\mu `$s interval of the record. Besides, on July 23 there were a few records of strong electromagnetic signals of different structure possibly correspondent to radiation from electric power line.

The results of short electromagnetic pulse observations for thunderstorms on July 23 and September 11, 2003. are summarized in the Tables 1 and 2 respectively. The time and number of event, delay of the trigger relative to radio pulse, and radio pulse arrival angles are given in the Tables. The delay is positive if radio pulses were observed before trigger and negative while after it. Arrival angle is given as an inclination angle from the vertical and an azimuth from the North of the direction to the radio source.

The summary of the short pulse arrival angles is shown in the Figures 8 and 9 as inclination–azimuth plots for July 23 and September 11, 2003 thunderstorms respectively. It is clearly seen that in both cases most radio pulses came from the horizon direction. Note, that on September 11, 2003 there were several pulses which came from the direction near the vertical.

Histograms of the delays between radio pulse arrival and EAS triggers are shown in the Figures 10 and 11 for July 23 and September 11, 2003 thunderstorms respectively. Positive delays correspond to radio pulses arrived before trigger pulse. It is seen that delays are near equally distributed over the recorded interval.

The radio pulses generated during two observed thunderstorms have similar main features. Below we will analyze in details mainly the thunderstorm 11 September. The trigger array at this storm was wider: the radio pulses I and II in Fig 1 for September 11 were used and only I for July 23. The array calibration at September 11 was totally completed and the time interval between triggers was stable (Fig. 2). Only three records show two short pulses at 100 $`\mu `$s interval and a few were observed after trigger.

## 5 Discussion

I. Radio pulses

RB - EAS combined action The process which is called RB - EAS combined action take place when high energy cosmic ray particle cross the thundercloud. The thundercloud is supposed to be in RB state, what means that the maximal electric field in thundercloud is close to $`E_c`$ : $`E_m1.01.2E_c`$. Exactly these values of maximal electric field are observed in thundercloud at balloon experiments . Large scale thundercloud electric field is directed close to vertical $`z`$ .

The secondary electrons of EAS are multiplied and accelerated in thundercloud due to RB process  \- . In the same time a very large number of thermal electrons is generated due to air moleculas ionization by high energy electrons . The thermal electrons move in the air under the action of thunderstorm electric field forming an electric current pulse. The current pulse direction is determined by the direction of thundercloud electric field: it could be positive or negative depending on field direction (up or down).

The current pulse  is unipolar. Its growth time $`\mathrm{\Delta }t`$ is determined by RB process $`\mathrm{\Delta }tl_a/c`$, where $`l_a`$ is characteristic scale of relativistic electron avalanche. At thundercloud heights not far from the Station (4 – 5) km $`l_a60100`$ m ,. The pulse decay is determined by thermal electrons disappearance due to triple collisions attachment in air. The characteristic time of triple attachment at heights 4 – 6 km $`\mathrm{\Delta }t_{at}100`$ ns ,.

The form of radio pulse, generated by unipolar current pulse is bipolar . Its characteristic time scales, determined by current pulse are $`\mathrm{\Delta }t_{at}`$100 ns and $`\mathrm{\Delta }tl_a/c200÷300`$ ns. The full width of a pulse $`0.5\mu `$s Exactly these time scales determine the frequency range and the form of observed in our experiment radio pulses (Fig. 7).

Amplitude distribution of the observed radio pulses is shown at Fig. 12 and Fig. 13. In  it was shown that the amplitude of radio signal $`a`$ is proportional to the number of secondary electrons $`n_s`$ and due to this is proportional to the energy $`\epsilon _p`$ of the cosmic ray particle creating EAS . One can see from Fig. 13 that according to our observational data the integral number $`N`$ of signals in high amplitude range is falling down with amplitude as $`a^\kappa `$, $`\kappa 2`$. It is in reasonable agreement with the well known cosmic ray integral distribution ($`N\epsilon ^{1.8}`$). That coincidence could be considered as direct indication that the observed radio pulses are generated by high energy particles crossing the thundercloud. Of course statistics of radio signals shown at Fig. 13 is very poor. We present the data to demonstrate only the new opportunities for cosmic ray studies using thunderstorm radio emission. From the integral distribution presented at Fig. 13 and its comparison with our numerical simulations one can deduce that the energies of cosmic ray particles generating observed pulses lays in the interval $`2\times 10^{14}eV<\epsilon _p<2\times 10^{15}eV`$.

Electric current

The main part of radio pulses had the bipolar form (Fig. 7) and is generated by the unipolar current $`J`$ directed vertically. The pulses are coming mainly from horizontal direction (Figs. 89). For the amplitude of radio pulse in this case it follows from

$$a=\frac{2J}{cR}$$
(1)
Here $`R`$ is the distance between current pulse generated in thundercloud and trigger array. As it will be shown below $`R1`$ km. From observations of pulses amplitudes presented in Fig. 12 one can deduce that the maximal value $`a_m`$ corresponding to cosmic ray particle energy $`\epsilon _p10^{15}`$ eV is

$$a_m60\text{mV}/\text{m},\epsilon _p10^{15}\text{eV}.$$
(2)
From (1), (2) it follows

$$J_a\frac{\epsilon _{15}}{R_{km}}A$$
(3)
In this expression current is given in Amperes, $`R_{km}`$ – the distance from trigger array to RB region in kilometers and $`\epsilon _{15}`$ – the energy of CR particle which generates EAS (in $`10^{15}`$ eV). Note that at the storm 23 July the radio pulses were in quite analogous amplitude range 20 – 40 mV/m but two pulses was extremely strong $`a200`$ mV/m

Full number of thermal electrons

According to  the current $`J`$ is produced by the motion of thermal electrons ($`\epsilon 1`$ eV) under the action of RB electric field $`EE_c2`$ kV/cm:

$$J=eV_dn_T/l_p$$
(4)
Here $`e`$ – is the electron charge, $`n_T`$ – full number of thermal electrons at the characteristic length of a pulse $`l_p`$ and $`V_d`$ – electron drift velocity:

$$V_d=\frac{eE}{m\nu }$$
(5)
where $`m`$ – electron mass, $`\nu `$ – electron collision frequency in air ,. For $`\epsilon =1`$ eV and $`N_m=1.7\times 10^{19}cm^3`$ (for thundercloud height 4 – 5 km)

$$\nu =10^{12}\left(\frac{N_m}{1.7\times 10^{19}cm^3}\right)s^1$$
(6)
From (3) – (6) it follows that

$$n_T=\frac{Jl_p}{eV_d}10^{16}\frac{\epsilon _{15}}{R_{km}}$$
(7)
Dissipated energy

Now let us take into account that for ionization of one air molecule (generation one free thermal electron) the average energy $`\overline{\epsilon }=20÷30`$ eV is needed . From (7) it follows that the energy needed to generate the observed radio pulse is

$$E\overline{\epsilon }n_T10^{17}eV$$
(8)
Note that analogous energy goes to $`\gamma `$ \- quanta generation. Thus it follows from (8) that the energy dissipated in RB - EAS interaction region at least about 100 times higher than the energy of CR particle, generating EAS! That is why we call the effect ”Giant Burst”.

Of course, the dissipated energy comes from the interaction of the high energy secondary RB - EAS electrons with thunderstorm electric field. Gamma quanta, generating electrons and positrons can give support of the effect also. The detailed theoretical discussion of basic processes which determine the Giant Burst is given in the accompanying paper .

II Thunderstorm Extensive Atmospheric Showers (TEAS)

Two main problems

The trigger array (system I, II) worked stable during thunderstorm with the same counting frequency as in quiet days: one shower was registered by the joint system I-II every 2.45 sec (Fig. 2). At the same time, the cases, when radio pulses were observed, were singled out by two very significant peculiarities:

1. The inclination angles for most observed radio pulses were large: $`\theta 60^o÷90^o`$ (Figs. 89, Tables 1, 2).

We note, that the previous detailed study of inclination angles of EAS generated by the same energy CR particles ($`\epsilon 10^{14}÷10^{15}`$) at Tien Shan Scientific Station gave quite a contrary statements. For example according to  the probability of inclination angles deviation from vertical direction $`\theta =0^o`$ is proportional to $`cos^6\theta `$. It means that the pulses should be concentrated near zenith. Analogous result is obtained in : the full number $`N`$ of EAS with $`\theta <30^o`$ was twice larger than for the $`\theta >30^o`$:

$$\frac{N(\theta <30^o)}{N(\theta >30^o)}2$$
(9)
From (9) it follows that the average angle density distribution $`\rho `$ for EAS with $`\theta <30^o`$ is fourteen times higher than for $`\theta >30^o`$:

$$\frac{\rho (\theta >30^o)}{\rho (\theta <30^o)}\frac{1}{2_{30^o}^{90^o}sin\theta 𝑑\theta }_0^{30^o}\mathrm{sin}\theta d\theta 0.07$$
(10)
Thus for large angles $`\theta 60^o÷90^o`$ the probability to observe EAS under quiet conditions is very low, practically close to zero.

2. The observed delay between a trigger array signal matching the EAS coming and the radio pulse arrival to radio antenna was large. EAS was usually late and the time delay was quasi-random distributed in the 5 – 80 $`\mu \text{s}`$ timescale (Fig. 11,Table 2).

We note that time delay between a EAS and radio signals should be no more than trigger window $`5\mu \text{s}`$, as both signals are propagating with the light velocity.

Thus observed during thunderstorm specific effects are in direct contradiction with the well established usual observations of EAS in quiet conditions.

Giant burst

The qualitative explanation of these effects could be given basing on the RB-EAS combined action analysis presented in the Part I of this Section. Analyzing the radio pulse production mechanism we have seen that when EAS pass the thundercloud under RB conditions fulfilled a large amount of energy is dissipated due to the combined action of thunderstorm electric field and cosmic ray secondaries. This energy lead to the production of a giant burst of electrons and $`\gamma `$-quanta. Electrons in their final thermal state generate the observed radio pulse ,. Gamma quanta determine both peculiarities observed in Thunderstorm Extensive Atmosphere Showers (TEAS).

Qualitatively it is shown at Fig. 14. Giant number of high energy electrons and $`\gamma `$-quanta generated in RB-EAS interaction region form a cascade wich is moving in all directions. In the cascade due to Compton inelastic scattering and $`e^+e^{}`$ pair production on electrons and air nucleas the $`\gamma `$-quanta energy is rapidly diminishing. As soon as it become less than 0.5 Mev $`(\epsilon _\gamma mc^2)`$ the motion of gamma quanta changes dramatically: the energy losses are small now and the main process become quasy-elastic Compton scattering. As a consequence the motion of $`\gamma `$-quanta with the energies $`\epsilon _{ph}<\epsilon _\gamma <mc^2`$ is diffusive (below $`\epsilon _{ph}=`$20 keV photo-ionization losses dominates). The diffusion is not a straightforward motion – it lead to the time delay $`\mathrm{\Delta }t`$:

$$\mathrm{\Delta }tR^2/2D_\tau ,D_\tau =\frac{l_\tau c}{3}$$
(11)
Here R is the distance $`\gamma `$-quanta passed by diffusion, $`l_\tau `$ – mean free path of $`\gamma `$-quanta determined by Thomson scattering. In air $`l_\tau =1/(\sigma _\tau ZN_m)50m`$, where $`\sigma _\tau `$ – Thomson crossection, $`Z`$ and $`N_m`$ – average number of electrons and number density in air molecula. It follows from (10) that for $`R1km`$ the time dealy is $`\mathrm{\Delta }t100\mu s`$.

The proposed model explains both main peculiarities of TEAS:

1. According the scheme at Fig. 14 the trigger array could see the pulses mainly at large inclination angles. The direct EAS could come to trigger array rarely.
2. It is natural to suppose in accordance with Fig. 14 that the array was excited by strong gamma fluxes from giant bursts in RB-EAS regions. It explaines the random time delay of shower array signals up to $`100\mu `$s (Figs. 1011), as the $`\gamma `$-flux motion was diffusive (11).
3. Analysis shows that we registered not EAS itself but a gumma burts generated by RB-EAS process. That is why the energy of the primary cosmic ray particle could be larger than it was supposed previously.

More detailed theoretical analysis of the process is presented in the accompanying paper .

Direct observations of EAS

In two events the time delay between trigger and radio signal is small enough: ($`0.2,3.3\mu s`$). It lays inside the trigger array accuracy ($`5\mu s`$). These events shown in Fig. 15 could be considered as the direct simultaneous observation of radio pulses and EAS. We see that the pulse records are practically the same as in usual case with time delay. The inclination angles for both events are large $`\theta 87^o,89^o`$, what means practically the horizontal motion of CR particle.

Near zenith direction

From the data presented in Fig. 9, Table 2 one can see that nine pulses are singled out – they are directed close to zenith: inclination angles from $`5^o`$ to $`18^o`$. They all are peculiar pulses being wider and having an unusual structure shown in Fig. 16. The radio pulses looks generated by double unipolar current separated in time at 0.5 $`\mu `$s or about 50 – 100 m in space. Note that one of the pulses presented in Fig. 16 is a direct EAS pulse triggered at time delay 2.7 $`\mu `$s which lies inside the 5$`\mu `$s trigger window. Again no difference is seen between records of direct and delayed pulses.

We have no clear understanding of the origin of these specific pulses. The following preliminary explanation could be proposed. The current in RB - EAS interference region has two components: the current of thermal electrons and current of relativistic runaway electrons. The thermal current is stronger and generates a wide spread radio emission. The relativistic current is weaker and its radio emission is concentrated near EAS direction . It is possible that in the close to zenith case we observe the emission generated by both current pulses what can explain its bipolar form. This speculation needs of course much more detailed theoretical and experimental study.

Lightning initiation

In previous work  the lightning initiation process was studied using analogous radio technique. The first lightning radio pulse was singled out and observed for large lightning number. It was established that the radio pulse is generated by unipolar current pulse and its form, width and amplitude agree with the predictions of RB - EAS interaction theory . Thus it was supposed that the lightning is initiated by combined action of RB - EAS generated by primary cosmic ray particles with energy $`\epsilon _p10^{16}`$ eV. In the present work the radio pulses and EAS are observed directly. The radio pulses have quite analogous form and frequency range as in . The difference in amplitudes is explained by smaller energy of cosmic ray primary particle ($`\epsilon 10^{15}`$ eV). Note that the observations of lightning first pulses from the faraway thunderstorm at Tien Shan gave the results quite analogous obtained in .

Thus one can state that the observations presented in the paper confirm the lightning RB - EAS initiation mechanism proposed in , .

Note that the radio pulses associated with extensive atmospheric showers were observed in a number of previous works  \- . The width and amplitude of observed pulses were quite different – from $`5\mu `$s  up to 100 $`\mu `$s , depending on conditions and technique used. In ,  the connection between the existence of radio pulses and weather conditions (clouds, thunderstorms) was indicated.

## 6 Conclusions

In conclusion we formulate briefly main results of the experiments described in the paper:

1. At Tien Shan Mountain Scientific Station a new EAS - Radio installation is constructed. The installation consists of a wide spread EAS trigger array and a high time resolution radio interferometer what allows to study the EAS and radio pulses emission simultaneously.
2. The absence of simultaneous EAS and radio pulses in a quiet (non-thunderstorm) atmosphere is established.
3. During two thunderstorms 150 simultaneous EAS and radio pulses were observed. The distribution of radio pulses on amplitude, form, inclination and azimuth direction was studied. The following main peculiarities are established:

1. The radio pulses came mainly from the horizontal direction.
    2. There is a significant random time delay (up to 100 $`\mu `$s) between trigger signal and radio pulse.
4. The analysis of the obtained experimental data allows to establish the existence of a new phenomena – giant electron-gamma bursts generated by high-energy cosmic ray particles in thunderclouds. The burst energy could be $`10^2÷10^3`$ times larger than the particle energy due to the RB-EAS interaction.
5. The burst generated by high energy cosmic ray particles $`\epsilon >10^{16}`$eV can serve for lightning initiation. That confirms the lightning initiation mechanism studied in .

The current generated in thunderclouds due to RB-EAS interaction is proportional to the cosmic ray particle energy $`\epsilon _p`$. It means that the radio pulse amplitude $`a`$ for $`\epsilon _p10^{18}÷10^{19}`$ eV should be extremely high. The pulse has definite characteristics quite analogous to the established in the present work. It could be observed from the immense distances – 1000 km and even more.

The further simultaneous observations of EAS together with giant bursts of X rays and radio pulses can serve for the development of a new effective method for the radio detection of high energy cosmic ray particles.

Acknowledgements

The authors are grateful to Prof V.L.Ginzburg, Prof. E.L.Feinberg and Dr. H.Carlson foe useful discussions. The work was supported by EOARD- ISTC grant #2236, ISTC grant #1480, by the President of Russian Federation Grant for Leading Scientific Schools Support and by the Russian Academy Fundamental Research Program ”Atmosphere Physics: Electric Processes, Radio Physics Methods”.


format: text

text:
# Kerr black holes as retro-MACHOs

## 1 Introduction

Gravitational lensing of electromagnetic waves is a well-understood phenomenon (see e.g. Schneider et al. sef (1992) for a comprehensive treatment) predicted by the theory of General Relativity. In the last years, gravitational lensing has been also used as an observational tool to investigate the mass distribution both in galaxy clusters (using distant QSOs as sources, Canizares canizares1982 (1982)) and in galactic halos. In the latter case, the phenomenon, which is called gravitational microlensing, led to the discovery of Massive Astrophysical Compact Halo Objects (MACHOs) in the halo of our Galaxy (Alcock et al. alcock1993 (1993), Aubourg et al. aubourg1993 (1993)).

As a dark object moves across the source-observer line of sight it acts as a gravitational lens leading to the formation of two unresolved images of the source and to an overall source light magnification. Despite the low probability of observing such events, Paczyńsky (pacz (1986)) showed that the continuous observation of at least $`10^6`$ stars implies a non-negligible chance of detecting a microlensing event towards the LMC.

Since the early 1990s several collaborations (MACHO (see e.g. Alcock et al. alcock1993 (1993)), EROS (Aubourg et al. aubourg1993 (1993)) and OGLE (Stanek et al. ogle (1995))) have been monitoring millions of stars towards the LMC, the SMC, as well as towards the galactic center detecting hundreds of microlensing events due to MACHOs along the line of sight. Unfortunately, the physical nature of MACHOs is still unknown, since the observed light curves of a microlensing event can be reproduced by models depending on a number of free parameters (mass, distance and transverse velocity of the lens). So, different populations of objects such as brown dwarfs, white dwarfs, main sequence stars and black holes may be consistent with microlensing observations. In particular, it is generally accepted that a population of black holes may exist in the Galaxy and it has been claimed that these compact objects have already been detected in microlensing surveys (Quinn et al. quinn (1999), Bennett et al. bennett (2002), Agol et al. agol (2002)). At least 6 extremely long events detected by the MACHO, GMAN and MPS collaborations towards the galactic bulge (Bennett et al. bennett (2002)) exhibit very strong microlensing parallax signals which lead to mass estimates up to $`10M_{}`$. Since the estimated upper limits on the absolute lens brightness are always less than $`1L_{}`$, this favors that the stellar black hole hypothesis (Mao et al. mao (2002)) and suggests that a substantial fraction of the galactic lenses may be massive stellar remnants.

Despite some experimental uncertainties, the theory of gravitational lensing, which was originally developed in the weak field approximation, has successfully explained all gravitational lensing observations from the giant blue luminous arcs in the rich cluster of galaxies (Lynds & Petrosian lp (1989), Soucail et al. soucail (1987)) to the first Einstein ring MG1131+0456 at redshift $`z1.13`$ (Hewitt et al. hewitt (1988)) to microlensing events detected towards the SMC, LMC, the galactic bulge and the M31 galaxy (see e.g. Alcock et al. alcock2000 (2000)).

Moreover, by using the next generation of high resolution imaging telescopes, it would be possible to test the theory of gravitational lensing in the presence of strong gravitational fields since the inspection of the images formed will allow us to investigate regions very close to the surface of compact massive objects such as neutron stars and black holes.

In a very interesting paper Holz & Wheeler (hw (2002)) proposed that a Schwarzschild black hole may act as retro-lens or retro-MACHO. This means that if the black hole is illuminated by a powerful light source, photons, due to its strong gravitational field, can be deflected by large bending angles. In the case of perfect alignment of source, black hole and Earth (with the Earth in the middle), the bending angles are odd multiples of $`\pi `$ leading to a series of circular rings.

However, black holes would be characterized by a non-zero intrinsic angular momentum, which breaks the Schwarzschild spherical symmetry and affects the gravitational field around the compact object. Thus, we expect a modification in the phenomenology of the gravitational lensing. Here, we investigate the retro-lensing phenomenon for slowly rotating Kerr black holes considering photon trajectories close to the black hole equatorial plane. Our treatment follows the second order expansion (given by Bozza bozza2 (2003)) of the large deflection angles of light rays winding around the retro-MACHO in the impact parameter $`b`$.

This paper is structured as follows. In Sect. 2 we briefly review the Schwarzschild retro-MACHO lensing results obtained by Holz & Wheeler (hw (2002)). In Sect. 3 we extend the Holz & Wheeler (hw (2002)) results to Kerr black holes . In Sect. 4 we investigate different geometries among the lens, the source and the observer and discuss the possibility that the next generation high resolution telescopes may have the capabilities necessary to detect such lensing events.

## 2 Schwarzschild retro-MACHOs

The theory of General Relativity predicts that the gravitational field of a massive object deflects light rays of a background source, if the impact parameter is sufficiently close to the massive object. This phenomenon is referred to as gravitational lensing and the massive object causing a detectable deflection is called a gravitational lens.

In gravitational microlensing, photons from a distant source suffer a very small angular deflection and the weak field approximation always works satisfactorily. However, strong field effects may be extremely important if the lens is a collapsed compact object. In fact, very close to a neutron star or a black hole, GR plays an important role and the study of the relativistic images which form (as a consequence of light deflection) allows us to investigate the regions close to the event horizon. The light bending angle in the presence of a static black hole is not limited to small angles but may reach large values ($`\pi `$ and odd multiples of $`\pi `$ if source, observer and retro-MACHO are perfectly aligned, see Holz & Wheeler hw (2002)) so that photons emerge in the direction of the source itself. Therefore, a Schwarzschild black hole, if illuminated by a powerful light source, shines back with a series of micro-arcsecond rings.

In the case of perfect alignment between the source, the observer and the lens the image of the source, in the plane orthogonal to the line of sight and containing the source, is an annulus with outer and inner radii corresponding to the appropriate impact parameters for photons coming from the top and bottom parts of the source.

However, in the general case the source, the observer and the black hole are not aligned. We consider a reference frame centered at the source $`S`$ and with axes $`X`$ and $`Y`$ on the Earth orbital plane (see Fig. 1). We choose $`Y`$ to contain the lens in the $`ZY`$ plane and define the misalignment angle $`\beta `$ as the angle between the observer-lens direction ($`\stackrel{}{OL}`$) and the $`Y`$ axis. Note that the case with $`\beta =0`$ and Earth in the orbit position $`O`$ corresponds to perfect alignment.

As the Earth moves on its orbit, the image rings deform into arcs whose angular extent is $`\mathrm{\Delta }\mathrm{\Theta }2\mathrm{tan}^1(R_S/D_{OS}\mathrm{sin}\theta ^{})`$, where $`\theta ^{}`$ is the angle between the source-observer and observer-lens ($`\stackrel{}{SO^{}}`$ and $`\stackrel{}{O^{}L}`$) directions, respectively.

The total area of each image is thus given by (Holz & Wheeler hw (2002))

$$A_I\pi (b_o^2b_i^2)\frac{\mathrm{\Delta }\mathrm{\Theta }}{2\pi }=(b_o^2b_i^2)\mathrm{tan}^1\left(\frac{R_S}{D_{OS}\mathrm{sin}\theta ^{}}\right),$$
(1)
where $`b_o`$ and $`b_i`$ represent the outer and inner impact parameters corresponding to bending angles

$$\alpha _{d,B}=(\pi \delta )\alpha ,\alpha _{d,S}=(\pi +\delta )\alpha $$
(2)
for the bigger and smaller image, respectively. Here $`\delta `$ is the angle between the $`\stackrel{}{SL}`$ and $`\stackrel{}{O^{}L}`$ directions 111Of course, since the observer moves around the Sun on a circular orbit of radius $`D_{O^{}S}1AU`$ with angular velocity $`\omega _E2\pi /T_E`$ (being $`T_E`$ the orbital period), the Earth position angle $`\varphi `$ and the angles $`\delta `$ and $`\theta ^{}`$ (see Fig. 1) are time dependent. (i.e. the deflection angle needed for a source photon to reach the observer).

The amplification of each image is given by the ratio between the area of the image given by eq. (1) to that of the source $`A_S=R_S^2/(4D_{OS}^2)`$, i.e.

$$\mu _{B,S}(b_o^2b_i^2)\mathrm{tan}^1\left(\frac{R_S}{D_{OS}\mathrm{sin}\theta ^{}}\right)\frac{D_{OS}^2}{\pi D_{OL}^2R_S^2}.$$
(3)
Finally, setting $`b_o`$ and $`b_i`$ to the appropriate impact parameters, the total image amplification is

$$\mu =\mu _B+\mu _S.$$
(4)
Obviously, in the evaluation of the total amplification $`\mu `$ we need the knowledge of the angles $`\theta ^{}`$ and $`\delta `$ defined above which, from inspection of Fig. 1, are given by

$`\theta ^{}(t)`$ $`=`$ $`\mathrm{arccos}\left({\displaystyle \frac{\stackrel{}{SO^{}}\stackrel{}{O^{}L}}{D_{O^{}S}D_{O^{}L}}}\right),`$
$`\delta (t)`$ $`=`$ $`\mathrm{arcsin}\left({\displaystyle \frac{D_{O^{}S}}{D_{O^{}L}}}\mathrm{sin}(\beta \gamma )\right),`$ (5)
where $`\gamma =\mathrm{arcsin}(D_{OS}\mathrm{sin}\beta /D_{SL})`$.

At this stage, we only need to evaluate the photon impact parameter values $`b_o`$ and $`b_i`$ so that light rays coming from the source are deflected toward the observer.

For a Schwarzschild black hole with mass $`M`$, Chandrasekhar (chandra (1983)) has found that, in the case of large deflection angle $`\alpha _d`$, the photon impact parameter, in geometrical units $`G=c=1`$, is given by

$$b(\alpha _d)=b_{1,S}+b_{2,S}e^{\alpha _d},$$
(6)
where

$$\begin{array}{c}b_{1,S}=3\sqrt{3}M\hfill \\ b_{2,S}=648\sqrt{3}e^\pi M(\sqrt{3}1)^2/(\sqrt{3}+1)^2.\hfill \end{array}$$
(7)
Consequently, for each position of the Earth around the Sun it is possible to evaluate the angles $`\theta ^{}(t)`$ and $`\delta (t)`$, the bending angles $`\alpha _{d,B}`$ and $`\alpha _{d,S}`$ and, through eqs. (3) and (6), the total image amplification (4) as seen by the observer.

The expected light curves (i.e. magnitudes $`m`$ corresponding to a certain amplification $`\mu `$ as a function of time) for a retro-lensing event due to a Schwarzschild black hole are shown in Fig. 2 for different misalignment angles $`\beta `$.

## 3 Kerr retro-MACHOs

The results by Holz & Wheeler (hw (2002)) on retro-lensing events hold in the case of Schwarzschild black holes. However, in general, a black hole is characterized by a non-zero intrinsic angular momentum which breaks the spherical symmetry and affects the gravitational field around the compact object.

Indeed, black holes are the ultimate stage of stellar evolution and may form by SNII explosions. Since stars rotate, angular momentum conservation implies that black holes also rotate. Black hole rotation is also confirmed by recent $`X`$-ray observations of Fe line width (Zakharov et al. zaklinee (2002)). Moreover, super massive black holes at the center of QSOs, AGNs and galaxies show beamed jet emission implying that they have non zero angular momentum.

Thus, due to black hole angular momentum $`J`$, we expect a modification in the phenomenology of both the usual gravitational lensing and the retro-MACHO microlensing. Such modifications are evaluated by a second order expansion of the photon deflection angles in the spin parameter $`a=J/(McR_{\mathrm{Sch}})`$. Here, $`R_{\mathrm{Sch}}=2GM/c^2`$ is the black hole Schwarzschild radius. 222 Note that the definition we use for the rotation parameter $`a`$ corresponds to half of the value in the usual notation $`a=2J/(McR_{\mathrm{Sch}})`$. Indeed, following Bozza (bozza2 (2003)) we measure distances in units of the Schwarzschild radius and not in units of the gravitational radius as is usually done (see e.g. Shapiro & Teukolsky st (1983)). We consider light trajectories on the black hole equatorial plane and a more general treatment will be addressed elsewhere (Zakharov et al. trajectories (2003)). However, we expect that the maximal effect due to the black hole spin occurs for photons lying in the equatorial plane while for partially inclined orientations the effect becomes closer to the Schwarzschild case (that is recovered for $`\theta =0`$).

A simple and reliable method to investigate the subject has been recently proposed by Bozza (bozza2 (2003)) (for an exhaustive review of the topics see also Bozza 2002a ) who revisited the Schwarzschild and Kerr black hole lensing in a strong gravitational field.

Considering light rays on the black hole equatorial plane ($`\theta =\pi /2`$), the Kerr line element in Boyer-Lindquist coordinates ($`t`$, $`x`$, $`\theta `$, $`\varphi `$) is given by

$$ds^2=A(x)dt^2B(x)dx^2C(x)d\varphi ^2+D(x)dtd\varphi ,$$
(8)
with

$$\begin{array}{ccc}A(x)=1\frac{1}{x}\hfill & & \\ B(x)=\left(1\frac{1}{x}+\frac{a^2}{x^2}\right)^1\hfill & & \\ C(x)=x^2+a^2+\frac{a^2}{x}\hfill & & \\ D(x)=2\frac{a}{x},\hfill & & \end{array}$$
(9)
where $`x`$ is the distance from the spinning black hole measured in units of the Schwarzschild radius in natural units $`R_{\mathrm{Sch}}=2M`$ and $`a=J/(2M^2)`$.

Following Bozza (bozza2 (2003)) the impact parameter $`b`$ of light rays approaching the black hole is

$$b=\frac{D_0+\sqrt{4A_0C_0+D_0^2}}{2A_0}.$$
(10)
Here and in the following all the metric functions with the subscript $`0`$ are evaluated at $`x=x_0`$ which is the photon minimum distance from the black hole.

After a straightforward calculation, the whole deflection angle (i.e. the bending angle for a photon with impact parameter $`b`$) is given by

$$\alpha _d=2_{x_0}^{\mathrm{}}\frac{d\mathrm{\Phi }}{dx}𝑑x\pi ,$$
(11)
where

$$\frac{d\mathrm{\Phi }}{dx}=P_1(x,x_0)P_2(x,x_0),$$
(12)
with

$$\begin{array}{c}P_1(x,x_0)=\frac{\sqrt{B}(2A_0Ab+A_0D)}{\sqrt{CA_0}\sqrt{4AC+D^2}}\hfill \\ P_2(x,x_0)=\left(A_0A\frac{C_0}{C}+\frac{b}{C}(AD_0A_0D)\right)^{1/2}.\hfill \end{array}$$
(13)
The deflection angle $`\alpha _d`$ grows as $`x_0`$ decreases diverging at the radius of the photosphere $`x_m`$ (to be defined below). In this case the photon is captured by the black hole. Defining the new variable $`z=(xx_0)/x`$ and hence $`x=x_0/(1z)`$, the integral in eq. (11) becomes

$$\alpha _d=_0^1R(z,x_0)f(z,x_0)𝑑z\pi ,$$
(14)
with

$`R(z,x_0)`$ $`=`$ $`2{\displaystyle \frac{(1y_0)}{A^{}(x)}}P_1(x,x_0)`$
$`f(z,x_0)`$ $`=`$ $`P_2(x,x_0).`$ (15)
The function $`R(z,x_0)`$ is regular for any values of $`z`$ while the function $`f(z,x_0)`$ diverges as $`z0`$. To find the divergence order, the argument of the square root in $`f(z,x_0)`$ may be expanded in powers of $`z`$ obtaining (Bozza bozza2 (2003))

$$f_0(z,x_0)=\frac{1}{\sqrt{c_1(x_0,a)z+c_2(x_0,a)z^2}},$$
(16)
where $`c_1(x_0,a)`$ and $`c_2(x_0,a)`$ represent the first and second order coefficients of the expansion. Therefore, when $`c_1(x_0,a)`$ is different from zero the leading order term in the integrand is $`z^{1/2}`$, so that the integral in eq. (14) converges. Instead, the integral diverges in the case $`c_1(x_0,a)=0`$, corresponding to photons being captured by the black hole. Accordingly, the photosphere radius $`x_m`$ is defined as the outermost solution of equation $`c_1(x_0,a)=0`$. This condition is equivalent to $`8a^2x_0(32x_0)^2=0`$ (for details see Bozza bozza2 (2003)). Of the three solutions of the cubic equation above, one is irrelevant as it lies inside the horizon while the other two correspond to the co-rotating and counter-rotating photons. The counter-rotating case may also be obtained by taking $`a<0`$. The roots are displayed in Fig. 3.

Once the photosphere radius $`x_m`$ (for a given value of $`a`$) is known, one can evaluate the deflection angle $`\alpha _d`$ from eq. (14) . Thus, retaining only terms up to $`O(x_0x_m)`$, one obtains

$$\alpha _d(x_0)=k_1\mathrm{ln}\left(\frac{x_0}{x_m}1\right)+k_2+O(x_0x_m),$$
(17)
where

$$\begin{array}{cc}k_1=\frac{R(0,x_m)}{\sqrt{c_2(x_m,a)}}\hfill & \\ k_2=\pi +k_1\mathrm{ln}\frac{2(1y_m)}{A^{}(x_m)x_m}+I_R\hfill & \\ I_R=_0^1[R(z,x_m)f(z,x_m)R(0,x_m)f_0(z,x_m)]𝑑z.\hfill & \end{array}$$
(18)
The impact parameter $`b`$ in eq. (10) can now be expanded in powers of $`(x_0x_m)`$ up to the second order term (the first order coefficient vanishes) obtaining

$$bb_0(x_m)+b_2(x_m)(x_0x_m)^2.$$
(19)
Eliminating the closest approach distance $`x_0`$ in eqs. (17) and (19) and rearranging the terms one finally gets (Bozza bozza2 (2003))

$$\alpha _d(b)=\overline{k}_1\mathrm{ln}\left(\frac{b}{b_0(x_m)}1\right)+\overline{k}_2+O(bb_0(x_m)),$$
(20)
with

$$\begin{array}{c}\overline{k}_1=\frac{k_1}{2},\hfill \\ \overline{k}_2=k_2+\frac{k_1}{2}\mathrm{ln}\frac{b_2(x_m)x_m^2}{b_0(x_m)}.\hfill \end{array}$$
(21)
In Fig. 4 we plot the coefficients $`\overline{k}_1`$ and $`\overline{k}_2`$ as a function of the black hole spin parameter. We note that these coefficients diverge for $`a0.5`$ implying that eq. (20) can be used only for low enough values of $`a`$ for which $`x_0x_m1`$ holds.

Therefore, by solving eq. (20) with respect to the impact parameter $`b`$, one finally obtains the generalization of the Chandrasekhar relation for Schwarzschild black holes given in eq. (6) to Kerr black holes

$$b(\alpha _d)b_{1,K}+b_{2,K}e^{\alpha _d/\overline{k}_1},$$
(22)
with

$$b_{1,K}=2b_0(x_m)M,b_{2,K}=2b_0(x_m)Me^{\overline{k}_2/\overline{k}_1}.$$
(23)
Note that the factor $`2`$ in the previous relations (23) has been introduced to make our expression comparable with that in eq. (6). It is also easy to verify that for $`a0`$, eq. (22) reduces to the result (6) given by Chandrasekhar (chandra (1983)) for Schwarzschild black holes.

We are now ready to apply this formalism to the retro-lensing phenomenon. In Fig. 5, assuming $`\alpha _d=\pi `$, we plot the impact parameter $`b(\pi )`$ as a function of the black hole spin parameter $`a`$. As one can see, for co-rotating photons ($`a>0`$) the impact parameters decreases with respect to the Schwarzschild case (i.e. photons may get closer to the black hole).

Let us now consider the Sun as the light source and a black hole of mass $`10M_{}`$ at a distance of $`0.01`$ pc. Since the solar luminosity in the optical band is $`L_{}4\times 10^{33}`$ erg $`s^1`$, corresponding to an optical visual magnitude $`m_{}=26.8`$, the expected magnitude $`m`$ of the retro-lensed image can be evaluated as

$$m=m_{}2.5\mathrm{log}\mu ,$$
(24)
$`\mu `$ being the total amplification given by eq. (4). In the case of a Schwarzschild black hole with $`a=0`$, the expected apparent magnitude is given in Fig. 2 as a function of the time $`t`$ (corresponding to the Earth position on its orbit), for different values of the misalignment angle $`\beta `$ (see also Holz & Wheleer hw (2002)).

In the case of a Kerr black hole we can evaluate the Sun image amplification by using eqs. (3) -(5), in the most favorable case of perfect alignment ($`\beta =0`$). For this purpose we only need to evaluate, through eq. (22), the impact parameters corresponding to the bending angles for which light rays emitted by the Sun may reach the observer after deflection of the angles $`\alpha _{d,B}`$ and $`\alpha _{d,S}`$ defined in Section 2. For a $`10M_{}`$ Kerr black hole with spin parameter $`a=0.1`$ at a distance of $`0.01`$ pc and misalignment angle $`\beta =0`$, the expected retro-lensing light curve is given in Fig. 6 (red solid line). In the same Figure we also give the light curve (dashed red line) corresponding to counter-rotating photons. As expected, for a Kerr black hole with spin parameter $`a`$, only a modification up to the second order in the retro-MACHO amplification appears, as is evident by inspecting Fig. 6 where the lightcurves for the Kerr black hole case are reported (dashed and solid red lines) together with the corresponding one for a Schwarzschild black hole with the same mass and $`\beta =0`$ (central black solid line). 333Co-rotating and counter-rotating photons should produce an image scintillation due to their relative time delay. However, the time scale for such scintillations goes from $`10^5`$ s to $`10`$ s for black holes with mass in the range $`1`$ M-$`10^6`$ M. Due to the short scintillation time scale with respect to the integration time required by observation, the net result is simply that of seeing the average, which coincides with the light curve for the Schwarzschild black hole. We also notice that we have neglected multiple rotations of the photons around the black hole. As emphasized by Holz & Wheeler (hw (2002)), photons may form a series of deformed rings corresponding to deflection angles that differ by $`2n\pi `$ (with $`n`$ as integer number). However, it is clear that the brightest image is that corresponding to the outermost ring (obtained for $`\alpha _d\pi `$) so that neglecting the other images does not change substantially the retro-lensing image brightness.

## 4 Discussion

It has been shown by Holz & Wheeler (hw (2002)) that the limiting distance $`D_L`$ at which the rings around a black hole of mass $`M`$ illuminated by our Sun can be observed is given by

$$D_L=0.02\mathrm{pc}\times \left[10^{(m30)/2.5}(M/10M_{})^2\right]^{1/3},$$
(25)
where the ring magnitude of the retro-lensing event has been chosen to be $`30`$ (see Fig. 2). For example, if the Sun is the light source and the black hole has mass $`10`$ M, the limiting distance at which the retro MACHO can be detected is of only $`0.02`$ pc for an instrument with $`m=30`$ (such as HST). On the other hand, it is very unlikely that black holes so close to the Earth exist. However, the most convenient geometry to get luminous rings around a black hole would be that of a very bright star close to a black hole at some distance $`D_L`$ from Earth. 444It has been shown that the best chance of observing retro-lensing images in the future is by looking towards the galactic center black hole in Sgr A around which a very bright star (named S2) with a mass of about $`15M_{}`$ is orbiting. The resulting magnitude of the retro-lensing images in the Schwarzschild case are in the range $`33.336.8`$ (depending on the star distance from the black hole) in the K-band, close to the limiting magnitude of the next generation of space-based telescopes (De Paolis et al. 2003a .) In this case, and for the most general geometrical configuration of source, lens and observer with relative distances $`D_{OL}`$, $`D_{LS}`$ and $`D_{OS}`$, the instrumental limiting magnitude $`\overline{m}`$ necessary to see the rings is thus given by

$$\overline{m}m_\lambda 2.5\mathrm{log}\mu ,$$
(26)
$`m_\lambda `$ being the magnitude of the source in a certain electromagnetic band and $`\mu `$ the ring amplification given by eq. (4). In the particular case of perfect alignment between observer, lens and source, the previous relation reads

$$\overline{m}m_\lambda 2.5\mathrm{log}\left[1.04\left(\frac{M}{M_{}}\right)^2\left(\frac{R_{}}{R_S}\right)\frac{D_{OS}^2}{D_{LS}D_{OL}^2}\right],$$
(27)
which reduces to eq. (25) taking the geometry adopted by Holz & Wheeler (hw (2002)). Notice that in the previous equation we are considering all the distances ($`D_{OS},D_{LS}`$ and $`D_{OL}`$) measured in cm.

As is clear from Fig. 6, the spin effect on the image brightness is a minor effect and so it cannot increase substantially the image amplification with respect to the Schwarzschild black hole case. Moreover, as already stated, our analysis is confined to photons lying in the black hole equatorial plane and the retro-lensing image magnitude is calculated by using eq. (4). This equation is rigorously valid only in the Schwarzschild case since it assumes that each part of the image is isoluminous (and this is not true for Kerr black holes). However, if we limit the estimate to only the total magnitude of the retro-lensing images, we do not expect a substantial change with respect to the Schwarzschild case since the increase in the luminosity of the corotating edge should approximately compensate the luminosity decrease of the counter-rotating edge.

We emphasize again that our treatment is valid only for slowly rotating black holes but it is clear that increasing the rotation parameter $`a`$, the brightness of the retro-lensing images also increases. A complete (purely numerical) treatment of this issue - valid for arbitrarily fast rotating black holes - will be published elsewhere (De Paolis et al. 2003b ).

Since in the slow rotation limit the light curves are almost indistinguishable from the Schwarzschild case, the most important effect of the black hole spin is that of deforming the ring shape with respect to the Schwarzschild case (Peter peter76 (1976), Bray bray (1986)). A detailed analysis of this effect will be presented elsewhere (Zakharov et al. trajectories (2003)). Detecting the shape of the ring images might give a unique possibility of investigating the black hole parameters (spin $`a`$, mass $`M`$ and relative position of observer-lens-black hole) in an active way.

To really observe the ring shape, the detecting instrument must have a high enough angular resolution in the range $`110`$ $`\mu `$arcsec for most cases of interest. Future space-based instruments like the Space Interferometry Mission (SIM) (see Shao sim (1998) for technical details), GAIA Satellite (see e.g. Belokurov & Evans gaia (2002)), FAME (Johnston et al. fame (2000)) and MAXIM (maxim (2003)) would have the angular resolution capabilities to resolve such images and therefore will be able to observe retro-MACHOs if the instrumental limiting magnitude satisfies eq. (26).


text:
# Spiral shock detection on eclipse maps: Simulations and ObservationsBased on observations made with the Isaac Newton and William Herschel telescopes operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias.

## 1 Introduction : spiral shocks in cataclysmic variables

Cataclysmic variables (CVs) are binaries where an accretion disk is formed around the white dwarf from gas escaping a companion red dwarf (Warner 1995). Their study is attractive for understanding binary formation and evolution and the physical processes occuring in the accretion flow. For example, the semi-regular radiation outbursts from the accretion disks around white dwarfs in the CV sub-class of the dwarf novae are thought to be driven by a thermal instability within the disk (Hameury 2003 and references therein). The optical radiation pattern on the accretion disk is asymmetric and is dominated either by the “bright spot” caused by the impact of the gas stream on the accretion disk in quiescence (Marsh 1989) or the “spiral structure” raised by tides from the companion star on the outer accretion disk during outburst (Steeghs 2001; Boffin 2001).

IP Pegasi is an eclipsing dwarf nova and one of the best studied CVs thus providing a paradigm for the accretion disk model. Spiral shocks in its accretion disk have been found throughout the various stages of an outburst (rise, maximum, decline; Steeghs et al. 1997; Harlaftis et al. 1999; and Morales-Rueda et al. 2000, respectively), using the Doppler tomography technique (Marsh & Horne 1988). The phase-resolved spectra in Harlaftis et al. (1999), which were used to reconstruct the Doppler maps - two-dimensional velocity maps - of the spiral shocks, were also used to make eclipse light curves in various emission lines. The light curves (integrated intensity versus phase) were then analysed using the eclipse mapping technique (Baptista 2001) in order to reconstruct the surface brightness distribution of the accretion disk and thus define the spatial extent of the spiral arms (Baptista et al. 2000). Direct comparison of the Doppler (velocity) and eclipse (spatial) maps showed that the outer disk - largely affected by the velocity field of the spiral shocks - in IP Pegasi is largely non-Keplerian and that the spiral arms extend between 0.2-0.6 $`R_{L_1}`$, contributing between 16-30 % of the total line flux.

Here, we address the parameters affecting the reconstruction of the spiral structure using the eclipse mapping technique. We undertake simulations in order to demonstrate the effect on the detectability of spiral arms in situations with an unfavorable combination of low S/N, low time resolution, and the dilution of the shock light by the disk. Hence, we define the observational requirements for future observations of spiral shocks using the eclipse mapping technique. Finally, in the light of the simulations performed we interpret the eclipse maps built from light curves of IP Pegasi during outburst (Morales-Rueda et al. 2000). This is an epoch where the spiral shocks have been excited in the accretion disk, as shown by Doppler tomography.

## 2 Simulations

Doppler tomography of IP Peg has clearly revealed the disk spiral structure on a number of occassions and showed its strength in imaging sub-structures within the accretion disk (Steeghs 2001). However, the disk spiral structure of IP Peg is, apparently, not imaged as clearly using the eclipse mapping technique (Baptista et al. 2000), even though the eclipse light curves are from the same outburst maximum (November 1996). There, the same phase-resolved spectra show one of the best examples of spiral shocks in an accretion disk using Doppler tomography as the reconstruction technique (Harlaftis et al. 1999). As a consequence, we have performed a series of simulations in order to better understand the structures reconstructed with eclipse mapping. The reconstructions were obtained using the PRIDA code (Baptista & Steiner 1993), which computes an eclipse map as an array of 51x51 pixels, with a size of 2 $`R_{L_1}`$ and centred on the white dwarf. Any additional uneclipsed component can also be modelled as a constant flux. We adopt the binary geometry of Wood & Crawford (1986) with ($`i,q`$)=($`81\mathrm{°}`$, 0.5), where $`i`$ is the inclination and $`q`$ is the mass ratio (see also section 3.1). For the reconstructions we adopted the default of limited azimuthal smearing (Rutten, van Paradijs & Tinbergen 1992), which is better suited for recovering asymmetric structures than the original default of full azimuthal smearing (e.g., Baptista 2001). For more details on the technique, see Baptista (2001). Using a model map with a symmetric disk and a spiral structure, we test the limits of the technique in resolving the spiral arms. The parameters we change are either defined by the observing technique, geometrical constraints or physical processes :

* spiral structure orientation
* phase resolution
* signal-to-noise ratio
* accurate zero phase of mid-eclipse
* spiral structure dilution by the symmetric disk component
* brightness difference between the two spiral arms

For the simulations, we use as input a model map consisting of two spirals superimposed on a symmetric disk. The spirals have a radial width of 0.08 $`R_{L_1}`$ and a radial range between 0.2–0.6 $`R_{L_1}`$, as found in Baptista et al. (2000). The trace of the spiral is defined by,

$`R_{\mathrm{spiral}}(\theta )=(\theta \theta _0)/\alpha `$

where $`\theta _0`$ is the angle between the tangent to the spiral pattern (extrapolated to the disk centre at $`R_{\mathrm{spiral}}=0R_{L_1}`$) and the line-of-centres and increases counter-clockwise in steps of $`30\mathrm{°}`$, using as a reference the observed rotation angle at $`46\mathrm{°}`$ (Baptista et al. 2000); $`\theta `$ is the angle between the tangent to the spiral pattern at $`0<R_{\mathrm{spiral}}<R_{L_1}`$ and $`\alpha (=180\mathrm{°})`$ is the opening angle of the spiral, the change in angle ($`\theta \theta _0`$) from $`R_{\mathrm{spiral}}=0`$ to $`R_{\mathrm{spiral}}=0.6R_{L_1}`$. Note that the above definition is for linear spiral arms (see Binney & Tremaine 1988 for logarithmic spiral arms). Fig. 1 shows graphically the spiral definition with $`\theta _0`$ and $`\alpha `$ as defined above. Setting a small (large) $`\alpha `$ results in a tightly wound (open) spiral arm. Our model map in Fig. 2 has an underlying symmetric disk containing 40 % of the total flux and a flat brightness distribution which drops sharply at R=0.6 $`R_{L_1}`$ following a $`exp[|(rr_0)/dr|^3]`$ law. The spirals are also cut off with the same radial dependency, both at the inner and outer radius.

### 2.1 Spiral structure orientation ($`\theta _0`$)

In Fig. 2, we present model maps of spiral arms at different rotation angles which we use as input images for the simulations, so that we compare the predicted eclipse light curves with the observed ones. The dashed lines in the left hand panels show the contribution of the two spiral arms to the light curves; the dash-dotted lines show the contribution of the symmetric disk component to the light curves; and the solid lines show the total light curve. The rotation angle $`\theta _0`$ of the spiral arm main axis is indicated in each panel. The model surface brightness distribution maps are shown in a logarithmic grayscale (dark regions are brighter). There, the dotted lines depict the primary Roche lobe, the gas stream trajectory, and the disk radii at 0.2 and 0.6 $`R_{L_1}`$. A cross marks the disk centre. The stars rotate counter-clockwise and the secondary is to the right of each panel. The two spiral arms, of equal brightness, are labeled ”A” and ”B” in two of the panels. The change in the shape of the eclipse is notable as the spiral arms rotate from 16$`\mathrm{°}`$ to 166$`\mathrm{°}`$. The wider eclipse ingress shifts to a wider egress of the spiral arms. The spiral arm light curves show best the change in eclipse shape with rotation angle. Reconstructions obtained from these light curves confirm the results of previous simulations (Baptista et al. 2000), namely, that the orientation and radial position of the arms are well recovered in the eclipse maps. The panels present spiral models for increasing steps of $`30\mathrm{°}`$ starting at $`16\mathrm{°}`$. The eclipse shape of the spiral model with rotation angle $`46\mathrm{°}`$ resembles best the observed eclipse maps presented here and in Baptista et al. (2000). This is clearer in the next figure which presents the simulation results with the phase resolution of the observed data. The eclipse shape of the spiral model with $`\theta _0=136\mathrm{°}`$ resembles the eclipse light curves of Baptista et al. (2002).

### 2.2 Phase resolution ($`\mathrm{\Delta }\varphi `$) and S/N ratio

Fig. 3 shows the effect of reducing the phase resolution of the light curve on the map reconstructions. These and all following simulations were performed with the adopted spiral model ($`\theta _0=46\mathrm{°}`$). The left hand panels show data light curves of S/N=50 for three phase resolutions: $`\mathrm{\Delta }\varphi `$= 0.005, 0.010 and 0.015 cycles. Gaussian noise has been added in the light curves. Note that the phase resolution of $`\mathrm{\Delta }\varphi =0.010`$ is similar to that of the average light curves extracted from the blue spectra of the August 1994 outburst (see section 3). The right-hand panels show the asymmetric part of the resulting eclipse maps in logarithmic greyscale for S/N=50 and 25, respectively. The asymmetric component of each eclipse map is computed as follows. The radial intensity profile is separated into 9 radial bins of equal width. We then sort the intensities in each bin, and fit a smooth spline function to the lower quartile of the intensities in each bin. The fit is then subtracted from each pixel and the resulting map is written as the asymmetric component. We practically remove the baseline of the radial profile, leaving all azimuthal structure in the asymmetric map. We further set to zero all negative intensities in the asymmetric map.

The quality of the reconstruction is sensitive to the phase resolution of the data. The light curve with the highest phase resolution and best S/N shows that the spiral structure is considerably blurred to a “butterfly” pattern, as a consequence of the azimuthal smearing effect of the eclipse mapping method. Poor phase resolution ($`\mathrm{\Delta }\varphi `$=0.015) results in spiral arms with poorer definition, since the gradient change in the eclipse light curve shape – which testifies the presence of the spirals – is poorly sampled. The noise masks the ingress/egress of the spiral arms (gradient changes in the eclipse light curve shape) and the eclipse mapping method produces a V-shaped fairly symmetric model light curve of constant ingress/egress slope which will drive the solution to a more symmetric configuration by smearing the spiral arms. Thus, the spirals cannot be resolved on the map if the light curves have low S/N (or are affected by significant disk flickering).

### 2.3 Accurate zero phase of mid-eclipse ($`\mathrm{\Delta }\varphi _0`$)

In Fig. 4 we apply additional phase offsets to the eclipse light curves in order to explore the effect on the spiral arms of the uncertainty on the exact zero phase (0.006 cycles for IP Peg; see section 3). Small errors in the mid-eclipse, $`\mathrm{\Delta }\varphi _00.005`$ cycles, lead to distortions in the eclipse map by changing the azimuthal position and extension of the spiral arms. Positive (negative) phase offsets displace both arms to the trailing (leading) side of the primary Roche lobe and also reduces the strength of the blue (red) arm, thus severely distorting the disk light distribution. A phase offset of -0.010 to the light curves is sufficient to push the light towards a single segment of the outer disk, in fact where the blue spiral arm lies at. This exercise emphasizes the importance of knowing the value of $`\varphi _0`$ with good accuracy.

### 2.4 Spiral structure dilution by the symmetric disk component

Fig. 5 presents the test results of the spiral structure light dilution caused by a symmetric disk component. The left-hand panels show model light curves (with phase resolution of $`\mathrm{\Delta }\varphi =0.01`$ and S/N=50) for disk models ($`\theta =46\mathrm{°}`$) in which the spiral structure contributes 40, 30, 20 and 10 per cent of the total flux. The right-hand panels show the asymmetric components \- after subtracting the symmetric map from the total intensity map - of the corresponding eclipse maps for S/N=50 and 25, respectively. At small fractional contributions (10 %), the spiral structure is significantly diluted by the superimposed symmetric disk light and the “butterfly”-like pattern (bottom panels) is smeared to a “crescent”-like pattern at the outer disk facing the $`L_1`$ point. Provided that the spiral structure contributes a significant fraction to the total light ($`f_{\mathrm{spiral}}/f_{\mathrm{total}}>0.2`$), the two-arm structure is then reasonably recovered as a “butterfly” pattern.

Estimating the contribution of the spirals from the ratio of the flux in the asymmetric map component and that of the total eclipse map may lead to systematically underestimated fractional contribution. The reason fot this is that part of the flux from the spiral arms goes into the symmetric component because of the azimuthal smearing effect of the eclipse mapping method. The error increases for light curves of lower phase resolution: the fractional contribution of the spiral arms is better recovered when the eclipse light curve shape is better sampled. Using our simulations, we determine that the inferred fractional contribution of the spirals (40%) may be underestimated by a factor of 1.6-2.4 for light curves (S/N=25) with phase resolution of $`\mathrm{\Delta }\varphi =0.0100.015`$ orbital cycles (as observed), respectively. Light curves with S/N=50 and $`\mathrm{\Delta }\varphi =0.005`$ can recover the fractional contibution of the spiral structure within 12%. At small fractional contributions ($`<`$10%), the spiral structure will start to be overestimated (30-40 %) by the technique. The error in the fractional contribution of the spirals also depends on the S/N of the light curve but is much smaller than that caused by low phase resolution. Thus, since we estimated a fractional contribution of 0.15–0.3 in Baptista et al. (2000), the true contribution of the spiral arms in those eclipse maps may be in the range 0.3–0.6. Hence, the models with $`f_{\mathrm{spiral}}/f_{\mathrm{total}}0.30.4`$ should be a good approximation of the real case. In conclusion, we can only place a lower limit to the fractional contribution of the spiral arms to the total light using eclipse maps.

### 2.5 Brightness difference between the two spiral arms

Finally, we undertake a test with a spiral structure whose arms are of unequal brightness. This is driven by the fact that observed Doppler maps show an intensity difference between the two spiral arms (Harlaftis et al. 1999; Steeghs 2001). Fig. 6 presents the effect of brightness difference between the two spiral arms on the map reconstructions. We build model light curves ($`\mathrm{\Delta }\varphi =0.01`$, S/N=50) from spiral structure/disk models, as described in Fig. 2 in which the red and blue arms have equal brightness (1:1, top-left panel), or the red arm is fainter - as observed - with a brightness ratio to the blue arm at 1:1.5 (middle-left panel) and finally at 1:2 (bottom-left panel). The other panels in the figure show the corresponding eclipse maps in logarithmic grayscale for three different S/N ratios (50, 30 and 20). It is notable that both arms appear in the map as a “butterfly” pattern when they have equal brightnesses regardless of the S/N ratio. However, the brightness difference between the spiral arms at low S/N ratio data ($`<30`$) leads to loss of the fainter spiral arm in the image reconstruction.

## 3 Spiral shock observations

### 3.1 Eclipse maps of the August 1994 outburst

Morales-Rueda et al. (2000) present Doppler tomography showing persistent spiral shocks in the accretion disk of IP Peg during the August 1994 outburst. The reader is referred to Morales-Rueda et al. (2000) for the details of the observations and data reduction procedure. In summary, 164 blue spectra were obtained at a spectral resolution of 100 km s<sup>-1</sup> at H$`\beta `$ and a time resolution of 220 seconds covering the emission lines H$`\beta `$, H$`\gamma `$, He II 4686 Å (4040–4983 Å). Here, we extend their analysis in order to include the eclipse maps for comparison to the Doppler maps. We combined the four eclipse light curves (E= 25158, 25159, 25164 and 25165) from the two consecutive nights the observations took place (30/31 August 1994) to produce average eclipse light curves with phase resolution of 0.010 cycles. In particular, we extracted light curves of the blue continuum and of the emission lines H$`\gamma `$, He II 4686 Å, the Bowen blend and the He I 4472 Å (Fig. 5 in Morales-Rueda et al. 2000) in order to reconstruct the images of the accretion disk.

We processed the light curves before applying the eclipse mapping procedure so that the continuum level before and after eclipse is the same; otherwise, artifacts are raised. These low-amplitude orbital variations were removed from the light curves by fitting a spline function to the phases outside eclipse, dividing the light curve by the fitted spline, and scaling the result to the spline function at phase zero. This procedure removes orbital modulations, such as the bright spot, with only minor effects on the eclipse shape itself (e.g., Baptista et al. 2000). The other factor that affects the eclipse mapping procedure (see simulation at Fig. 4) is an accurate ephemeris for the eclipse minimum. Here, we use the linear ephemeris by Wolf et al. (1993) after correcting the phase zero for mid-eclipse rather than the white dwarf egress (the white dwarf eclipse width is $`\mathrm{\Delta }\varphi =0.086`$ cycles; Wood & Crawford 1986)

$`T_{\mathrm{mid}}(HJD)=2445615.4156+0.15820616E`$

where $`T_{\mathrm{mid}}`$ is the inferior conjunction of the red dwarf. IP Peg shows cyclical orbital period changes of amplitude 0.002 days (or 0.013 orbital cycles; Beekman et al. 2002). Furthemore, the scatter of eclipse timings for a given epoch is also about 0.001 days (0.006 cycles). Clearly, even a sinusoidal ephemeris cannot predict eclipse timings with precision better than 0.01-0.02 cycles for IP Peg (Wolf et al. 1993). Therefore, the best approach in order to properly phase the light curves is to center the white dwarf eclipse using contemporary quiescent eclipse timings. The August 1994 data (eclipse cycle $`E25160`$) is luckily bracketed by two HST sets of observations of IP Peg in quiescence (Baptista et al. 1994). From these timings, we infer that the white dwarf mid-eclipse at that epoch occurred 0.008 orbital cycles before the prediction of the linear ephemeris of Wolf et al (1993). Hence, we added +0.008 cycles to the phases of the light curves to make phase zero coincident with the inferior conjunction of the red dwarf. The minimum of the eclipse in the lines is displaced towards earlier phases in comparison to the continuum light curves. Eclipse maps were then made with three different geometries: ($`i,q`$)=($`81\mathrm{°}`$, 0.5) (Wood & Crawford 1986), ($`i,q`$)=($`79.4\mathrm{°}`$, 0.58) (Marsh 1988) and ($`i,q`$)=($`86.5\mathrm{°}`$, 0.31) (Beekman et al. 2000), where $`i`$ is the inclination and $`q`$ is the mass ratio. These tests do not show any difference caused by a slight change in the adopted orbital parameters. Thus, we adopt in the following the geometry of Wood & Crawford (1986).

All maps in Fig. 7 show an asymmetric disk. In order to resolve the asymmetric part of the disk we subtract the symmetric part of the disk and plot the residual on the right panels of the figure. The continuum eclipse is skewed towards the egress and maps into an asymmetric disk which is brighter towards the gas stream and $`R_{L1}`$. The Bowen blend and He II high-ionization lines show a relatively narrow, asymmetric eclipse shape which maps into asymmetric brightness distributions consistent with spiral structure and similar to those in Baptista et al. (2000, see their Fig. 2). However, the Balmer and He I line maps are markedly different. The strong asymmetric Balmer emission is located close to and even outside the tidal radius, with the maximum emission at an angle of $`30\mathrm{°}`$ with respect to the line-of-centres. The centre of the He I eclipse occurs before that of the Balmer lines, which shifts the asymmetric emission towards a direction opposite to the gas stream (blue-shifted velocities of the disk), and beyond the tidal radius. We apparently see only one of the arms in He I and there is almost no evidence of the spirals in the Balmer lines – probably because their emission is significantly diluted by other light sources in these lines (see following section on velocity-resolved light-curves). The “crescent”-shape of the brightness distribution of the Balmer lines is reminiscent of the distributions obtained for simulations where the spirals contribute only 10 per cent of the light at the corresponding wavelength (see Fig. 5). In comparison, Morales-Rueda et al. (2000) use the out-of-eclipse spectra from the same orbital cycles to easily resolve the spiral shocks using Doppler tomography.

### 3.2 Eclipse maps from July 1998 outburst

We also used WHT service time to obtain phase-resolved spectrophotometry with the ISIS spectrograph of an eclipse of IP Peg during its July 1998 outburst decline, four days after the peak of the outburst (18 July 1998). The 23 blue spectra (and 31 red spectra) covered the range 4407.6–4847.7 Å (and 6310.7–6717.6 Å) at a spectral resolution of 0.2 Å per pixel (and 0.8 Å, respectively) and with an exposure of 120 seconds. The spectra covered the eclipse between 0.85 to 1.2 cycles at a phase resolution of 0.011 cycles for the red spectra and 0.014 cycles for the blue spectra. The reconstructed eclipse maps show similar distributions to the August 1994 maps, but with inferior quality which is most likely due to the poorer phase resolution and the later outburst stage (closer to quiescence which should result in a fainter spiral structure). In summary, the eclipse maps are consistent with the reconstructions from the August 1994 data (previous section) and the August 1996 data (Baptista et al. 2000).

### 3.3 Velocity-resolved eclipse maps : Revealing the spiral shocks

We investigate here velocity-resolved eclipse maps since it is a well established fact that there are non-disk emission components to the total line fluxes, for example from the donor star during outburst (Harlaftis 1999). During outburst, another contaminating line source to the disk emission is a low-velocity emission component (Harlaftis et al. 1999; Steeghs et al. 1996). Since the spiral arms and the secondary star contribute to the line profile at different Doppler velocities, it is possible to separate the contribution of these distinct light sources by suitable extraction of velocity-resolved eclipse light curves. Baptista et al. (2002) performed velocity-resolved eclipse mapping to separate the He I emission from each of the spiral arms. Here, we use this technique as a useful tool in order to separate the contribution of the low-emission components (secondary star, centre-of-binary emission; Harlaftis et al. 1999) from the disk emission and, thus, make the spirals more easily detectable in the line eclipse maps. For example, we can extract eclipse light curves from the blue (red) velocities alone. The application of the eclipse mapping technique results in the blue (red) part of the accretion disk (iso-velocity surfaces on the disk due to Keplerian rotation). The sum of the blue and red parts is indeed, as expected, consistent with the total image as built by the velocity-integrated light curve. After this test, we proceed with the aim to isolate the inner/outer disk by extracting eclipse light curves from different velocity ranges. For this purpose, we extracted Balmer and He I light curves for the line core (-250 to +250 km s<sup>-1</sup>) and for the line wings (-1300 to -700 and +700 to +1300 km s<sup>-1</sup>). The eclipse light curves from the core of the Balmer lines show a broad eclipse $`\mathrm{\Delta }\varphi >0.2`$ whose minimum is displaced towards earlier phases with respect to the continuum mid-eclipse. The eclipse light curves from the wings of the Balmer lines are centred at phase zero and show the same asymmetric eclipse shape to the He II and Bowen blend light curves with a clear step during ingress.

Fig. 8 illustrates the effects of low-velocity emission components in diluting the spiral structure embedded in the accretion disk. The figure presents the results of eclipse mapping of the H$`\gamma `$ line, from top to bottom, for the velocity-integrated total profiles (full line), the sum of the core/wings velocity regions, the core of the profiles, and the wings of the profiles. The reconstruction from the “core” light curve has a broad asymmetry, towards the $`L_1`$ point, in the outer regions of the eclipse map. The eclipse map from the “wings” light curve clearly shows the two arc-shaped spiral arms at the same orientation and azimuthal extent as are the ones in the Bowen blend and He II eclipse maps (Fig. 7). The “core/wings” light curve (and corresponding eclipse map) is similar to the velocity-integrated light curve (and eclipse map). Similar results are found for the H$`\beta `$ and He I lines. This phase offset from mid-eclipse and the resulting asymmetries in the eclipse maps are seen in the Balmer and He I lines (which also have a substantial contribution from the secondary star at outburst), but are not present in the He II and Bowen blend light curves and eclipse maps (the contribution of the secondary star to these lines is weak or negligible). The eclipse mapping technique, in order to account for the wide and offset eclipse, tries to build a brightness distribution in the disc regions being “eclipsed”, irrespective if the physical origin lies on the secondary star or a disk wind, thereby producing a fake, broad asymmetry in the outer regions of the eclipse map displaced towards the $`L_1`$ point.

On the right-hand panels, we present the S/N ratio of the maps built. For each eclipse light curve, a set of 15 artificial curves was generated, in which the data points were independently and randomly varied according to a Gaussian distribution with standard deviation equal to the uncertainty at that point. The artificial curves were fitted using the eclipse mapping technique to produce a set of eclipse maps, based on random data. These, then, were combined to produce a “residuals” map by taking the pixel-to-pixel standard deviation with respect to the mean intensity. This yields the statistical uncertainty at each pixel. The uncertainty of the standard deviation for N=15 is $`<`$ 20% which is sufficient for an illustration of the confidence limits of the eclipse map (more than an order of magnitude of Monte Carlo simulations are needed in order to reduce the uncertainty of the standard deviation down to 7 %). The result of dividing the eclipse maps of the observed data (or the mean map of the Monte Carlo simulations) by the “residuals” map is shown in the right-hand panels of Fig. 8, as maps of the inverse of the relative error or S/N ratio maps. Darker regions have smaller uncertainties (greyscale is such that pixels with S/N $`5`$ are shown in white color, while all pixels with S/N$`>`$100 are shown in black color). The spiral arms in the H$`\gamma `$ “wing” map have typical S/N = 12–20. The H$`\gamma `$ map (full line) has the best S/N because it was integrated over the whole line profile.

There may be various sources of low-velocity emission in the system to produce the above effects on the eclipse shape. One such emission component of unknown origin but close to the binary centre and slightly towards the gas stream has been seen at outburst maximum in IP Peg (Harlaftis et al. 1999; Steeghs et al. 1996). Wind-related emission has also been suggested by Baptista et al. (1995) in order to accommodate excess light that was not eclipsed by the companion star. The velocity of the red star, as measured from different investigators, is also low either at 280, or 300 or 330 km s<sup>-1</sup> (Wolf et al. 1998; Marsh & Horne 1990; Beekman et al. 2000, respectively). Irradiation of the hemisphere of the star facing the white dwarf during outburst causes a significant contribution to the Balmer and He I total emission (Morales-Rueda et al. 2000), and even in He II at outburst maximum (Harlaftis et al. 1999). The implication is that if the emission pattern from the Roche lobe is not symmetric (possibly due to shielding by the bright spot and gas stream), the centre of the self-eclipse of the secondary star line emission happens earlier than the disc eclipse at phase zero. However, the data from Morales-Rueda et al. (2000) do not show any solid evidence for an asymmetric emission distribution on the Roche lobe of the companion star.

## 4 Discussion

We have performed a thorough investigation to determine the conditions under which it is possible to detect spiral structure in accretion disks of cataclysmic variables, and then compared the results with observed data. The simulations (see Table 1 for a summary) show that the phase resolution combined with the S/N ratio are critical in resolving the features caused by the spiral shocks in the eclipse light curve shape. Indeed, poor phase resolution ($`\mathrm{\Delta }\varphi >0.01`$) and low S/N ratio ($`<25`$) can distort the spiral structure beyond visible recognition (bottom panels in Fig. 3). In addition, the eclipse mapping technique assumes a uniform, flat accretion disk and this smears the spiral pattern into a “butterfly” pattern (top panels in Fig. 3).

Aside for the above observational parameters (time resolution and S/N) and technical effects (maximum entropy method and model assumptions), the emissivity parameters of the spiral structure are equally important in order to reveal it using the eclipse mapping technique. Dilution by the accretion disk light makes it difficult to resolve the spiral structure. The dilution of light ($`>`$ 70%) can distort, through the eclipse mapping technique, the smeared spiral structure (“butterfly” pattern) to a “crescent”-like structure (bottom panels in Fig. 5). In addition to the above, the comparative emissivity between the two spiral arms may not be the same, as indeed observations indicate (“red”-arm weaker; Harlaftis et al. 1999, Steeghs 2001). This also contributes to further attenuation of the reconstructed “red” arm (bottom panels in Fig. 6). Poor phase resolution itself may result in loss of the red spiral arm even if its brightness is equal to that of the other arm. Having good phase resolution is even more important than having high S/N light curves. The observed maps presented here from the early decline of the August 1994 outburst (but also similar results from spectra at late decline of the July 1998 outburst) can be more easily read given the simulations. They are consistent with reconstructing the spiral structure as a smeared and distorted spiral structure, and frequently one-arm pattern in the accretion disk. Similarly, the maps presented by Baptista et al. (2001) at outburst maximum and by Baptista, Haswell & Thomas (2002) at outburst decline can be interpreted in the same fashion. The WHT observations (July 1998 outburst) show an application for the case where the data does not match the requirements put forward by the simulations, and the INT observations (August 1994 outburst) is an example of a data set that meets the observational requirements.

The presence of spiral arms are inferred by small asymmetries in the shape of the eclipse light curve. When the noise in the light curve becomes comparable to the depth of these asymmetries, the method looses the ability to distinguish (or disentangle) the spiral structure from the noise. Also, the characterization of the spiral structure in the eclipse maps demands that the corresponding asymmetries in eclipse shape are properly sampled in phase. For light curves of poor phase resolution (or, more specifically, when the number of data points sampling the asymmetry caused by the spiral arm is small) it is hard to distinguish the true asymmetry from random deviations in the eclipse profile caused by noise. These effects become more important if the contribution of the spirals is diluted by other sources such as a symmetric disc component, because the relative depth (and significance) of the asymmetries decrease for increasing brightness of the additional source.

In the absence of dilution by other light sources, as with the Bowen blend light curve (Fig. 7), the spiral structures may be well recovered even for light curves of relatively poor phase resolution and incomplete phase coverage such as those of Baptista et al. (2000; see their Fig.2). In the presence of light sources diluting the spiral structure, optimal extraction of eclipse light curves in narrow passbands matching the velocity range of the spiral arms is an interesting way to optimize the ability of the eclipse mapping method in recovering these asymmetric structures. Indeed, this is the case with the eclipse light curves of the Balmer and helium lines (Fig. 8) where an unknown source of asymmetric light shifts the mid-eclipse of the emission lines resulting in distorted eclipse maps (“crescent”-shape or appearance of a “single spiral-arm”) with the surface brightness distribution towards the companion star. Suitable velocity-range selection from the line profiles can avoid any contamination to the spiral structure and thus resolve clearly the spiral structure in the accretion disk. Such is the case, for example, with the eclipse light curve extracted from the “wings” of the emission-line H$`\gamma `$ profiles and the corresponding eclipse map in Fig. 8. Given the many parameters involved, it is not a surprise that the spiral structure was first discovered with the use of Doppler tomography. Doppler tomography defines very clearly the spirals since it provides more constraints from phase-resolved spectroscopy (velocity, intensity, phase).

We conclude that the eclipse mapping method is able to properly reproduce the spiral structure provided that the observations are tuned (S/N$`>`$25, and $`\mathrm{\Delta }\varphi `$=0.01 or a time resolution of at least 2 minutes for IP Pegasi), the data analysis is careful (avoiding any contaminating low-velocity emission by selection of suitable velocity-resolved light curves), and the physical characteristics of the spiral arms are sufficient for detection (contribution of $`30`$ % of the total disk light and minor brightness difference between the two spiral arms).

The analysis we have undertaken using eclipse maps from real and synthetic data demonstrate how future observation runs can be tailored in order to deduce in full detail the significance of the spiral structure in fast-rotating cataclysmic variable disks. The presence of spiral structure suggests that the application of eclipse mapping techniques could benefit from the use of a default model disk which smears structures along spiral patterns instead of the usual azimuthal direction (Harlaftis et al.,2004, in preparation).

###### Acknowledgements.

The INT and WHT telescopes are operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. The WHT spectra of IP Pegasi were obtained through the SERVICE programme of the Isaac Newton Group. We thank an anonymous referee for his comments.