text:
# The HERA-B Ring Imaging Čerenkov Counter

## 1 Introduction

HERA-B, a fixed target experiment (see Fig. 1) at the HERA storage ring at DESY, was designed  to measure rare processes in the decays of $`B`$ mesons. The $`B`$ mesons are produced in collisions of 920 GeV/c protons with a fixed target, which consists of 8 wires which can be individually inserted into the halo of the proton beam in order not to disturb experiments measuring $`ep`$ collisions. One of the essential components of the spectrometer is the Ring Imaging Cherenkov counter (RICH) . The main purpose of the RICH counter is the identification of charged hadrons, in particular kaons from decays of $`B`$ mesons. Identifying charged kaons essentially means separating them from pions in the momentum range between 3 GeV/c and about 50 GeV/c at an interaction rate of up to 40 MHz.

This paper is organized as follows. We shall first discuss the design criteria and their implementation. The components of the counter will be presented, as well as the monitoring and analysis programs. Finally, results of measurements made with the counter will be discussed, and the performance of the RICH as an identification device, and also as a tracking device, will be presented.

## 2 Design Criteria

In a fixed target experiment such as HERA-B where one has to deal with an intense flux of charged particles, the design choices of a ring imaging Čerenkov counter are governed by the following criteria. To achieve the necessary performance, a sufficient number of Čerenkov photons (at least 20) has to be detected per ring image of a $`\beta =1`$ particle . This requirement fixes the length of a gas radiator to a few meters. The high rate capability of the photon detector requires a detector with pad readout, which limits the resolution to a few millimeters. This, in turn, together with the required high resolution in the measurement of the Čerenkov angle for single photons (better than 1 mrad) sets the mirror focal length value to several meters. Last but not least, the requirement that single photons are efficiently detected with low background, forces the photon detector to be kept out of the solid angle for charged particles and conversion products.

## 3 Implementation

### 3.1 Radiator

Perfluorobutane gas ($`\mathrm{C}_4\mathrm{F}_{10}`$) was chosen for the HERA-B RICH radiator, since it combines a relatively high refractive index, $`n=1.00137`$, and low dispersion (5% variation in $`n1`$ over the spectral acceptance from 300 nm to 600 nm ). In this gas the Čerenkov radiation threshold momenta for pions and kaons are 2.7 GeV/c and 9.6 GeV/c, respectively. For $`\beta =1`$ particles, the Čerenkov angle is 52.4 mrad, while the $`\pi \mathrm{K}`$ difference in Čerenkov angle is 0.9 mrad at 50 GeV/c. The r.m.s. angular spread due to dispersion amounts to 0.33 mrad. Since some freons are known to scintillate considerably, (e.g. CF4), we have studied C4F10 scintillation  and found that this contribution to background is less than 0.2 detected photons per charged particle.

The radiator vessel, made from stainless steel plates with 1 mm aluminum particle entrance and exit windows (see Fig. 2), is placed about 8.5 m downstream of the target (see Fig. 1). Two beam shrouds made of carbon fiber-reinforced epoxy resin CFK close the gas volume around the two beam pipes for protons and electrons. The Čerenkov light exits the vessel through 2 mm thick UV grade Plexiglass windows . The vessel is filled with 108 m3 of C4F10 gas, which is being circulated in a closed system with liquefaction stages for cleaning and buffering . Requirements on impurities are moderate due to operation in the visible and near UV part of the Čerenkov spectrum. The radiator gas is kept at 2.5 mbar overpressure relative to ambient pressure. Additionally, a bubbler with a 10 mbar pressure difference threshold is attached to the vessel.

### 3.2 Mirrors

The main imaging device is a spherical mirror placed inside the radiator vessel with the center of the sphere near the target and a radius of curvature of 11.4 m. The mirror, a 6 m by 4 m rectangular cutout from the sphere, consists of 80 full or partial hexagons made from 7 mm thick Pyrex glass, coated with 200 nm of aluminum and 30 nm of MgF2 (see Fig. 3). To be able to place the focal surface outside the particle flux ($`\pm `$160 mrad vertically), the mirror is split horizontally and both halves are tilted by 9 away from the beam-line. A set of two planar mirrors, composed of 18 rectangular elements each, translates the focal surface to the photon detector area above and below the radiator vessel (see Fig. 2). The planar mirrors are made of float glass, thus being significantly cheaper than Pyrex mirrors at the required optical quality. All 80 spherical and 36 planar mirrors are mounted on a rigid, low mass support structure inside the radiator volume and can be individually adjusted by stepper motors from the outside.

The mirror quality was determined upon delivery by measuring for each segment the radius of curvature and the fraction of reflected light, and by recording a Ronchi image to check the homogeneity of the mirror surface . The reflectivity was required to exceed 85% in the wavelength interval from 250 nm to 600 nm. From the data gathered, the mirror segments were grouped in a tiling scheme according to their optical quality and resolution requirements .

All mirrors were aligned to better than the required precision by surveying them inside the vessel. An off-line data based algorithm was developed for monitoring of the positions of the mirror segments during the measurement .

### 3.3 Photon detector

Two types of wire chamber based photon detectors were initially considered, a CsI photocathode in a MWPC, and a TMAE detector with 10 cm deep, 8 x 8 mm2 unit cells. After considerable success of on-the-bench and beam tests , the detectors were tested in a high rate environment, as expected in the HERA-B experiment. Both detectors had to be abandoned; the TMAE detector showed a prohibitive decrease of avalanche gain due to aging effects , while the CsI photocathode could not be routinely produced and maintained with sufficiently high quantum efficiency, in addition to problems with rates in excess of a few kHz per pixel .

The actual photon detector consists of Hamamatsu multi-anode R5900-00-M16 and R5900-03-M4 photomultiplier tubes (PMTs). In what follows they are denoted by M16 and M4. The M16 version has 16 pads of 4.5 x 4.5 mm2 each, with a 12-stage, metal-foil dynode system . The M4 version has 4 pads of 9 x 9 mm2 each, and 10 dynodes. The quantum efficiency of the M16 photocathode with borosilicate window has a broad plateau in the wavelength region between 300 nm to 500 nm and a maximum value about 20%. The M4 tubes have a UV transparent window, which shifts the low wavelength cut-off to about 250 nm. The other PMT characteristics, such as the required cathode high voltage ($``$ 1000 V), current amplification (107), dark current ($``$ 1 nA), pulse rise time (0.8 ns), transit time spread (0.3 ns) are also satisfactory . In the initial set of on-the-bench measurements , the single photon counting properties of the photomultiplier tube were investigated, in particular the efficiency for single photon detection as well as the background count rate. It was established that the PMT allows for a good single photoelectron detection efficiency (above $`98\%`$), with very small cross talk (below $`0.2\%`$), low background rate (few Hz per channel, Fig. 4), and acceptable uniformity as shown in Figs. 5 and 6 .

The PMTs have the outer dimension of $`28\times 28`$ mm2 with an active area of $`18\times 18`$ mm2. They are positioned on a $`36\times 36`$ mm2 grid. To increase the fraction of the active area, a two lens demagnification system (2:1) was designed, as shown in Fig. 7 . The lenses are made of UVT perspex  with high transparency over most of the wavelength region where the photocathode is sensitive (Fig. 8). The angular acceptance of the optical system is also satisfactory and is uniform for incident angles below about 110 mrad . The increase in the photodetector active area is achieved at the expense of increasing the pixel size from the 4.5 x 4.5 mm2 PMT pad size to 9 x 9 mm2 pixel in the central detector region. The resulting angular measurement error (0.46 mrad) slightly exceeds the spread due to dispersion (0.33 mrad). For the outer detector region, with lower track densities and lower typical track momentum, the photon flux is smaller, and the resolution requirements looser. In order to reduce the detector cost, this detector part uses the M4 version of the tube with two times larger pads (9 x 9 mm2). The same lens system thus results in 18 x 18 mm2 pixels which correspond to an angular measurement error of 0.92 mrad.

In order to reduce the contribution of spherical aberration to the overall resolution of the Čerenkov angle, an optimal surface of the Čerenkov photon detector was calculated . Each half-detector (upper and lower) consists of 7 flat supermodules placed to approximate the optimal surface, which is close to the shape of a flattened (ellipsoidal) cylinder. Such an arrangement also ensures better acceptance for the Čerenkov photons, which should be incident onto the flat supermodules at angles below 110 mrad. The supermodules are rectangular 1.1 $`\times `$ 0.4 m2 boxes containing a grid made from 1 mm thick soft iron sheets that serve as magnetic shield and mounting structure for the PMT base-boards and the light collection system. The two types of photomultipliers were arranged according to occupancy and reconstruction requirements as indicated in Fig. 9. Altogether 1488 M16 PMTs and 752 M4 PMTs are needed to cover the detector surface, thus totaling 26816 read-out channels.

Prior to the installation in the photon detector tests of all the 2305 acquired PMTs (1543 M16’s and 762 M4’s) were made in the laboratory . Čerenkov photons produced in a 1 cm thick quartz window by $`\beta `$ particles from a 90Sr source were used as a stable light source . Four PMTs were tested simultaneously, of which one was used as a reference. For each photomultiplier tube, the source and background rates were recorded as a function of high voltage and threshold setting. On the basis of these tests, the photomultipliers were grouped according to similar high voltage characteristics, allowing all PMTs within a group to be connected to the same high voltage, thus minimizing the number of independent HV channels . The results obtained for the optimal high voltage as well as for the relative PMT sensitivity were compared to values provided by the manufacturer, and good agreement was found . The count rates of the reference M4 and M16 photomultipliers, recorded during the quality assessment tests and also continually measured after the tests, were used in order to obtain an estimate of the long term stability. The count rate decrease in two years of operation is consistent with the known decay rate of the 90Sr source.

### 3.4 Read-out and monitoring

To mount, power, and read out the photomultipliers, two types of base-boards were developed, one for M4 and one for M16 PMTs. Each base-board is a light-tight, four-layer (M16) or two-layer (M4) circuit-card, 70 $`\times `$ 70 mm2 in size, equipped with surface mount components. One side of the board holds four custom designed sockets and voltage divider chains, appropriate for the particular type of PMT; the other side contains one (M4) or four (M16) connectors for the 16-channel readout-cards, attenuation circuitry, and two daisy-chainable high voltage connectors. The voltage divider chain is composed of a sequence of 1 M$`\mathrm{\Omega }`$ and 560 k$`\mathrm{\Omega }`$ resistors that was optimized for single photon counting at high rates . The base-boards are mounted in such a way that the PMTs are positioned at centers of the $`36\times 36`$ mm2 grid (see Fig. 10).

Between two and five base-boards, with their PMTs selected for similar optimal operating voltage (ranging from 750 V to 890 V), are powered by a single output of a 160-channel CAEN SY527 high-voltage supply .

The front end electronics board employed in the system is a 16-channel board  based on the ASD8, an amplifier, shaper, and discriminator chip . The board was developed for the Outer Tracker (OTR), the main tracking system of the HERA-B spectrometer. Since the photomultiplier tube signals are typically considerably higher than the pulses for drift chambers, a charge divider was added in front of the ASD8 board (1:10 for M16, and 1:5 for M4 PMTs). The front end electronics boards were arranged according to their optimal threshold voltage by using test data which were obtained by varying pulse height and threshold voltages . The same test was also used to reject faulty boards.

The low voltage supply system for the ASD8 boards is powered by four KEPCO ATE 6-100M supplies, delivering voltages of $`+3\mathrm{V}`$ and $`3\mathrm{V}`$ separately to the upper and lower photon detector halves . Analog outputs of 16 photon detector channels are transmitted to the electronics hut for diagnostic purposes.

The signals from the ASD8 cards are transmitted over 7.5 meter long 16 channel twisted pair cables to the Front End Drivers (FEDs). For each FED set there are 4 daughter cards with 16 cable connections, and one mother card, such that each FED set can accomodate up to 1024 channels . The RICH uses 28 such sets. The FED mother board provides the interface to the HERA-B Fast Control System, which provides a clock signal for synchronization as well as triggers to indicate event forwarding. The mother board also provides the interface to the HERA-B DAQ System via the SHARC protocol , as well as a ring buffer memory capable of storing 128 events, pending the First Level Trigger decision.

A sizable fraction of the operating parameters of the RICH counter (low voltage and threshold for the front end electronics boards, temperatures at the photon detector and power supplies) are controlled by a Field Point  based system. Two separate systems are used for monitoring the high voltage  and steering the gas system . The overall control of these systems is implemented with a server-user interface system. The server process is running on a VME computer connected to the Field Point system through a RS232 cable. In the case of a malfunction (if the temperature, threshold, or low voltage are out of range) the server switches off the low voltage power supplies. The parameters can be monitored and set via a tcl/tk and BLT based graphical user interface program running on a linux PC and communicating with the server through a TCP-IP based protocol. The RICH low voltage and high voltage slow control systems are integral parts of the overall HERA-B slow control system .

As a triggerable light source for testing the photon detector in the periods without beam induced reactions, a system of light emitting diodes (LED) was installed at the sides of the photon detector windows . The system provides good illumination over the whole photon detector area.

### 3.5 Expected performance

The number of detected photons on a ring measured by a Čerenkov detector is often parametrized as $`N=N_0L\mathrm{sin}^2\theta _c`$, where $`L=2.82`$ m is the length of the radiator, $`\theta _c`$ is Čerenkov angle, and $`N_0`$ is the detector response parameter. From the data available on the quantum efficiency , mirror reflectivities , window and optical system transmissions , one calculates the merit factor $`N_0=41`$ cm-1. Accordingly, the expected number of photons for particles approaching the velocity of light ($`\theta _c=52.4`$ mrad) amounts to 32.

Individual contributions to the expected single photon resolution are summarized in Table 1. The optical error includes contributions from spherical aberration, mirror quality and alignment of mirrors, while photon detector resolution is dominated by detector granularity with a small contribution of the alignment of photon detector components. The contribution of multiple scattering in the RICH counter  depends on momentum and is given by $`\frac{3.5\mathrm{mrad}}{p(\mathrm{GeV}/\mathrm{c})}`$. The resulting single photon resolution, $`0.65\mathrm{mrad}\frac{3.5\mathrm{mrad}}{\mathrm{p}(\mathrm{GeV}/\mathrm{c})}`$ ($`1.02\mathrm{mrad}\frac{3.5\mathrm{mrad}}{\mathrm{p}(\mathrm{GeV}/\mathrm{c})}`$) for the regions covered by M16 (M14) PMTs, does not include the contribution from the track direction as given by the tracking system, and will therefore be referred to as intrinsic resolution.

## 4 Data Analysis

Several computer programs were developed for the analysis and monitoring of the data collected by the RICH counter. In the initial phase of the commissioning of the HERA-B spectrometer, when major parts of the main tracking system were not yet available, a stand-alone ring search algorithm was employed . The program searches for rings with a given ring radius, starting at the maximum value and scanning down to $`70\%`$ of the maximal ring radius. Only rings with a sufficiently high signal over background ratio are stored, and the corresponding hits are excluded from further ring search. As described in section 6, the reconstructed ring center gives information that can be used in tracking. The algorithm is routinely used to monitor the current maximal value of the Čerenkov angle  which is, in turn, used as the input for all particle identification algorithms. Small variations of refractive index and thus of the Čerenkov angle are caused by pressure and temperature variations of the radiator gas. We note that several alternative ring search algorithms were also tested on the data .

For highly populated events, typical for the reactions recorded in the HERA-B spectrometer, additional information helps to improve the particle identification capabilities. The information from the tracking system, track direction and momentum, is used to determine the expected position of the ring center and the expected ring radius for various particle hypotheses. To evaluate the likelihood for each individual hypothesis, extended maximum likelihood methods  are applied, either directly  or combined with an expectation-maximization algorithm . The latter is described in more detail in section 5.

Another way to take the track information into account is to perform a seeded ring search in regions given by the track direction and momentum . Since only restricted regions corresponding to five hypotheses are scanned, this approach is quicker than the stand-alone ring search. Recently, a new particle identification method was developed which maximizes the likelihood value by allowing the track direction to vary within errors given by the track fit .

## 5 Measurements and Results

### 5.1 Photon detector performance

In the first step of the commissioning of the system, some basic parameters were investigated. Fig. 11 gives the count rate versus high voltage for a representative photomultiplier. It is shown that the curves measured in-situ with the HERA proton beam agree nicely with the 90Sr source measurements. The time averaged occupancy is shown in Fig. 12, where the region occupied by M16 PMTs is clearly distinguished from the region occupied by M4 PMTs. The fraction of dead channels is about $`2\%`$, and is partly due to missing PMTs and partly due to PMTs which were damaged during the installation phase. Only about 0.3% of the 26816 channels were found to be noisy, and were excluded from the analysis. As a result, very clean rings can be observed in low multiplicity events such as the one shown in Fig. 13. The absence of random hits in the event also confirms the earlier result that the background due to scintillation in the radiator gas is negligible.

### 5.2 Number of photons

In order to obtain the basic RICH parameters from the data, initially only events with a small number of Čerenkov rings were analyzed . For each ring, the number of photon hits was counted and the Čerenkov angle was determined, obtaining a point in Fig. 14. It is seen that the observed rings follow the expected dependence of the number of photons versus the Čerenkov angle. The largest measured Čerenkov angle of 52.4 mrad also corresponds to the expectation for $`\beta =1`$ particles ($`\mathrm{cos}\theta _c=1/n`$). The parameter $`N_0`$ was calculated for each of the isolated rings and the distribution is plotted in Fig. 15, from which it is seen that the mean value is 42 cm-1. This value corresponds to 33 detected Čerenkov photons for $`\beta =1`$ particles.

### 5.3 Angular resolution

To determine the Čerenkov angle uncertainty given by a single Čerenkov photon, tracks with known identity are selected. For each track within a given momentum interval, one fills a histogram with the photon hits at their reconstructed Čerenkov angle, so that the hits due to Čerenkov rings should show up as a peak, as shown on Fig. 16. The distribution was fitted with a Gaussian function and polynomial background. The resulting momentum dependence of the single photon resolution is shown in Fig. 17. As expected, the single photon resolution depends on the particle momentum. For the high momentum tracks the resolutions of 0.81 mrad and 1.0 mrad were measured for the regions covered by M16 and M4 PMTs, respectively. The difference between the measured single photon resolution and the expected intrinsic resolution (as shown in Fig. 17) is attributed to the error in the measurement of charged particle trajectory, which was not taken into account in the calculation of the expected intrinsic Čerenkov angle uncertainty (Table 1).

The intrinsic resolution of the RICH counter was studied by the stand-alone ring search analysis , and turned out to be consistent with the expected values. The intrinsic resolution in Čerenkov angle given by a single Čerenkov photon was also obtained with the aid of hits in the electromagnetic calorimeter . For the selected set of tracks with cluster energy above 50 GeV, the center of the Čerenkov ring is accurately determined from the straight line connecting the center of the cluster with the wire target (the magnetic field was off for this measurement). The single Čerenkov photon resolutions obtained correspond to $`\sigma `$ = 0.7 mrad ($`\sigma `$ = 1.0 mrad) for the M16 (M4) PMTs, in good agreement with expected values given in Table 1.

The particle identification capabilities are affected by the high multiplicity of rings in a single event. Typically 50 overlapping rings are present in an event of which about 1/3 could be associated with measured tracks, while the rest belong to tracks coming from secondary interactions. As a result, channel occupancies in some regions are as large as 25%. We note that in such an environment the resolution in the Čerenkov angle measurement of a track does not simply scale with the inverse of the square root of the number of detected photons.

## 6 Particle Identification Performance

To illustrate the particle identification capabilities, a maximum likelihood fit was performed on Čerenkov angle distributions for individual tracks. Fig. 18 shows the resulting bands for different particles in the plot of reconstructed Čerenkov angle versus momentum, as obtained by analyzing a sample of measured data.

Two methods are routinely used for particle identification, a ring search based and a track based algorithm. The ring search based particle identification combines the stand-alone and seeded approach discussed in Section 3.5. In the first step, the stand alone reconstruction is attempted and rings found are matched with tracks. In the next step, tracks found by the vertex detector and the tracking detectors, that were not matched with a RICH ring, are used as seeds in a seeded ring search.

For the track based particle identification the extended likelihood method is used , combined with the expectation-maximization algorithm. The method is described in detail in . The first step in the method is the calculation of Čerenkov angle (as described in ) for all track-photon pairs with the Čerenkov angle smaller than 70 mrad. The Čerenkov angles of the pairs are stored in a list, together with a probability (weight, $`w_i`$) that the photon from a pair is emitted by the track of that pair. Initially the weights are set to $`1/N_{track}`$, where $`N_{track}`$ is the number of tracks corresponding to a specific photon hit. In the next step the expectation-maximization algorithm is applied to the list of pairs. The result of the algorithm is a set of new values of weights $`w_i`$. In this way the track-photon pair which is more likely to be the right one, gets a larger weight.

In the last step the extended likelihood probabilities are calculated for 6 possible hypotheses for the identity of a track: electron, muon, pion, kaon, proton and other. The last one represents the case, where the distribution of Čerenkov angles of photons for a given track is consistent with the background. The resulting likelihood probabilities have the range of values between 0 and 1 and are normalized so their sum equals to one. We note that the background level is determined on an event-by-event basis for each individual track.

The selection of tracks belonging to a particular particle type is made by applying a cut on the appropriate likelihood. For easier use, three levels of selection are defined: soft, medium and hard. The levels are set to 0.05 (soft), 0.50 (medium) and 0.95 (hard) for pion selection and 0.05 (soft), 0.30 (medium) and 0.95 (hard) for kaon or proton selection.

To determine the efficiencies and misidentification probabilities from the measured data, an a priori knowledge of particle types is required. The following reconstructed decays were used: $`(1)`$ $`K_S^0\pi ^+\pi ^{}`$, as a source of pions, $`(2)`$ $`\mathrm{\Lambda }p\pi ^{}`$ and $`\overline{\mathrm{\Lambda }}\overline{p}\pi ^+`$ as a source of protons, antiprotons and pions, and $`(3)`$ $`\varphi (1020)K^+K^{}`$, as a source of kaons. For the first two decays a very clear signal in the invariant mass plot is obtained by cutting on the secondary vertex distance and removing the reflections ($`\mathrm{\Lambda }`$, $`\overline{\mathrm{\Lambda }}`$ in $`K_S^0`$, and similarly for $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$). The $`\varphi `$ decay products come from the primary vertex together with roughly 10 other particles. The invariant mass plot shows a huge combinatorial background. To reduce the background, one of the particles was used to tag the decay by identifying it as a kaon. The other one was used for efficiency and misidentification evaluation. In all cases, the number of particles of a given type surviving the selection criteria, and the number of all particles of the same type in a selected momentum bin were obtained by fitting a Gaussian plus a linear function to the invariant mass plots.

The results are presented in Fig. 19. At medium selection criteria, kaon identification efficiency for momenta between 10 GeV/c and 60 GeV/c lies in a range from 60% to 95%, with typically about 5% pion mis-identification probability. The corresponding proton identification efficiency for momenta from 20 GeV/c to 60 GeV/c lies in a range between 60% and 80% with less than 5% pion mis-identification probability. In the region below kaon and proton Čerenkov threshold, the efficiencies are 70% for kaons (from 5 GeV/c to 10 GeV/c), and from 60% to 80% for protons (momentum region between 5 GeV/c and 20 GeV/c) with less than 20% pion mis-identification. Below 5 GeV/c the performance is degraded due to track direction uncertainty caused by multiple Coulomb scattering, and a small number of Čerenkov photons for pion tracks.

Examples of performance and impact on physics analyses are shown in Figs. 20-22. From Figs. 20 and 21, it is seen that the inclusion of the requirement that both tracks are identified improves the signal to background ratio in the $`\varphi K^+K^{}`$ and $`\mathrm{\Lambda }\pi ^{}p`$ decays. By using the particle identification as provided by the RICH counter, a measurement of hadron fractions was performed . The results shown in Fig. 22, agree with Monte Carlo predictions.

## 7 RICH as a Tracking Device

The RICH can also serve as a tracking device . This turned out to be very useful in the initial phase of the HERA-B experiment when the main tracking system was not available. From the centers of the rings found on the photon detector by the stand-alone ring search algorithm, track directions were deduced. In the non-bending plane they were matched with track candidates from the silicon vertex detector. From the deflection in the bending plane the momentum of the particle could then be determined. With both the Čerenkov angle and momentum deduced from the data, an excellent particle identification was possible even without the main tracking system, as illustrated in Fig. 23.

## 8 Summary

The HERA-B ring imaging Čerenkov detector behaves as expected in all respects. The detector parameters, 33 detected photons for $`\beta =1`$ particle, $`N_0=42`$ photons/cm and single photon resolutions of $`\sigma `$ = 0.7 mrad (1.0 mrad) for regions with finer (coarser) granularity, are in very good agreement with design values. Estimates of the basic parameters obtained from measurements of physical properties of RICH components agree well with the direct measurements. Therefore, we conclude that the functioning of the counter is well understood. Results on particle identification of highly populated events are also satisfactory and meet the requirements. Finally, it is worth adding that during the 4 years the detector has been in operation, not the least degradation of performance has been observed.

## 9 Acknowledgments

We want to express our thanks to the DESY laboratory for the strong support since the conception of the HERA-B experiment. We are especially grateful to the following persons and their groups: M. Bosteels (CERN), J. Eržen (IJS Ljubljana), K. Ludwig (DESY), E. Michel (DESY), H.B. Peters (DESY), P. Pietsch (DESY), M. Pohl (DESY Zeuthen), J. Spengler (MPI Heidelberg), and T. Stoye (DESY). We acknowledge the support from the U.S. Department of Energy (DOE), Ministry of Education, Science and Sports of Slovenia, Fundação para a Ciência e Tecnologia of Portugal, Spanish CICYT under the contract AEN99-0483, Internationales Büro, Jülich, and the Texas Advanced Research Program.


format: text

text:
# Imaging Studies of Characteristics for a Slab of a Lossy Left-handed Material

## I Introduction

Materials that have simultaneously a negative electric permittivity and a negative magnetic permeability have attracted much attention recently r1 ; r2 ; r3 ; r4 . These materials have extraordinary properties, such as the negative refraction, inverse Doppler shift and backward-directed Cherenkov radiation cone, etc veselago . The electric field, the magnetic field and the wave vector of an electromagnetic wave propagating in such a material obey the left-hand rule and thus it is called the left-handed material (LHM). Due to its effect of the negative refraction, a LHM planar slab can act as a lens, and focusing waves from a point source veselago . Recently, Pendry predicted that a LHM slab can amplify evanescent waves and thus behaves like a perfect lens pendry . It is well-known that a conventional imaging system is limited in resolution by the wavelength of the used light, since the information finer than the wavelength is carried by evanescent waves, which are drastically decayed before they reach the image plane. A LHM slab appears attractive in provideing an effective approach to achieve a subwavelength imaging. However, the LHM that Pendry considered is quite ideal and has no loss. Garcia et al garcia  pointed out later that an actual LHM should be lossy and the LHM loss may completely suppress the amplification of evanescent waves, and consequently a LHM slab can not act as a perfect lens. It is interesting to find out answers to the following questions: (1) What is the characteristic of a realistic LHM imaging system? (2) How does the resolution of the system depend on the LHM loss and other slab parameters? (3) Is it possible to use a realistic LHM imaging system to make a subwavelength imaging? If yes, then how to realize it? In the present paper, we will study these problems theoretically.

## II  Transfer function of a LHM imaging system

Consider a LHM slab imaging system with co-ordinates as shown in Fig. 1. We assume that the system is invariant along the y-axis direction and only the one-dimensional imaging will be considered here for simplicity. In the imaging system, the space can be divided into three regions: region 1 ($`z<0`$), where the objective plane is located; region 2 (the LHM slab); region 3 ($`zd`$), where the image plane is located. The electric permittivity and magnetic permeability in the three regions are denoted by $`ϵ_l`$ and $`\mu _l`$ ($`l=1,2,3`$), respectively. Here $`ϵ_1=ϵ_30`$, $`\mu _1=\mu _30`$, $`Reϵ_2<0`$, and $`Re\mu _2<0`$. This imaging system has the properties of being linear and invariant under translation along the $`x`$ direction, and thus can be described by a transfer function ($``$) or impulse response ($`h`$) (the latter is just the Fourier transform of the former).

In order to calculate the transfer function of the LHM imaging system, we consider a plane wave propagating along the z direction from the objective plane $`z=a`$ ($`a>0`$). The plane wave will refract and reflect simultaneously at the interfaces $`z=0`$ and $`z=d`$ before it reaches the image plane $`z=d+b`$ (at a distance of b away from the LHM slab). The electromagnetic fields can be decomposed into the $`E`$\- and $`H`$-polarization modes. Let’s first consider the $`E`$-polarization case, where the electric field has only $`y`$ component (i.e., $`\stackrel{}{E}=E\widehat{y}`$), and the magnetic field has the form of $`\stackrel{}{H}=H_x\widehat{x}+H_z\widehat{z}`$. All the field components are assumed to have the time dependence of form $`\mathrm{exp}(i\omega t)`$. The electric field $`E`$ of the wave in the three regions can be written in the following form of

$`E_1(x,z)`$ $`=`$ $`\mathrm{exp}(ik_xx)\left\{A_1\mathrm{exp}(ik_z(z+a))+B_1\mathrm{exp}(ik_zz)\right\},`$ ()
$`E_2(x,z)`$ $`=`$ $`\mathrm{exp}(ik_xx)\left\{A_2\mathrm{exp}(iqz)+B_2\mathrm{exp}(iqz)\right\},`$ ()
$`E_3(x,z)`$ $`=`$ $`A_3\mathrm{exp}(ik_xx)\mathrm{exp}(ik_z(zd)),`$ ()
where $`k_x`$ is the $`x`$-component of the wave vector, $`k_z`$ is the z-component in regions 1 and 3, and $`q`$ is the $`z`$-component of the wave vector in region 2. Obviously, one has $`k_z=\sqrt{n_1^2k_0^2k_x^2}`$ for propagating waves (i.e., $`k_x<n_1k_0`$, where $`n_1=\sqrt{ϵ_1\mu _1}`$), and $`k_z=i\sqrt{k_x^2n_1^2k_0^2}`$ for evanescent waves (i.e., $`k_xn_1k_0`$). Similarly, $`q=\sqrt{ϵ_2\mu _2k_0^2k_x^2}`$ for propagating waves and $`q=i\sqrt{k_x^2ϵ_2\mu _2k_0^2}`$ for evanescent waves. The coefficients for the forward and backward (concerning the energy flow) propagating components are denoted by $`A`$ and $`B`$, respectively. Note that the energy flow for the term $`A_2\mathrm{exp}(iqz)`$ in Eq.(1.2) is along the $`+z`$ direction since the refraction index of the LHM slab is negative.

The transfer function of the imaging system is defined as $`(k_x)=E_3(x,d+b)/A_1\mathrm{exp}(ik_xx)`$, i.e., $`(k_x)=(A_3/A_1)\mathrm{exp}(ik_zb)`$. From Maxwell’s equations, one has $`H_x=(i/\omega \mu _0\mu )(E/z)`$. From the continuity condition of the field components of $`E`$ and $`H_x`$ at the interfaces $`z=0`$ and $`z=d`$, one can derive formulas for $`A_3/A_1`$ and $`(k_x)`$. In the present paper we can write the transfer function $`(k_x)`$ in the following product form (with each term having a clear physical interpretation)

$`(k_x)`$ $`=`$ $`t_1Pt_2\overline{P},`$ ()
where

$`t_1`$ $`=`$ $`{\displaystyle \frac{2}{(1\xi )+(1+\xi )P^2}}{\displaystyle \frac{1}{t_2}},`$ ()
$`t_2`$ $`=`$ $`2(1{\displaystyle \frac{\mu _2k_z}{\mu _1q}})^1,`$ ()
$`P`$ $`=`$ $`\mathrm{exp}(iqd),`$ ()
$`\overline{P}`$ $`=`$ $`\mathrm{exp}\{ik_z(a+b)\},`$ ()
$`\xi `$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\mu _1q}{\mu _2k_z}}+{\displaystyle \frac{\mu _2k_z}{\mu _1q}}\right).`$ ()
Here $`t_1`$ represents the effective transmission coefficient for a plane wave at the interface $`z=0`$ (i.e., $`t_1=A_2/(A_1exp(ik_za))`$), while $`t_2`$ is the simple transmission coefficient at $`z=d`$ (i.e., $`t_2=A_3/(A_2exp(iqd))`$). $`\xi `$ is related to the impedance matching of the LHM with the surrounding medium ($`\xi =1`$ when the impedance is matched). The second term $`(1+\xi )P^2`$ in the denominator of expression (3.1) for $`t_1`$ is attributed to the effect of the reflection of the wave at the interface $`z=d`$. In fact, this term ($`(1+\xi )P^2`$) vanishes if the LHM slab is completely matched in the impedance with the medium in region 3 (i.e. $`ϵ_2=ϵ_3`$, $`\mu _2=\mu _3`$) regardless what the value of $`ϵ_1`$ or $`\mu _1`$ is. $`P`$ is the phase change factor or amplitude amplification factor as the wave passes the LHM slab, and $`\overline{P}`$ is the total phase change factor or amplitude decaying factor as the wave goes from the objective plane to the left surface of the LHM slab and from the right surface of the LHM slab to the image plane.

For an ideal LHM imaging system (i.e., $`ϵ_2=ϵ_1`$, $`\mu _2=\mu _1`$ and $`b=da`$), one has $`\xi =1`$, $`t_1=1`$, $`t_2=1`$, and $`(k_x)=1`$ for all $`k_x`$, this means the imaging system is perfect.

## III Influence of the loss in the LHM

It is known that an actual LHM should be lossy (i.e. $`ϵ_2`$ or $`\mu _2`$ has an imaginary part) since any LHM must be dispersive smith . Here we analyze the characteristics of the transfer function for an imaging system formed by a lossy LHM slab. Consider first the case when the loss comes only from the imaginary part of $`ϵ_2`$, i.e.,

$$ϵ_2=ϵ_1(1i\delta ),\mu _2=\mu _1,$$

where we assume $`0<\delta 1`$. As $``$ is an even function of $`k_x`$, we only needs to consider the interval $`[0,\mathrm{})`$ for $`k_x`$. For the convenience of analysis, we further divide this interval into two regions, namely, $`k_x[0,n_1k_0)`$ and $`k_x[n_1k_0,+\mathrm{})`$.

1. The case when $`k_x<n_1k_0`$. This case corresponds to propagating waves with $`k_z=\sqrt{n_1^2k_0^2k_x^2}`$ being a real number. Then one has $`q=\sqrt{ϵ_2\mu _2k_0^2k_x^2}=\sqrt{k_z^2in_1^2k_0^2\delta }`$.

(i) Consider the special situation when the plane waves are propagating nearly along the z axis, i.e., $`k_x\sqrt{1\delta }n_1k_0`$ (or equivalently, $`k_z\sqrt{\delta }n_1k_0`$). In this case, one has $`\xi 1+O(\delta ^2)`$, which implies that the media are approximately matched and $`t_11/t_21`$. The transfer function can then be approximated by

$`\left(k_x\right)\mathrm{exp}\left({\displaystyle \frac{n_1k_0}{2\sqrt{n_1^2k_0^2k_x^2}}}D\delta \right),`$ ()
where $`D=n_1k_0d`$. At the point $`k_x=0`$ (corresponding to the plane wave normally incident on the LHM slab), one has $`(0)=\mathrm{exp}(n_1k_0d\delta /2)`$. This indicates that the small loss of the LHM makes $`H(0)`$ a bit below $`1`$.

(ii) Then consider the other situation when $`k_xn_1k_0`$ (i.e., $`k_z0`$). In this situation, one finds that $`k_z\sqrt{\delta }n_1k_0`$ and $`\xi =\sqrt{i\delta }n_1k_0/(2k_z)`$ differs from $`1`$ completely. Therefore, a small loss of the LHM causes a serious mismatch problem for these waves. The transfer function can be approximated by

$`(k_x){\displaystyle \frac{2\sqrt{n_1^2k_0^2k_x^2}}{n_1k_0D\delta }}\mathrm{exp}\left\{D\sqrt{{\displaystyle \frac{\delta }{2}}}(1+i)\right\},`$ ()
Obviously, $`(k_x)0`$ when $`k_xn_1k_0`$. Therefore, the existence of a small loss in the LHM changes the value of the $``$ function completely (dropping to 0 from 1; see Fig. 2(c) below) in the neighborhood of $`k_x=n_1k_0`$.

2. The case when $`k_xn_1k_0`$. This case corresponds to evanescent waves with $`k_z=i\sqrt{k_x^2n_1^2k_0^2}`$. Then one has $`q=i\sqrt{k_x^2ϵ_2\mu _2k_0^2}=i\sqrt{|k_z|^2+i\delta n_1^2k_0^2}`$.

(i) Consider first the situation when $`k_xn_1k_0`$ (i.e., $`k_z0`$). In this situation, one finds that $`|k_z|\sqrt{\delta }n_1k_0`$ and $`\xi \sqrt{i\delta }n_1k_0/(2|k_z|)`$ differs very much from $`1`$. Thus, one sees that a small loss of the LHM causes a serious mismatch problem for these evanescent waves. One can then approximate the transfer function as

$`(k_x)i{\displaystyle \frac{2\sqrt{k_x^2n_1k_0^2}}{n_1k_0D\delta }}\mathrm{exp}\left\{D\sqrt{{\displaystyle \frac{\delta }{2}}}(1+i)\right\}.`$ ()
Obviously, $`0`$ when $`k_xn_1k_0`$. We find that the real part ($`_r`$) of $``$ is far smaller than the imaginary part ($`_i`$) in this case. $`_i`$ decreases as $`k_x`$ decreases to $`n_1k_0`$. On the other hand, $`_i`$ should approach zero when $`k_xn_1k_0`$ garcia . Thus, there exists a peak of $`_i`$ on the right-hand side of the point $`k_x=n_1k_0`$ (see Fig. 2(d) below).

(ii) Then consider the other situation when $`k_x\sqrt{1+\delta }n_1k_0`$ (or equivalently, $`|k_z|\sqrt{\delta }n_1k_0`$). These evanescent waves are very important for a subwavelength imaging. For this case our analysis shows that $`\xi =1+(n_1k_0/k_z)^4\delta ^2/8`$, which indicates that $`\delta `$ causes only a second-order perturbation to $`\xi `$ (the LHM slab seems still well-matched with the surrounding media for these evanescent waves). However, $`|P|exp(\sqrt{k_x^2n_1^2k_0^2}d)`$ increases exponentially as $`k_x`$ increases. Thus, unlike the case of propagating waves, the second term $`(1+\xi )P^2`$ in the denominator of the expression (3.1) for $`t_1`$ can not be neglected. In other words, a tiny $`\delta `$ may cause a significant change of the effective transmission coefficient $`t_1`$ for the evanescent wave at the interface $`z=0`$ (note that $`\delta `$ causes only a small perturbation to $`t_2`$ and $`P`$), and consequently leads to a quick decrease of $`||`$. For a slab of usual medium (with or without loss), the effective transmission coefficient $`t_1`$ must be close to 1 if the impedance is almost matched (regardless the value of $`k_x`$). However, for a LHM slab $`t_1`$ can be very small (consequently the total transmission coefficient of the LHM slab can be very small) even the the impedance is almost matched. The main reason is that the evanescent wave ”reflected” at the back interface of the LHM slab may be effectively amplified (before it comes out of the slab) and thus significantly suppress the transmission of the evanescent wave according to Eq. (3.1). Therefore, the conventional concept of impedance match is insufficient to describe the medium match for evanescent waves for a LHM slab. We introduce a critical value ($`k_x^c`$) for $`k_x`$, at which $`Re1/2`$. Apparently, $`k_x^c`$ is closely related to the resolution of the imaging system. After a perturbation analysis for Eqs. (2)-(3.5), one can conclude that the critical value $`k_x^c`$ can be approximately determined from the condition $`Re[(\xi +1)P^2]=2`$ (physically this condition can be understood as a balance condition for the two terms in the denominator of expression (3.1) for $`t_1`$). When $`k_xk_x^c`$, the term $`(\xi +1)P^2`$ becomes dominant in the denominator of expression (3.1) for $`t_1`$, and $`t_1`$ decreases rapidly to 0 (consequently $`0`$) as $`k_x`$ increases. The transfer function can then be approximated as

$`(k_x)16\left({\displaystyle \frac{k_x^2n_1^2k_0^2}{n_1^2k_0^2\delta }}\right)^2\mathrm{exp}(2\sqrt{k_x^2n_1^2k_0^2}d),k_xk_x^c.`$ ()
Figs. 2(a) and 2(b) show the real and imaginary parts of $`(k_x)`$ (calculated with Eq. (2)), respectively, where we have chosen $`d=\lambda `$ ($`\lambda =2\pi /n_1k_0`$) and $`\delta =10^4`$. Figs. 2(c) and 2(d) give the corresponding enlarged views of $`_r(k_x)`$ and $`_i(k_x)`$ in the neighborhood of $`n_1k_0`$. From Figs. 2(b)and 2(d) one sees that $`_i`$ has a peak around $`k_x=n_1k_0`$ (a bit larger than $`n_1k_0`$) as discussed before. Nevertheless, $`_i`$ can be approximated to zero everywhere (except in a very narrow region around $`n_1k_0`$; this non-zero region is so narrow that it can be neglected in an imaging system for the image in the physical space). Our numerical analysis shows that $`k_x^c2.3n_1k_0`$ for this example. Fig. 2(a) shows that $``$ decreases to zero rapidly as $`k_x`$ increases in the region $`k_xk_x^c`$. From these figures one sees that the characteristics of the transfer function $``$ (calculated from Eq. (2)) agree very well with our earlier analysis.

The value of $`k_x^c`$ is very important for an imaging system since $`k_x^c`$ directly relates the smallest scale in which information can be correctly transferred by the system. To see clearly what influences this parameter, we wish to derive an analytic formula for $`k_x^c`$. From $`Re[(\xi +1)P^2]=2`$, one obtains $`\gamma e^\gamma =p`$ (i.e., $`\mathrm{ln}\gamma \gamma =\mathrm{ln}p`$), where $`\gamma =D\sqrt{(k_x^c)^2n_1^2k_0^2}/(2n_1k_0)`$ and $`p=D\sqrt{\delta }/4`$. We assume that loss parameter $`\delta `$ is so small that $`p1`$. Then one has $`e^\gamma \gamma >1`$, i.e., $`\gamma \mathrm{ln}\gamma >0`$. Thus it follows that $`\gamma =\mathrm{ln}p+\mathrm{ln}\gamma \mathrm{ln}p`$. This approximation can be improved as $`\gamma \mathrm{ln}p+\mathrm{ln}\mathrm{ln}(1/p)`$. One can finally obtains the following analytic formula

$`k_x^c`$ $``$ $`\left\{\left[{\displaystyle \frac{2}{D}}\mathrm{ln}\left({\displaystyle \frac{4}{D\sqrt{\delta }}}\right)+{\displaystyle \frac{2}{D}}\mathrm{ln}\mathrm{ln}\left({\displaystyle \frac{4}{D\sqrt{\delta }}}\right)\right]^2+1\right\}^{\frac{1}{2}}\left(n_1k_0\right).`$ ()
Note that $`D=n_1k_0d`$.

Eq. (8) shows that for a given loss parameter$`\delta `$, $`k_x^c`$ decreases as the thickness of the LHM slab increases (actually $`k_x^c`$ is almost inverse proportional to $`d`$ as one can see from $`\gamma \mathrm{ln}p`$). On the other hand, $`k_x^c`$ is almost linearly proportional to $`\mathrm{ln}\delta `$ (since $`\left(2/D\right)\mathrm{ln}\left[4/\left(D\sqrt{\delta }\right)\right]=\mathrm{ln}\left(\delta \right)/D+2\mathrm{ln}\left(16/D\right)/D`$ and the first term is dominant over the second term). Thus, $`k_x^c`$ does not increase noticeably even if $`\delta `$ decrease by one order of magnitude. To realize the subwavelength imaging, e.g., $`k_x10k_0`$, it is required that $`\delta <10^{52}`$ when $`d=2\lambda `$, and $`\delta <2\times 10^{25}`$ when $`d=\lambda `$. It is doubtful whether the loss of a realistic LHM could be reduced to such a level. Therefore, the thickness of the LHM slab should be smaller than the wavelength in order to realize the subwavelength imaging for such a system.

From Eq. (2)-(3.5) (or Eq. (8)) one can find that the characteristics of the transfer function (as well as $`k_x^c/(n_1k_0)`$) depend only on $`\delta `$ and $`d`$/$`\lambda `$. Fig. 3(a) gives the contour plot of $`k_x^c/(n_1k_0)`$ (calculated with Eq. (2)) as a function of $`\delta `$ and $`d`$/$`\lambda `$. From this figure one sees that $`k_x^c`$ is very sensitive to the thickness $`d`$. Our numerical calculation shows that the requirement for $`k_x^c=10n_1k_0`$ (for a subwavelength imaging) is $`\delta =2\times 10^{25}`$ when $`d=\lambda `$, $`\delta =9\times 10^{11}`$ when $`d=\lambda /2`$, and $`\delta =6\times 10^5`$ when $`d=\lambda /4`$. Fig. 3(b) shows $`k_x^c`$ obtained from Eq.(2) and the analytic expression (8), respectively, as the loss parameter varies. One sees that they agree very well.

Consider now the impulse response for a lossy LHM imaging system. The impulse response $`h(x)`$ is defined as the field distribution on the image plane produced by a point source on the objective plane. Obviously, it is just the Fourier transform of the transfer function $`(k_x)`$. Our analysis shows that $`(k_x)`$ is approximately a step function, i.e., $`(k_x)1`$ for $`|k_x|<k_x^c`$ (except the abrupt change in the narrow neighborhood of $`k_x=n_1k_0`$, which is narrow enough to be neglected in calculating the impulse response by the Fourier transform), while $`(k_x)0`$ elsewhere. Thus, we have approximately $`h(x)=2k_x^csin(k_x^cx)/(k_x^cx)`$. According to the definition of the resolution (denoted by $`\mathrm{\Delta }`$) for an imaging system, we find that $`\mathrm{\Delta }=\pi /k_x^c`$. Figs. 4(a) and 4(b) show the real and imaginary parts of the impulse response (normalized by $`h(0)`$) of a LHM imaging system, respectively, with $`\delta =10^4`$ and $`d=\lambda `$. One sees that the imaginary part of the impulse response is very small and can be neglected as compared with the real part. From Fig. 4(a) one sees that the real part of $`h`$ calculated from the Fourier transform of Eq.(2) (the solid line) agrees very well with that from our analytic expression (the dotted line) in the region of $`|x|<2\pi /k_x^c`$ (note that our analytic expression gives a real value of $`h`$; the imaginary part of $`h`$ is always very small, as show in Fig. 4(b)). This indicates that the resolution expression $`\mathrm{\Delta }=\pi /k_x^c`$ is quite reliable. Thus, one can use the following explicit resolution for a lossy LHM imaging system

$`\mathrm{\Delta }\left\{\left[{\displaystyle \frac{2}{D}}\mathrm{ln}\left({\displaystyle \frac{4}{D\sqrt{\delta }}}\right)+{\displaystyle \frac{2}{D}}\mathrm{ln}\mathrm{ln}\left({\displaystyle \frac{4}{D\sqrt{\delta }}}\right)\right]^2+1\right\}^{\frac{1}{2}}\left({\displaystyle \frac{\pi }{n_1k_0}}\right).`$ ()
In the above analysis, we have assumed that the loss of the LHM comes from the imaginary part of $`ϵ_2`$. We can also discuss the case in which $`\mu _2`$ gives rise to the loss, i.e. $`ϵ_2=ϵ_1`$, $`\mu _2=\mu _1(1i\delta )`$. Our analysis and numerical calculation have shown that the characteristics of the transfer function for this case is similar to the previous case that we have discussed. In this case one has $`k_x^c\left(n_1k_0\right)\sqrt{\left[\left(1/D\right)\mathrm{ln}\left(2/\delta \right)\right]^2+1}`$ and $`\mathrm{\Delta }\left(\pi /n_1k_0\right)/\sqrt{\left[\left(1/D\right)\mathrm{ln}\left(2/\delta \right)\right]^2+1}`$. For the same $`d`$ and $`\delta `$, $`k_x^c`$ for this case is a bit smaller than that for the previous case (see Fig.5).

In the above analysis the electromagnetic wave is assumed to be $`E`$ polarized. For the $`H`$polarized case, one can obtain the transfer function from Eq.(2) directly by interchanging $`ϵ`$ and $`\mu `$ due to the symmetry of Maxwell’s equations. Thus, similar analytic formulas for $`k_x^c`$ and $`\mathrm{\Delta }`$ can be found.

## IV Conclusion and discussion

In this paper we have analyzed the characteristics of the transfer function for an imaging system formed by a lossy LHM slab. The small loss in the LHM may completely change the transfer function in the following aspects: (1) the transfer function drops to zero abruptly near the turning point $`k_x=n_1k_0`$ for the propagating and evanescent waves; (2) $``$ decreases to zero rapidly when $`k_x`$ increases from a critical $`k_x^c`$ (at which $`=1/2`$). Our analysis shows that $`k_x^c`$ is almost linearly proportional to the logarithm of the loss parameter $`\delta `$ (related to the imaginary part of $`ϵ_2`$ or $`\mu _2`$) of the LHM. Thus, $`k_x^c`$ does not increase noticeably even if the loss decreases by one order of magnitude. Analytical formulas for $`k_x^c`$ and the resolution of the imaging system have been derived. It has been shown that the thickness of the LHM slab should be smaller than the wavelength in order to realize the subwavelength imaging for such a system.

Our analysis has also shown that a tiny loss of the LHM may suppress the transmission of evanescent waves through the LHM slab due to the amplification of the evanescent waves in the LHM slab. Such a suppression can be effectively relaxed by reducing the thickness of the LHM slab. We have also tried to eliminate or reduce the destructive influence of the LHM loss to the imaging resolution through manupilating the real and imaginary parts of $`\mu _2`$ and $`ϵ_2`$. However, we found this is not possible (see the appendix).

###### Acknowledgements.

The partial support of Natural Science Foundation of China under a project (60277018) and a key project (90101024) is gratefully acknowledged.

## Appendix: Possibility of reducing the destructive influence of the LHM loss of $`\mu _2`$ and $`ϵ_2`$

In this appendix, we study the possibility of eliminating or reducing the destructive influence of the LHM loss to the image through manupilating the real and imaginary parts of $`\mu _2`$ and $`ϵ_2`$. Consider the case when the electric permittivity and magnetic permeability of the LHM slab have the following forms

$`\epsilon _2`$ $`=`$ $`\epsilon _1(1\delta ^{^{}}),`$
$`\mu _2`$ $`=`$ $`\mu _1(1\delta ^{^{\prime \prime }}),`$
where $`\delta ^{^{}}=\delta _r^{^{}}+i\delta _i^{^{}}`$, $`\delta ^{^{\prime \prime }}=\delta _r^{^{\prime \prime }}+i\delta _i^{^{\prime \prime }}`$ ($`\delta _r^{^{}}`$, $`\delta _i^{^{}}`$, $`\delta _r^{^{\prime \prime }}`$, and $`\delta _i^{^{\prime \prime }}`$ are all real). We assume that $`\delta _i^{^{}},\delta _i^{^{\prime \prime }}>0`$, and $`|\delta ^{^{}}|,|\delta ^{^{\prime \prime }}|1`$. In the following analysis for the LHM slab, our attention will be focused on the evanescent waves for the subwavelength imaging.

From Eq. (3.5), one has

$`\xi =1+{\displaystyle \frac{1}{2}}\left(\delta ^{^{\prime \prime }}+{\displaystyle \frac{\overline{\delta }}{\beta ^2}}\right)^2+O(\overline{\delta }^3),`$
where $`\overline{\delta }=(\delta ^{^{}}+\delta ^{^{\prime \prime }})/2`$, $`\beta =\sqrt{k_x^2n_1^2k_0^2}/(n_1k_0)`$ ($`k_x>n_1k_0`$ for the evanescent waves). Then the transfer function becomes

$$\frac{1}{1+0.25(\delta ^{^{\prime \prime }}+\overline{\delta }/\beta )^2\mathrm{exp}(2D\beta )}.$$

Therefore, the requirement for reducing significantly the influence of the LHM loss is

$$\delta ^{^{\prime \prime }}=0,\overline{\delta }=0.$$

Obviously the above requirement is equivalent to that of $`\delta ^{^{}}=0`$ and $`\delta ^{^{\prime \prime }}=0`$, which in fact corresponds to the case of an ideal LHM without any loss. Thus, it is impossible to effectively reduce the destructive influence of the LHM loss to the subwavelength imaging by manupilating the real and imaginary parts of $`\mu _2`$ and $`ϵ_2`$.


text:
# 1 Normalized correlation of intertrade time intervals in tick time. Solid line shows a powerlike fit of the form 0.15/𝑑^0.29. Upper and lower lines show the standard deviation of the measured correlation characterizing its variation on the day-to-day basis.

## Appendix A Derivation of Eq. (3)

In this Appendix we sketch, following , the derivation of Eq. (3). The probability of having $`N`$ trades in the interval $`[0,T]`$ is conviniently written as a convolution of a probability of having $`N`$ trades in the interval $`[0,t]`$ (the $`N`$-th trade falling into the infinitesimal vicinity of $`t`$) $`w_N(t)`$ and a probability of having no trades in the remaining interval $`[t,T]`$ $`w_0(Tt)`$:

$$p_N(T)=_0^T𝑑tw_0(Tt)w_N(t)$$
(A.1)
so that the Laplace transform of $`p_N(T)`$ reads

$$\stackrel{~}{p}_N(\lambda )=𝑑T\mathrm{e}^{\lambda T}p_N(T)=\stackrel{~}{w}_0(\lambda )\stackrel{~}{w}_N(\lambda )=\frac{1\stackrel{~}{\psi }(\lambda )}{\lambda }\stackrel{~}{\psi }(\lambda )^{N1},$$
(A.2)
where we took into account that in the considered unit memory depth case all probabilities in question are fully determined by the intertrade (waiting) time distribution $`\psi (\tau )`$, so that 555Here we assume that there is always a first trade at the beginning of the time interval in question $`\stackrel{~}{w}_N(\lambda )=\stackrel{~}{w}(\lambda )^{N1}`$ and, for $`w_0`$ we have, recalling that $`w_0(\tau )=_\tau ^{\mathrm{}}𝑑\tau ^{}\psi (\tau ^{})`$, $`\stackrel{~}{w}_0(\lambda )=(1\stackrel{~}{\psi }(\lambda ))/\lambda `$.

To compute the moments of trade’s multiplicity distribution it is convenient to introduce a generating function

$$G(\theta |T)=\underset{k=1}{\overset{\mathrm{}}{}}\mathrm{e}^{k\theta }p_k(T)=\mathrm{\hspace{0.17em}1}+\underset{k=1}{\overset{\mathrm{}}{}}\frac{\theta ^k}{k!}N^k$$
(A.3)
so that the moments are given by its derivatives at the origin

$$N^k_T=\left(\frac{^kG}{\theta ^k}\right)_{\theta =0}.$$
(A.4)
Combining Eqs (A.2,A.3) one gets the following expression for the generating function:

$$G(\theta |T)=\frac{1}{2\pi i}_{ci\mathrm{}}^{c+i\mathrm{}}𝑑\lambda \mathrm{e}^{\lambda T}\frac{1}{\lambda }\frac{1\psi (\lambda )}{\mathrm{e}^\theta +\psi (\lambda )}$$
(A.5)
Equation (3) follows now by substituting the expansion of $`\psi (\lambda )`$ in the moments of $`\psi (\tau )`$

$$\psi (\lambda )=_0^{\mathrm{}}𝑑\tau \mathrm{e}^{\lambda \tau }\psi (\tau )=\mathrm{\hspace{0.17em}1}\psi _1\lambda +\frac{1}{2}\psi _2\lambda ^2+\mathrm{}$$
(A.6)
into Eq. (A.5) and computing, through Eq. (A.4), the leading contributions to $`N_T`$ and $`N_T^2`$ in the large $`T`$ limit

$`N_T`$ $``$ $`{\displaystyle \frac{T}{\psi _1}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\psi _2}{\psi _1^2}}`$ (A.7)
$`N^2_T`$ $``$ $`{\displaystyle \frac{T^2}{\psi _1^2}}+2T{\displaystyle \frac{\psi _2}{\psi _1^3}}{\displaystyle \frac{T}{\psi _1}}`$ (A.8)


text:
# External shock model for the large-scale, relativistic X-ray jets from the microquasar XTE J1550-564

## 1. Introduction

Black-hole X-ray binaries with relativistic jets resemble, on a much smaller scale, many of the phenomena seen in quasars and are therefore called microquasars (Mirabel & Radriguez 1999). Apparently superluminal radio jets are observed from at least the two best known microquasars GRS 1915+105 ( Mirabel & Rodriguez 1994; Fender et al. 1999) and GRO J1655-40 (Tingay et al. 1995; Hijellming & Rupen 1995), with the actual jet velocities greater than 0.9c. Their motions are consistent with being purely ballistic, up to a projected separation of 0.08pc from the compact object on a maximum time scale of 4 months after the ejection for GRS 1915+105. With the help of Chandra Observatory, Corbel et al. (2002) recently discovered large-scale (with projected separation more than a half pc from the compact object), relativistically moving and decelerating radio and X-ray emitting jets from the microquasar XTE J1550-564. It is the first time that an X-ray jet proper-motion measurements spanning several years are obtained for accretion-powered Galactic sources. Hence, XTE J1550-564 provides a good opportunity to study the dynamical evolution of relativistic jets.

In this paper, we propose that the dynamical evolution and radiation of the large-scale X-ray jets from XTE J1550-564 can be understood as the interaction between the jet and the surrounding interstellar medium (ISM), quite similar to the external shock model for afterglows of gamma-ray bursts (GRBs) (Rees & Mészáros 1992; Mészáros & Rees 1997). The main differences lie in the initial Lorentz factor and the kinetic energy of the jets. The microquasar jets are inferred to be at least mildly relativistic with initial Lorentz factor $`\mathrm{\Gamma }_025`$ (Mirabel & Radriguez 1999)111Recently, Fender (2003) argued that the two-sided jet proper motions observed from microquasars can only allow placing lower limits on the initial Lorentz factors. , in contrast to the ultra-relativistic jets ($`\mathrm{\Gamma }_0>100`$) in GRBs. When the relativistic ejecta from the microquasar is significantly decelerated by the ISM, a relativistic forward shock expands into the ISM and a reverse shock moves into and heats the original ejecta. The shocked ambient and ejecta materials are in pressure balance and separated by a contact discontinuity. The forward shock continuously heats fresh ISM and accelerates electrons, while the reverse shock operates only once and after that the shocked gas in the ejecta expands and cools adiabatically. As more and more ISM matter are swept-up, the ejecta and the shocked ISM — we shall call them ’jet’— should be decelerated more and more and finally transit to the non-relativistic motion phase. We find that this dynamic model can fit well the observed proper-motion evolution of the large scale eastern X-ray jet from XTE J1550-564.

In both shocks, the kinetic energy are converted to the internal energy. In the standard theory of GRB afterglows, it is assumed that shocked electrons and the magnetic field acquire constant fractions of the total shock energy, denoted by $`ϵ_e`$ and $`ϵ_B`$ for electrons and magnetic field respectively. The inferred values from a few afterglows are in the range: $`ϵ_e0.10.6`$ and $`ϵ_B10^60.1`$ ( e.g. Granot et al. 1999; Wijers & Galama 1999; Panaitescu & Kumar 2002). If we assume that the same processes also occur in the large-scale, decelerating jet we consider here, we are able to know the emission from the forward shock and reverse shock regions. We find that forward shock emission decays too slowly to be consistent with the X-ray emission from the eastern jet of XTE J1550-564, but the emission from the post-shocked adiabatically expanding ejecta can give a reasonable fit.

First, we give a brief review of the observations of the large scale jets from XTE J1550-564 in section 2. We present the dynamic model fit in section 3 and interpret the radiation in section 4. Finally, we give conclusions and discussions.

## 2. Observations of the large-scale jets from XTE J1550-564

The X-ray transient XTE J1550-564 was discovered by the All-Sky Monitor on aborad the Rossi X-ray Timing Explorer on 7 September 1998 (Smith 1998). Optical observations of the source in quiescence showed that the mass of the compact object is near $`10M_{}`$, indicating that the compact object is a black hole and revealed that the binary companion to be low-mass star (Hannikainen et al. 2001). The distance ($`d`$) to the source is constrained to be in the range 2.8-7.6 kpc with a favored value of 5.3 kpc (Orosz et al. 2002). Soon after the discovery of the source, an extremely strong X-ray flare was observed on 20 September 1998 (Sobczak et al. 2000; Homan et al. 2001), and radio jets with apparent superluminal velocities (the initial proper motion was greater than $`57\mathrm{mas}\mathrm{day}^1`$) was observed beginning 24 September 1998 (Hannikainen et al. 2001). During the 2002 X-ray outburst, radio observations were made with the Australia Telescope Compact Array (ATCA). The detection of the large-scale radio jet $`22\mathrm{arcsec}`$ to the west of XTE J1550-564 led to a re-analysis of the archival Chandra data and discovery of an X-ray jet to the east of XTE J1550-564 (Corbel et al. 2002; Tomsick et al. 2003). It is thought that both jets are connected with the 20 September 1998 ejection event, based on the detection of superluminal jets following an extremely large X-ray flare and the absence of any other X-ray flare of similar magnitude in continual X-ray monitoring from 1996 to 2002.

According to the archival Chandra data, the field of view of XTE J1550-564 was imaged by Chandra on 9 June, 21 August, and 11 September 2000. These observations present for the first time the proper-motion measurement for an X-ray jet from microquasars, which show direct evidence for gradual deceleration. Radio observations with ATCA show a decaying and moving radio source consistent with the position of the eastern X-ray jet. Remarkably, the overall radio flux detected on 1 June 2002 for the eastern jet is consistent with an extrapolation of the X-ray spectrum with a single power law of spectral index of $`\alpha 0.6`$ ($`F_\nu \nu ^\alpha `$) (Tomsick et al. 2003), which provides evidence for a synchrotron X-ray emission mechanism222 The broadband spectral energy distribution of the western jet around 11 March 2002, which is consistent with a single power law of spectral index of $`0.660\pm 0.005`$, strengthened the synchrotron radiation origin of the X-ray emission (Corbel et al. 2002) .. The field of view of XTE J1550-564 was further observed by Chandra on 11 March 2002. We summarize the observations for both the eastern and western jets in Table 1.

For the western jet, Kaaret et al. (2003) also examined the archival Chandra data on XTE J1550-564 at June, August and September 2000, but found no evidence for X-ray emission in any of the archival observations with an upper bound on the absorbed flux of $`1.1\times 10^{14}\mathrm{ergcm}^2\mathrm{s}^1`$. The facts that no X-ray emission detected from the western jet during 2000 and that the eastern source apparently moves faster than the western one are consistent with the interpretation in which the eastern jet is the approaching one and the western jet is the receding one. The brightening of the western jet at the late time argues against symmetric jet propagation and may reflect the non-uniformity in the ISM.

## 3. The dynamic model

As we only know the proper motion measurements and the light curve of the eastern jet, we focus on modelling the dynamical evolution and radiation of this jet.

We envision a beamed outflow with kinetic energy $`E_0`$ and Lorentz factor $`\mathrm{\Gamma }_0`$ ejected from the microquasar XTE J1550-564 expands conically with a half opening angle $`\theta _j`$ into the ambient medium with a constant number density $`n`$. The interaction between the relativistic ejecta and the surrounding medium is analogous to GRB external shock, but with quite different $`E_0`$ and $`\mathrm{\Gamma }_0`$. The supersonic motion of the ejecta should drive a blast wave propagating into the ISM. We shall assume that the radiation loss of the shock wave is a negligible fraction of the total energy, so the dynamic is adiabatic throughout. From the view of the energy conservation, the dynamic equation can be simplified as (see also Huang, Dai & Lu 1999)

$$(\mathrm{\Gamma }1)M_0c^2+\sigma (\mathrm{\Gamma }_{\mathrm{sh}}^21)m_{\mathrm{sw}}c^2=E_0,$$
(1)
where $`\mathrm{\Gamma }`$ and $`\mathrm{\Gamma }_{\mathrm{sh}}`$ are the Lorentz factors of the jet material and the shock front respectively, $`m_{\mathrm{sw}}=(4/3)\pi R^3m_pn(\theta _j^2/4)`$ is the mass of the swept-up ISM (where $`R`$ is the shock radius, $`m_p`$ is the mass of the proton) and $`M_0`$ is the mass of the original ejecta. The first term on the left of the equation is the kinetic energy of the ejecta and the second term is the internal energy of the shock. For ultra-relativistic shocks $`\sigma =6/17`$, while $`\sigma =0.73`$ for non-relativistic shocks (Blandford & Mckee 1976) 333Strictly speaking, these values are corresponding to point explosions with shock structure described by the similarity solution. . For simplicity, we approximate this term as $`0.7(\mathrm{\Gamma }^21)m_{\mathrm{sw}}c^2`$ in the following calculation (Note that at observation time the jet has transited to the sub-relativistic motion phase and that a slight difference of $`\sigma `$ may reflect it in the value of $`E_0`$, which is a free parameter here).

The kinematic equation of the approaching jet is

$$\frac{dR}{dt}=\frac{\beta (\mathrm{\Gamma })c}{1\beta (\mathrm{\Gamma })\mathrm{cos}\theta },$$
(2)
where $`v=\beta c`$ is the bulk velocity of the jet with $`\beta (\mathrm{\Gamma })=(1\mathrm{\Gamma }^2)^{1/2}`$, $`t`$ is the observer time and $`\theta `$ is the jet inclination angle to the line of sight.

Given the initial condition at $`t=0`$, which is chosen to be $`R_0=10^7\mathrm{cm}`$ and $`\mathrm{\Gamma }_0=3`$, the two equations can be solved and we can get the relation between the proper motion $`\mu `$ ($`\mu =R\mathrm{sin}\theta /d=R\mathrm{sin}\theta /5.3\mathrm{kpc}`$) and time $`t`$.

We find that the following combination of the parameters fits the observed proper-motion data reasonably well: $`E_0=3.6\times 10^{44}\mathrm{erg}`$, $`n=1.5\times 10^4\mathrm{cm}^3`$, $`\theta _j=1.5^{}`$ and $`\theta =50^{}`$.444$`\theta =50^{}`$ and $`\mathrm{\Gamma }_0=3`$ are consistent with the initial proper motion $`>57\mathrm{mas}\mathrm{day}^1`$ for $`d=5.3\mathrm{kpc}`$.  The model fit is plotted in Figure 1 as the solid line. The late time behavior approaches the well-known Sedov solution $`Rt^{2/5}`$ (see the dashed line in Fig. 1). This is expected since the second term on the left of Eq.(1) becomes dominant at the late time, so $`\beta ^2R^3=\mathrm{constant}`$. In addition, $`\frac{dR}{dt}\beta c`$ when $`\beta \mathrm{cos}\theta 1`$. The earlier phase can be regarded as the transition regime from the mildly relativistic motion to the non-relativistic motion. By comparison, we also plot the $`Rt^{0.5}`$ and $`Rt^{0.3}`$ cases in Figure 2, which show clear deviations from the data.

The inferred value of the ISM density is surprisingly low. Such low densities with $`n10^3\mathrm{cm}^3`$ have been inferred around another two microquasars GRS 1915+105 and GRO J1655-40 by Heinz (2002) from the fact that jets move with constant velocities up to distance $`0.04`$ pc. As argued by Heinz (2002), this implies that either the sources are located in regions occupied by the hot ISM phase or previous activities of the jets have created evacuated bubbles around the sources.

## 4. The radiation model

### 4.1. The forward shock emission

First, we try to fit the X-ray emission of the eastern jet using the forward shock model, the mechanism believed to be responsible for GRB afterglows. In the standard picture of GRBs, an afterglow is generally believed to be produced by the synchrotron radiation or inverse Compton emission of the shock-accelerated electrons in an ultra-relativistic shock wave expanding into the ambient medium. As more and more ambient matter is swept up, the shock gradually decelerates while the emission from such a shock fades down. The microquasar jet is similar to the GRB remnant at the time that its Lorentz factor has decreased to $`\mathrm{\Gamma }13`$, usually months to years after the burst. So, we expect that similar emission processes should also occur in the case of the microquasar decelerating jet.

If the distribution of the shock-accelerated electrons takes a power-law form with the number density given by $`n(\gamma _e)d\gamma _e=K\gamma _e^pd\gamma _e`$ for $`\gamma _m<\gamma _e<\gamma _M`$, the volume emissivity at the frequency $`\nu ^{}`$ in the comoving frame of the shocked gas is (Rybicki & Lightman 1979)

$$j_\nu ^{}=\frac{\sqrt{3}q^3}{2m_ec^2}\left(\frac{4\pi m_ec\nu ^{}}{3q}\right)^{\frac{1p}{2}}B_{}^{\frac{P+1}{2}}KF_1(\nu ^{},\nu _m^{},\nu _M^{}),$$
(3)
where $`q`$ and $`m_e`$ are respectively the charge and mass of the electron, $`B_{}`$ is the strength of the component of magnetic field perpendicular to the electron velocity, $`\nu _m^{}`$ and $`\nu _M^{}`$ are the characteristic frequencies for electrons with $`\gamma _m`$ and $`\gamma _M`$ respectively, and

$$F_1(\nu ^{},\nu _m^{},\nu _M^{})=_{\nu ^{}/\nu _M^{}}^{\nu ^{}/\nu _m^{}}F(x)x^{(p3)/2}𝑑x$$
(4)
with $`F(x)=x_x^+\mathrm{}K_{5/3}(t)𝑑t`$ ($`K_{5/3}(t)`$ is the Bessel function).

The physical quantities in the pre-shock and post-shock ISM are connected by the jump conditions:

$$n^{}=\frac{\widehat{\gamma }\mathrm{\Gamma }+1}{\widehat{\gamma }1}n,e^{}=\frac{\widehat{\gamma }\mathrm{\Gamma }+1}{\widehat{\gamma }1}(\mathrm{\Gamma }1)nm_pc^2,$$
(5)
$$\mathrm{\Gamma }_{\mathrm{sh}}^2=\frac{(\mathrm{\Gamma }+1)[\widehat{\gamma }(\mathrm{\Gamma }1)+1]^2}{\widehat{\gamma }(2\widehat{\gamma })(\mathrm{\Gamma }1)+2},$$
(6)
where $`e^{}`$ and $`n^{}`$ are the energy and the number densities of the shocked gas in its comoving frame and $`\widehat{\gamma }`$ is the adiabatic index, which equals $`4/3`$ for ultra-relativistic shocks and $`5/3`$ for sub-relativistic shocks. A simple interpolation between these two limits $`\widehat{\gamma }=(4\mathrm{\Gamma }+1)/3\mathrm{\Gamma }`$ gives a valid approximation for trans-relativistic shocks (Dai, Huang & Lu 1999).

Assuming that shocked electrons and the magnetic field acquire constant fractions ($`ϵ_e`$ and $`ϵ_B`$) of the total shock energy, we get

$$\gamma _m=ϵ_e\frac{p2}{p1}\frac{m_p}{m_e}(\mathrm{\Gamma }1),B_{}=\sqrt{8\pi ϵ_Be^{}}$$
(7)
and

$$K=(p1)n^{}\gamma _m^{p1}$$
(8)
for $`p>2`$.555 The expressions of $`\gamma _m`$ and $`K`$ are different if $`p<2`$ (see e.g. Dai & Cheng 2001). It is reasonable to believe that $`\nu _M^{}`$ is well above the X-ray band throughout the observations, because the forward shock continuously heats fresh ISM and accelerates electrons. The observer frequency $`\nu `$ relates to the frequency $`\nu ^{}`$ in the comoving frame by $`\nu =D\nu ^{}`$, where $`D=1/\mathrm{\Gamma }(1\beta \mathrm{cos}\theta )`$ is the Doppler factor. The observed flux density at $`\nu `$ and the X-ray flux in the band 0.3-8keV are respectively given by

$$F_\nu =\frac{\theta _j^2}{4}(\frac{R}{d})^2\mathrm{\Delta }RD^3j_\nu ^{};F(0.38\mathrm{k}\mathrm{e}\mathrm{V})=_{\nu _1}^{\nu _2}F_\nu 𝑑\nu ,$$
(9)
where $`\mathrm{\Delta }R`$ is the width of the shock region and is assumed to be $`\mathrm{\Delta }R=R/10`$ in the calculation. The model fit to the light curve of the eastern jet is presented in Figure 3. Clearly, the model light curve decays too slowly to fit the observed data. This is expected from the following analytic study. For $`\nu _m^{}\nu ^{}\nu _M^{}`$, $`F_\nu D^3B_{}^{(p+1)/2}KR^3D^3B_{}^{(p+1)/2}n^{}\gamma _m^{p1}R^3`$. When $`\beta 1`$, $`\beta t^{3/5}`$, $`Rt^{2/5}`$, $`B\beta `$, $`\gamma _m\beta ^2`$ and $`n^{}=4n=\mathrm{constant}`$. So, we get $`F_\nu t^{(15p21)/10}t^{1.2}`$ for $`p=2.2`$. We also consider the case in which the jet spreads laterally with the sound velocity as discussed in GRBs(Rhoads 1999), but find that it makes little difference.

### 4.2. The reverse shock emission

Another probable emitting region for X-rays is the adiabatically expanding ejecta itself. A reverse shock wave that moves back through the original ejecta becomes important when the swept-up matter mass equals $`1/\mathrm{\Gamma }_0`$ of the ejecta mass. The reverse shock accelerates electrons in the ejecta and may amplify the original magnetic field to be close to equipartition. The emission from the non-thermal ($`N(\gamma _e)\gamma _e^p`$) relativistic electrons in the adiabatically expanding ejecta with radius $`R`$ is described by the van der Laan (1966) model, where the flux density is given by $`F_\nu R^{2p}\nu ^{(1p)/2}`$ in the optically-thin regime. So, when $`\beta 1`$, $`F_\nu t^{4p/5}t^{1.76}`$ for $`p=2.2`$. This decay is still too slow to fit the observed light curve of the X-ray jet for which the power-law fit gives a decay index of $`3.7\pm 0.7`$ (Kaaret et al. 2003). However, different from the continuous forward shock, the reverse shock operates only once, so the electrons, including those with the maximum energy in the ejecta, all cool adiabatically. So it is likely that, at some time, the characteristic frequency of the electrons with the maximum energy $`\gamma _Mm_ec^2`$ may fall close to the X-ray band and the X-ray flux would decay quite rapidly since then.

The maximum energy of the power-law distribution electrons just after the reverse shock crosses the ejecta is determined by the shock acceleration process, which is however not well understood. If the electrons in the high energy end of the power law cool faster than the dynamical time scale, the real maximum energy of the electrons in the ejecta is limited by the synchrotron cooling timescale; electrons with energy greater than this energy would have cooled down within the dynamical timescale. So, $`E_M^0=\mathrm{min}(E_{M,acc},E_{M,cool})`$, where $`E_{M,acc}`$ and $`E_{M,cool}`$ denote the maximum energies allowed by the shock acceleration process and the cooling process respectively, the superscript 0 in $`E_M^0`$ denotes the value at $`t_0`$— the time when the reverse shock heats the ejecta. If the synchrotron radiation dominates the cooling of the electrons, $`E_{M,cool}=6\pi m_e^2c^3/(\sigma _TB_{}^2t_0^{})`$, where $`t_0^{}=t_0/D(t_0)`$ is the dynamical time in the comoving frame, $`D(t_0)`$ is the Doppler factor of the X-ray jet at the time $`t_0`$, $`\sigma _T`$ is the Thomson cross section.

The physical quantities in the adiabatically expanding ejecta with radius $`R`$ relate with their initial value at $`t_0`$ by ( van der Laan 1966)

$$\gamma _m=\gamma _m(t_0)\frac{R_0}{R},\gamma _M=\gamma _M(t_0)\frac{R_0}{R};$$
(10)
$$K=K(t_0)\left(\frac{R}{R_0}\right)^{(2+p)},B_{}=B_{}(t_0)\left(\frac{R}{R_0}\right)^2$$
(11)
if the synchrotron emission cooling is negligible, where $`R_0`$ is the radius of the ejecta at time $`t_0`$ and $`\gamma _M(t_0)=E_M^0/m_ec^2`$. These relations are derived from the assumption that total number of the electrons is conserved and that the magnetic field is frozen to the plasma fluid666The same relation for the magnetic field holds if we assume that the energy in the magnetic field constitutes a constant fraction of the internal energy in the ejecta during the whole adiabatically expanding phase. .

From the dynamic model in section 3, we get $`t_0=121`$ days and $`R_0=0.71\times 10^{18}\mathrm{cm}`$ for the eastern X-ray jet from XTE J1550-564. Given the initial condition for $`K(t_0)`$, $`B_{}(t_0)`$, $`\gamma _m(t_0)`$ and $`\gamma _M(t_0)`$, we can obtain the model light curves of the emissions from the expanding ejecta using Eqs.(3), (4), (9) and $`R(t)`$. In the calculation, we have assumed $`\gamma _M(t_0)=E_{M,cool}/m_ec^2=6\pi m_ec/(\sigma _TB_{}(t_0)^2t_0^{})`$. We find that the following combination of the initial values can fit the flux data well (see Figure 4): $`B_{}(t_0)=0.5\mathrm{mG}`$ and $`K(t_0)=0.8\mathrm{cm}^3`$ and $`\gamma _m(t_0)=100`$.777The model result is very insensitive to the value of $`\gamma _m(t_0)`$ For these values, $`\gamma _M(t_0)=3.55\times 10^8`$.

Noting that $`K=ϵ_e(p2)e^{}\gamma _m^{p2}`$ and $`B_{}=\sqrt{8\pi ϵ_Be^{}}`$, the above value for $`K(t_0)`$ and $`B_{}(t_0)`$ imply that the equipartition parameters for energies in electrons and magnetic field in the reverse shock are $`ϵ_e=0.6`$ and $`ϵ_B=0.44\times 10^2`$ respectively, if the internal energy density in the reverse shock is equal to that in the forward shock region at time $`t_0`$. Surprisingly, these equipartition values are close to the typical values inferred for GRB afterglows (Granot et al. 1999, Wijers & Galama 1999; Wang, Dai & Lu 2000). Please also note that here the emitting electrons are in slow cooling regime as $`\nu _m(t_0)10^7\mathrm{Hz}\nu _c`$, where $`\nu _c`$ is the cooling frequency, and that at the X-ray band the synchrotron emission dominates over the inverse Compton emission (Sari & Esin 2001; Panaitescu & Kumar 2000).

We further calculate the model spectrum for the eastern jet on 1 June 2000 using the above parameter values and plot the fit of the observed data in Figure 5. Clearly the model fits quite well the energy spectrum on 1 June 2000, 621 days after the ejection of the eastern jet from XTE J1550-564. At this time, the characteristic frequency of the electrons with $`\gamma _M`$ is just near the X-ray band, so the X-ray spectrum doesn’t become steeper. The model predicts that there should be significant steepness of X-ray spectrum at March 2002, but the relatively few counts do not give a reliable spectral index (Kaaret et al. 2003).

We have shown that the emission from the shocked ejecta provides a viable mechanism for the rapidly decaying X-ray flux from the eastern jet of XTE J1550-564. To guarantee the dominance of the emission from the shocked ejecta over that from the forward shock region during the period of the observations, the forward shock emission should be below the upper limit of the observation made on June 19 2002. Fig.3 tells us that the forward shock emission for $`ϵ_e=0.1`$ and $`ϵ_B=10^4`$ is inferred to be below the limit. Because $`F_\nu \gamma _m^{p1}B^{(p+1)/2}ϵ_e^{p1}ϵ_B^{(p+1)/4}`$, the magnetic field and electron energy fractions in the forward shock region must satisfy (see Fig. 3)

$$\left(\frac{ϵ_e}{0.1}\right)^{p1}\left(\frac{ϵ_B}{10^4}\right)^{(p+1)/4}1.$$
(12)
The high value for $`ϵ_B`$ in the reverse shock region relative to that in the forward shock region can be accounted for if the ejecta has already been magnetized before the further shock compression888The high value of $`ϵ_B`$ inferred for GRB afterglows can, however, be attributed to particular environment (e.g. Königl & Granot 2002) around the burst or a particular magnetic field amplification mechanism (e.g. Thompson & Madau 2000).. This is reasonable as the ejecta from microquasars are suggested to originate from the inner accretion disks (Mirabel & Rodriguez 1999). A pair-rich outflow from the microquasar probably account for the comparatively large value of $`ϵ_e`$ in the reverse shock. Future observations of large-scale jets from microquasars may provide better understanding of the shock physics and physical condition in relativistic jets.

## 5. Conclusions and Discussions

It is believed that relativistic jets exist in accreting systems ranging from Galactic X-ray binaries, gamma-ray bursts and active galactic nuclei (AGNs). Jets in Galactic X-ray binaries such as XTEJ1550-564 evolves much more rapidly than AGN jets and therefore offer a good opportunity to study the dynamical evolution of relativistic jets on time scales inaccessible for AGNs. Although afterglows in GRBs evolves also rapidly, their cosmological distances make the direct measurements of the proper-motion impossible and their dynamics can be studies only indirectly.

The discovery of the extended radio and X-ray emission from the microquasar XTE J550-564 (Corber et al. 2002; Tomsick et al. 2003; Kaaret et al. 2003) represents the first detection of large-scale relativistic jets from a Galactic black hole candidate in both radio and X-rays. These large-scale jets appear to arise from a relatively brief ejection event and, therefore, offer a unique opportunity to study the large-scale evolution of an impulsive jet. We find that the dynamical evolution of the observed eastern jet from the XTE J550-564 is consistent with the well-known Sedov evolutionary phase, during which the energy in the jet is conserved and $`Rt^{2/5}`$. The apparent superluminal motion observed at the very early epoch implies that the initial motion of the jet is at least mildly relativistic. As more and more ISM matter is swept up, the jet decelerates and finally transits to the non-relativistic phase. A trans-relativistic external shock dynamical model is shown to be able to fit the observed proper motion data reasonably well. The inferred ISM density around the jet is $`n1.5\times 10^4\mathrm{cm}^3`$, well below the canonical value. Such low ISM density gains support from the inferred value $`n10^3\mathrm{cm}^3`$ around another two microquasars GRS 1915+105 and GRO J1655-40 by Heinz (2002), from the fact that jets move with constant velocities (i.e. no slowing down) up to distance $`0.04`$ pc. As suggested by Heinz (2002), this implies either that the sources are located in regions occupied by the hot ISM phase or that previous frequent activities of the jets have created evacuated bubbles around the sources.

We first try to fit the X-ray light curve of the eastern jet with the emission from the shocked ISM. However, it is found that this predicts a decay too slow to fit the observations. The model predicts $`F_\nu t^{(15p21)/10}t^{1.2}`$ during the non-relativistic phase, while the power law fit of the X-ray flux data gives $`F_\nu t^{3.7\pm 0.7}`$ (Kaaret et al. 2003). We then turn to another likely emission region—the adiabatically expanding ejecta heated by the reverse shock, quite similar to the mechanism suggested to be responsible for the optical flash and radio flare from GRB990123 (Sari & Piran 1999). Different from the shocked ISM, all electrons in the ejecta cool by adiabatic expansion, so the maximum energy of the electrons in the ejecta decreases as well999 However, for the forward shock, there are always new electrons being shocked. $`\gamma _M`$ is determined by the balance between the acceleration timescale and the cooling timescale: $`\gamma _MB^{1/2}`$ (e.g. Mészáros et al. 1993) . So, for the forward shock, $`\gamma _M`$ even increases with time, since the magnetic field decreases with time. . Once the characteristic synchrotron radiation frequency of these electrons falls close to the X-ray band, the X-ray flux from the ejecta should decay quite rapidly (drops exponentially with time) since then. Using this model, we fitted both the X-ray light curve data and the energy spectrum on 1 June 2000 of the eastern jet successfully (see Figures 4 and 5). One prediction we can make here is a flattening of the light curve at late time, once the forward shock emission, if still above the detection limit, overtakes the reverse shock emission.

The western (receding) jet from XTE J1550-564 was detected in radio and X-rays in 2002 while archival Chandra data on this source from June, August and September 2000 only give upper limits. The non-detection in 2000 is consistent with the deduction that the western source is the receding jet. Unlike the smoothly decaying eastern jet, the western jet brightens at the late time. The brightening might be caused by the inhomogeneities in the ISM (Kaaret et al. 2003) or internal shocks produced by a faster jet overtaking a slower one (Kaiser et al. 2000), and need further careful study. The western jet moved by $`0.52\pm 0.13`$ arc sec between 11 March and 19 June 2002 with a mean apparent speed significantly less than the average apparent speed from 1998 to early 2002. Interestingly, we find that the decay of the X-ray flux of the western jet between March and June 2002 is also consistent with that predicted by the reverse shock emission $`F_\nu t^{4/5p}t^{1.86}`$ for $`p=2.32`$ of the western jet.

Besides the relativistically moving, decelerating jets from XTE J1550-564, large-scale X-ray jets and radio lobes up to $`40\mathrm{arcmin}`$ size have been observed from SS433 (Brinkmann et al. 1996; Dubner et al. 1998). The radiation is suggested to come from the termination shock which results from the interaction of the mass outflow with the nebula W50. Recently, reheating of baryonic material in X-ray jets of SS433 is inferred to take place within $`10^{17}\mathrm{cm}`$ from the core, based on the observed iron emission lines (Migliari, Fender & Méndez 2002). We think that external shock is a possible mechanism for such reheating. As suggested for TeV neutrino emission from microquasar jets (through internal shocks) by Levinson & Waxman (2002), external shocks of microquasar jets may also accelerate the protons in both the shocked ISM and the shocked ejecta. So, they are also potential sources of cosmic-rays (Heinz & Sunyaev 2002), high-energy neutrinos and high-energy gamma-rays.

In summary, we developed a model for the dynamical evolution and radiation of the large-scale X-ray jets from the microquasar XTE J1550-564 analogous to the external shock model for GRB afterglows. In this model, the observed jet emission is due to interaction between the jets and external ISM. Future continuous, follow-up multi-wavelength observations of new ejection events from microqusasrs up to the significant deceleration phase should provide more valuable insights into the nature and physical condition (e.g. shock and particle acceleration physics) of relativistic jets. Owing to the proximity of the Galactic X-ray binaries, further studies on them also offer an exciting way for a better understanding of relativistic jets seen elsewhere in the Universe.

We are grateful to the anonymous referee for his/her valuable suggestions. This work was supported by the National Natural Science Foundation of China under grants 19973003, 19825109 and 10233010, and the National 973 project.