text:
# Extended Air Showers and Muon Interactions.

## I Introduction

When an ultra high energy astroparticle interacts with an atom of the Earth’s atmosphere, it produces a shower of secondary particles that continues interacting and generating more secondary particles. The study of the characteristics of air showers initiated by ultra high energy cosmic rays is of central importance. This is due to the fact that in our days such primary particles cannot be detected directly; instead, they must be studied from different measurements of the air showers they produce.

We have been studying the physics of air showers for several years. We started working on the topic of the electromagnetic processes in air showers analyzing the modifications in the shower development due to the reduction of the electron bremsstrahlung and electron pair production by the Landau-Pomeranchuk-Migdal effect and the dielectric suppression ; and we also studied the influence of the geomagnetic field in an air shower .

The main goal of this work is to analyze other radiative processes that take place during the development of an ultra high energy air shower. We have studied the processes of muon bremsstrahlung, muonic pair production (electron and positron) and muon-nucleus interaction. At high energies these processes become important and dominate the energy losses of the energetic muons that are present in an air showers. The mentioned mechanisms are characterized by small cross sections, hard spectra, large energy fluctuations and generation of electromagnetic sub-showers for the case of muon bremsstrahlung and muonic pair production, and hadronic sub-showers for the case of muon-nucleus interaction. As a consequence, the treatment of such energy losses as uniform and continuous processes is for many purposes inadequate .

We have studied the three mentioned processes concluding that the muon nucleus interaction has less probability than the other ones to produce hard events, and therefore, to generate sub-showers and introduce significant modifications in the air shower development. For this reason, this work will be primarily focused on studying the consequences of the purely electromagnetic processes, namely, muon bremsstrahlung and muonic pair production.

In order to analyze the influence of muon bremsstrahlung and muonic pair production on the air shower observables, we have developed new procedures for these mechanisms and have incorporated them in AIRES (AIRshower Extended Simulations) . Using the data generated with the AIRES code, we have studied the changes introduced by those processes on the different physical quantities.

This work is organized as follows: in section II we briefly review the theory of muon bremsstrahlung, muonic pair production and muon-nucleus interaction. At the end of this section we compare the three effects and analyze under which conditions they can modify the development of air showers. In section III we show the results of our simulations. Finally we present our conclusions and comments in section IV.

## II Theory

### A Theory of Muon Bremsstrahlung

The first approach to the muon bremsstrahlung (MBR) theory was due to Bethe and Heitler . Their results can be reproduced by the standard method of QED  similarly as in the case of electron bremsstrahlung. Bethe and Heitler also considered in their calculation the screening of the atomic electrons.

After this first formulation some corrections were introduced. Kelner, Kokoulin and Petrukhin  also considered the interactions with the atomic electrons. The nuclear form factor was investigated by Christy and Kusaka  for the first time and then by Erlykin . Petrukhin and Shestakov  found that the influence of the nuclear form factor is more important than the predictions of previous papers. These last results have been confirmed by Andreev et al.  who also considered the excitation of the nucleus. In the following paragraphs we give some details of these different approaches for MBR theory.

#### 1 MBR with the effect of the screening by the atomic electrons

It is possible to reproduce the results found by Bethe and Heitler  (in the case of no screening) performing the calculations of the Feynman diagrams (see figure 1) at first order of perturbation theory. For energies that are large compared with the muon mass, the MBR differential cross section integrated over final muon and photon angles takes the form:

$$\sigma (v,E)dv=\alpha \left[\frac{2Zr_0m_e}{m_\mu }\right]^2\left[(2v+v^2)\frac{2}{3}(1v)\right]\frac{dv}{v}$$
(1)
where $`E`$ is the primary energy of the muon, $`k`$ is the photon energy, $`v=k/E`$ is the fraction energy transferred to the photon, $`r_0`$ is the classical electron radius ($`r_0=2.81794092\times 10^{15}`$ m), $`m_e`$ ($`m_\mu `$) is the electron (muon) mass and $`\alpha `$ is the fine structure constant ($`\mathrm{}=c=1`$ throughout this paper).

When the atomic state involved is not changed, the effect of atomic electrons (screening) is taken into account by introducing the elastic atomic form factor in the cross section for bremsstrahlung under the effect of a Coulomb center .

As we have just mentioned before, Bethe and Heitler , took into account in their calculation the influence of the screening. The atomic electrons change the Coulomb potential of the nucleus in the following way :

$$V=\frac{1}{q^2}\left[ZF(q)\right]^2$$
(2)
where $`Z`$ is the charge of the nucleus, $`q`$ is the momentum transferred to the nucleus and $`F`$ is the atomic form factor, that is (in spherical coordinates):

$$F(q)=\rho (r)e^{iqr}𝑑\mathrm{\Omega }$$
(3)
where $`\rho (r)`$ is the density of the atomic electrons at the distance $`r`$ of the nucleus. Bethe and Heitler assume the Fermi distribution for this density that is,

$$\rho (r)=\frac{\rho _0}{1+exp\left[\frac{r\alpha }{\beta }\right]}$$
(4)
where the constants $`\alpha `$ and $`\beta `$ are different for each of the elements. This distribution describes adequately all elements with $`Z10`$. The Fermi radius of the atom is given by:

$$a(Z)=a_0Z^{1/3}$$
(5)
where $`a_0`$ is the Bohr radius.
The screening effect becomes important when $`F`$ in equation (2) is comparable with $`Z`$. This occurs when $`q`$ is of the order (or smaller than) the reciprocal atomic radius, that is:

$$q\frac{Z^{1/3}}{a_0}=\alpha m_eZ^{1/3}$$
(6)
In this case the phase ($`qr`$) in equation (3), and thus $`F(q)`$, are small.

Due to the fact that the differential cross section is proportional to $`1/q^2`$, the largest contribution to the radiation cross section originates from the region where the momentum transferred, $`q`$, is small. Let $`q_{min}`$ be the minimum of $`q`$. Using equation (6), the condition for the screening to be effective reads

$$q_{min}q\alpha m_eZ^{1/3}$$
(7)
The minimum value of $`q`$ occurs when the momentum of the muon is parallel to the emitted photon,

$$q_{min}=\delta =p_1p_2p_k$$
(8)
where $`p_1`$ ($`p_2`$) is the initial (final) momentum of the muon and $`p_k`$ is the momentum of the emitted photon. When the energies considered are larger compared with the muon mass, the last equation reduces to

$$\delta =\frac{m_\mu ^2}{2E}\frac{v}{1v},$$
(9)
and from equations (7), (8) and (9) we can write

$$\frac{E(Ek)}{k}\frac{m_\mu ^2Z^{1/3}}{2m_e\alpha }$$
(10)
It is common to use the ratio between the atomic shell radius (5) and the distance from the nucleus R, as a parameter that gives a quantitative estimation of the importance of the screening effect:

$$\gamma =\frac{a(Z)}{R}.$$
(11)
One can then estimate the distance from the nucleus using the uncertainty relation ($`R=1/\delta `$) and equation (5), so $`\gamma `$ can be written as

$$\gamma =\frac{1}{\alpha m_eZ^{1/3}}\frac{m_\mu ^2}{2E}\frac{v}{1v}$$
(12)
The limits $`\gamma 0`$ and $`\gamma 1`$ correspond respectively to the cases of appreciable and negligible screening effect.

When the influence of the atomic electrons is taken into account, the integration over angles in the differential cross section of muon bremsstrahlung  can only be carried out numerically. Accordingly to the calculation of Bethe and Heitler  the MBR cross section can be written in the following way (for $`Em_\mu `$)

$$\sigma (v,E)dv=\alpha \left[\frac{2Zr_0m_e}{m_\mu }\right]^2\left[(22v+v^2)\varphi _1(\delta )\frac{2}{3}(1v)\varphi _2(\delta )\right]\frac{dv}{v}$$
(13)
$`\varphi _1(\delta )`$ and $`\varphi _2(\delta )`$ are the well known functions displayed in figure 1 of Bethe and Heitler’s paper . Since the relative difference $`|\varphi _1(\delta )\varphi _2(\delta )|/\varphi _1(\delta )`$ remains less than 3% for all $`\delta `$, the approximation $`\varphi _1(\delta )=\varphi _2(\delta )\varphi (\delta )`$ is justified and equation (13) can be put in the following way:

$$\sigma (v,E)dv=\alpha \left[\frac{2Zr_0m_e}{m_\mu }\right]^2\left[\left(\frac{4}{3}\frac{4}{3}v+v^2\right)\varphi (\delta )\right]\frac{dv}{v}$$
(14)
Petruhkin and Shestakov  give an analytical expression for $`\varphi (\delta )`$ for any degree of screening. Their result can be expressed as

$$\varphi (\delta )=\mathrm{ln}\left\{\frac{189m_\mu }{m_e}Z^{\frac{1}{3}}\left[1+\frac{189\sqrt{e}}{m_e}\delta Z^{\frac{1}{3}}\right]^1\right\}$$
(15)
This $`\varphi (\delta )`$ expression differs from the original $`\varphi _1`$ and $`\varphi _2`$ of Bethe and Heitler theory less than 3.3 % for $`\delta =0`$ and about 10 % when $`\delta =2\alpha meZ^{\frac{1}{3}}`$ (this corresponds to $`\gamma =2`$).

#### 2 Correction due to MBR by the atomic electrons.

Another contribution that must be taken into account is the MBR with the atomic electrons . This correction is especially important for light nuclei.

It has been found  that the usual transformation $`Z^2`$ replaced by $`Z(Z+1)`$ in the differential cross section is not accurate enough to take into account adequately the influence of the above mentioned process. Instead, the inelastic atomic form factor must be included  in order to take into account the electron binding within the atom. If $`d\sigma _0(E,k,q)`$ is the differential cross section for MBR by free electrons then the MBR by the atomic electrons is given by :

$$d\sigma (E,k)=Z_0^{q_{max}}F_a^{in}(q)𝑑\sigma (E,k,q)$$
(16)
where $`F_a^{in}(q)`$ is the inelastic atomic form factor.

In the work of Kelner et. al.  the inelastic form factor was calculated according to the Thomas-Fermi model . Those authors have found the following formula that approximates the differential MBR cross section for muon scattering by atomic electrons

$$\sigma (v,E)dv=\alpha \left[\frac{2Zr_0m_e}{m_\mu }\right]^2\left[\left(\frac{4}{3}\frac{4}{3}v+v^2\right)\varphi _e(\delta )\right]\frac{dv}{v}$$
(17)
where

$$\varphi _e(\delta )=\mathrm{ln}\left[\frac{m_\mu /\delta }{m_\mu \delta /m_e^2+\sqrt{e}}\right]\mathrm{ln}\left[1+\frac{m_e}{\delta B^{}Z^{2/3}\sqrt{e}}\right],$$
(18)
and $`e=2.718`$, $`B^{}=1429`$.

#### 3 Correction due to nucleus form factor and nucleus excitation.

The influence due to the nuclear form factor is usually taken into account as a correction to $`\varphi (\delta )`$. Petrukhin and Shestakov  have found (in the case of nucleus with $`Z>10`$) that the modification due to the nuclear size can be accounted for by changing the equation corresponding to point nucleus, that is equation (15), by

$$\varphi (\delta )=\mathrm{ln}\left\{\frac{189m_\mu }{km_e}Z^{\frac{2}{3}}\left[1+\frac{189\sqrt{e}}{m_e}\delta Z^{\frac{1}{3}}\right]^1\right\}$$
(19)
where $`k=3/2`$.

Petrukhin and Shestakov  proved that the inclusion of the elastic nuclear form factor decreases appreciably the differential bremsstrahlung cross section (by approximately by 10-15 % when $`Z>10`$). This result is in agreement with the one of Andreev et. al. .

It is also possible to take into account an additional correction due to the nuclear level excitations, but this contribution amounts to only about 1 % of the elastic one in the case of nuclei with $`Z10`$ .

#### 4 Final result taking into account the different corrections analyzed above.

To define the algorithms to be used in our calculations we have considered the correction due to atomic screening, the MBR by atomic electrons and the nuclear form factor. The final result for the MBR differential cross section can be written as :

$$\sigma (v,E)dv=\alpha \left[\frac{2r_0m_e}{m_\mu }\right]^2Z(Z\varphi +\varphi _e)\left[\frac{4}{3}\frac{4}{3}v+v^2\right]\frac{dv}{v}$$
(20)
where $`\varphi `$ is given by equation (19) and $`\varphi _e`$ is given by equation (18).

The expression for the differential cross section diverges when the photon energy tends to zero (infrared divergence). In order to overcome this mathematical problem, it is usual to put a cut off in the photon energy, $`k_c`$. Therefore, the total cross section for a photon emitted with energy bigger than $`k_c`$ is calculated by

$$\sigma (Z,T,k_c)=_{k_c}^T\frac{d\sigma (Z,T,k)}{dk}𝑑k$$
(21)
Below the cut off $`k_c`$, the mean energy loss by the muon due to bremsstrahlung is given by

$$E_{loss}^{brem}(Z,T,k_c)=_0^{k_c}k\frac{d\sigma (Z,T,k)}{dk}𝑑k$$
(22)
It is worthwhile mentioning that, as in the case of electron bremsstrahlung, the cut off in the photon energy can be naturally introduced if the Landau-Pomeranchuk-Migdal effect and dielectric suppression  are taken into account .

### B Theory of Muonic Pair Production

In the lowest significant order of perturbation theory, the muonic pair production (MPP) is a 4th order process in QED. Two types of diagrams are present, respectively labelled $`\mu `$ and $`e`$ in figure 2. The main contribution to the total cross section and to the energy loss of muons comes from the $`e`$-diagrams. The $`\mu `$-diagrams have to be taken into account if the energy fraction transferred is large .

Similarly as in the case of MBR, there are several corrections that need to be taken into account in the MPP calculations. Such corrections are of the same kind that the ones introduced for MBR. For brevity we are not going to review them in full detail. Instead, we are going to describe the final results which include all the relevant corrections.

Racah  calculated for the first time the MPP cross section in the relativistic region, but without taking into account the atomic and the nuclear form factor. Thereafter, Kelner  included the correction due to the screening of the atomic electrons. The analytical expression for any degree of screening was introduced by Kokoulin and Petrukhin . Those authors also take into account  the correction due to the nuclear form factor. We wish to emphasize that the influence of the nuclear size is more important when the energy transferred to the pair is large . This last case is important for ultra high energy air showers and therefore the nuclear size effect needs to be included in the algorithms used for the simulations.

If the atomic and nuclear form factors are taken into account, the MPP differential cross section can be expressed as :

$$\frac{d^2\sigma }{d\nu d\rho }=\frac{2}{3\pi }Z(Z+1)(\alpha r_0)^2\left(\frac{1\nu }{\nu }\right)\left[\varphi _e+(\frac{m_e}{m_\mu })^2\varphi _\mu \right]$$
(23)
where

$$\nu =\frac{E^++E^{}}{E},$$
(24)
$$\rho =\frac{E^+E^{}}{E^++E^{}},$$
(25)
$`E^\pm `$ is the total energy of the $`e^\pm `$, and

$$\varphi _{e,\mu }=B_{e,\mu }L_{e,\mu },$$
(26)
with

$$B_e=\left[(2+\rho ^2)(1+\beta )+ϵ(3+\rho ^2)\right]\mathrm{ln}\left(1+\frac{1}{ϵ}\right)+\frac{1\rho ^2\beta }{1+ϵ}(3+\rho ^2),$$
(27)
$$B_\mu =\left[(1+\rho ^2)(1+\frac{3\beta }{2})\frac{1}{ϵ}(1+2\beta )(1\rho ^2)\right]\mathrm{ln}(1+ϵ)+\frac{ϵ(1\rho ^2\beta )}{1+ϵ}(1+2\beta )(1\rho ^2),$$
(28)
$$L_e=\mathrm{ln}\left[\frac{CZ^{1/3}\sqrt{(1+ϵ)(1+Y_e)}}{1+\frac{2m_e\sqrt{e}CZ^{1/3}(1+ϵ)(1+Y_e)}{E\nu (1\rho ^2)}}\right]\frac{1}{2}\mathrm{ln}\left[1+\left(\frac{3m_eZ^{1/3}}{2m_\mu }\right)^2(1+ϵ)(1+Y_e)\right],$$
(29)
$$L_\mu =\mathrm{ln}\left[\frac{\frac{2}{3}\frac{m_\mu }{m_e}CZ^{2/3}}{1+\frac{2m_e\sqrt{e}CZ^{1/3}(1+ϵ)(1+Y_\mu )}{E\nu (1\rho ^2)}}\right],$$
(30)
$$Y_e=\frac{5\rho ^2+4\beta (1+\rho ^2)}{2(1+3\beta )\mathrm{ln}(3+1/ϵ)\rho ^22\beta (2\rho ^2)},$$
(31)
$$Y_\mu =\frac{4+\rho ^2+3\beta (1+\rho ^2)}{(1+\rho ^2)(\frac{3}{2}+2\beta )\mathrm{ln}(3+ϵ)+1\frac{3}{2}\rho ^2},$$
(32)
$$\beta =\frac{\nu ^2}{2(1\nu )},$$
(33)
$$ϵ=\left(\frac{m_\mu \nu }{2m_e}\right)^2\frac{1\rho ^2}{1\nu }.$$
(34)
C is equal to 189, and the kinematic ranges of $`\nu `$ and $`\rho `$ are:

$$\frac{4m_e}{E}=\nu _{min}\nu \nu _{max}=10.75\sqrt{e}\frac{m_\mu }{E}Z^{\frac{1}{3}}$$
(35)
$$0=\rho _{min}\rho (\nu )\rho _{max}=\left[1\frac{6m_\mu ^2}{E^2(1\nu )}\right]\sqrt{1\frac{\nu _{min}}{\nu }}$$
(36)
The total cross section for the emission of an $`e^+e^{}`$ pair is:

$$\sigma (Z,T,E_c)=2_{\nu _c}^{\nu _{max}}𝑑\nu _0^{\rho _{max}(\nu )}𝑑\rho \frac{d^2\sigma }{d\nu d\rho }$$
(37)
where $`E_c`$ is the energy cut-off that has to be introduced to overcome the infrared divergence of equation (23). Below the energy $`E_c`$ the process can be treated as a continuous energy loss. The mean value of the energy lost by the incident muon due to $`e^+e^{}`$ pair production with energy below $`E_c`$ is

$$E_{loss}^{pair}(Z,T,E_c)=2E_{\nu _{min}}^{\nu _c}\nu 𝑑\nu _0^{\rho _{max(\nu )}}𝑑\rho \frac{d^2\sigma }{d\nu d\rho }$$
(38)
### C Muon-nucleus interaction

The nuclear interaction of high energy muons is theoretically much less understood than the purely electromagnetic processes studied in sections A and B.

Borog and Petrukhin  calculated a formula for the differential cross section of this process based on Hand’s Formalism  for inelastic muon scattering, and semi-phenomenological inelastic form factor; and includes nuclear shadowing effect. Their final result is given by:

$$d\sigma (E,k)=\mathrm{\Psi }(k)\mathrm{\Phi }(E,v)$$
(39)
where

$$\mathrm{\Psi }(k)=\frac{\alpha }{\pi }A_{\mathrm{eff}}\sigma _{\gamma N}(k)\frac{1}{k},$$
(40)
$$\mathrm{\Phi }(E,v)=v1+\left[1v+\frac{v^2}{2}\left(1+\frac{2\mu ^2}{\mathrm{\Lambda }^2}\right)\right]\mathrm{ln}\left[\frac{\frac{E^2(1v)}{\mu ^2}\left[1+\frac{\mu ^2v^2}{\mathrm{\Lambda }^2(1v)}\right]}{1+\frac{Ev}{\mathrm{\Lambda }}\left[1+\frac{\mathrm{\Lambda }}{2M}+\frac{Ev}{\mathrm{\Lambda }}\right]}\right],$$
(41)
$`k`$ is the energy lost by the muon, $`v=k/E`$, $`m_\mu `$ ($`M`$) is the muon (proton) mass and $`\mathrm{\Lambda }^2=0.4`$ $`\mathrm{GeV}^2`$. The nuclear shadowing effect is taken into account in $`A_{\mathrm{eff}}`$ according to the parameterization of Brodsky

$$A_{\mathrm{eff}}=0.22A+0.87A^{0.89},$$
(42)
where $`A`$ is the atomic mass. The photo nuclear cross section, $`\sigma _{\gamma N}`$, can be approximated by the Caldwell parameterization  on the basis of experimental data on photo-production by real photons:

$$\sigma _{\gamma N}(k)=49.2+11.1\mathrm{ln}(k)+\frac{151.8}{\sqrt{k}}$$
(43)
### D Analysis of the influence of the muonic events for different conditions

In order to estimate the order of magnitude of the different processes in the case of air showers, we have calculated the mean free path (MFP), $`\lambda `$, of the different effects. Each MFP is inversely proportional to the corresponding cross section. In fact $`\lambda =m/\sigma `$ where $`m`$ is the mass of the target nuclei.

Figure 3 shows the MFP’s in $`\mathrm{g}/\mathrm{cm}^2`$ as a function of the kinetic energy of the muon for the cases of MBR, MPP, emission of knock-on (KNO) electrons ($`\delta `$ rays), and muon-nucleus interaction. One can compare the MFP’s of the different muonic events with the depth of the atmosphere ($`1000`$ $`\mathrm{g}/\mathrm{cm}^2`$ for vertical showers and $`9000`$ $`\mathrm{g}/\mathrm{cm}^2`$ for showers with zenith angle $`85^{}`$). The influence of such processes in the development of the shower will be more important for large zenith angles where the total depth of the shower is bigger and the muonic events have more probability to take place. There is another reason to expect that the influence of the effects will be more appreciable for large zenith angles: The muonic component of the showers at ground level becomes very important for zenith angles larger than $`60^{}`$ (see, for example, figure 2 in reference ).

Due to the fact that the MFP’s decrease when the initial energy of the muon is large, it is expected that the influence of the MBR, MPP and muon-nucleus interactions will be more important under these conditions.

Figure 3 also illustrates that the emission of knock-on electrons (KNO) is the most probable process for energies smaller than $`10^6`$ GeV, while for energies larger than this value, the MPP dominates. The MFP of this last process is about three orders of magnitude smaller than the ones of MBR and muon nucleus interaction.

In the case of MBR, the MFP is 1 or 2 orders of magnitude (for very inclined or vertical showers, respectively) greater than the depth of the atmosphere. Therefore, the total probability of the process during the entire shower path remains always small.

The muon-nucleus interaction competes with the MBR, but, as it is shown in figure 4 (where the differential cross section of both processes is plotted versus $`v`$), the muon-nucleus interaction has less probability to generate hard events and then to produce sub-showers. Moreover, the average energy loss of the muon-nucleus interaction, that is, the integral of $`v`$ times the differential cross section (43), is only about 5% of the total energy loss. Due to these facts, the influence of the muon-nucleus interaction in the shower development will be less important than the one of MBR. In consequence, it can safely be assumed that this interaction will not appreciably affect the air shower development and therefore we have not taken it into account our simulations.

With the aim of analyzing the modifications that MBR and MPP may induce in the shower development we have incorporated in AIRES  new procedures for both processes. The corresponding algorithms emulate the formulations described in this section, including all the details that are relevant for the case of air showers. The technique used to implement such algorithms employs in first term very fast approximate calculations of MFP’s, using adequate parameterizations of the corresponding theoretical quantities. Then, the exactness of the procedure is ensured by means of acceptance-rejection tests performed after primary acceptance of the interactions. As a result, the procedures developed for AIRES do not increase significantly the computer time required by the general propagating procedures, while ensuring that the interactions are treated properly.

## III Air Shower Simulations

In order to analyze the influence of MBR and MPP in the development of air showers initiated by ultra high astroparticles we have performed simulations using the AIRES program  with different initial conditions: primary particles (protons, iron nuclei, muons), primary energies (from $`10^{11.5}`$ eV to $`10^{20.5}`$ eV) and zenith angles (from $`0^{}`$ to $`85^{}`$). We compare results of simulations where the MBR and MPP have been taken into account with results obtained in identical conditions but not considering those muonic processes. Notice that both, the emission of knock-on electrons and the muon decay, are always taken into account.

Unless otherwise specified, the geomagnetic field was not taken into account in the simulations, in order to avoid large muon deflections that are present in quasi-horizontal showers.

Hadronic interactions with primary energy greater than 140 GeV were processed using the QGSJET model , while for energies below that threshold, a modified version of the Hillas splitting algorithm  tuned to match QGSJET predictions at 100 GeV, was used.

Due to the fact that the number of particles in an ultra high energy simulation is very large (for example a $`10^{20}`$ eV shower contains about $`10^{11}`$ particles) it is necessary, from the computational point of view, to introduce a sampling technique in order to reduce the number of particles actually simulated. An extension of the so called thinning algorithm, originally introduced by Hillas , is used in AIRES . This technique allows to propagate all particles whose energy is larger than a fixed energy, called thinning energy, $`E_{thin}`$; and only a small representative fraction of the total number of particles is followed below this energy. A statistical weight is assigned to the accepted particles, which is adjusted to ensure that the sampling method is unbiased. The thinning algorithm of AIRES is controlled by two parameters, namely, the thinning energy and the weight limiting factor, $`W_f`$. The quality of the sampling improves when these parameters diminish. The AIRES thinning technique is explained in detail elsewhere . The thinning energy is usually expressed in units of the shower primary energy, and in this case it is named relative thinning.

### A Evolution of single muons

Let us consider first the case of the evolution of a single muon eventually produced during the development of a given shower. We have simulated such muon initiated showers in a representative case: Primary energy $`10^{14}`$ eV and zenith angle $`85^{}`$. One can observe that due to MBR and MPP a muon of such energies may generate secondary showers. This effect is clearly illustrated in figure 5 (a) where the number of electrons and positrons $`(e^+e^{})`$ is plotted versus the slant depth, $`X_s`$. Due to the processes of MBR and MPP the number of electrons and positrons is enlarged with respect to the no MBR-MPP case. When these effects are not taken into account, only KNO and muon decay can affect the propagation of the muon during all its path. Notice that the average number of muons is virtually equal to one during the entire shower development. This can be explained taking into account that: (1) At energies of the order of $`10^{14}`$ eV the muon decay probability is quite small, about 1 % (the mean free path for decay is approximately $`6\times 10^8`$ m, while the length of the atmosphere along a $`85^{}`$ inclined axis is about $`10^6`$ m). (2) The probability of generating additional muons via decay of pions coming from photo-nuclear reactions involving secondary gammas is vanishingly small.

In figure 5 (b) the fraction of energy accumulated by secondary particles relative to the primary energy, that is

$$f_{E_s}=\frac{E_{prim}E_\mu }{E_{prim}}$$
(44)
is plotted versus $`X_s`$. When the MBR and MPP effects are not taken into account the muon almost does not loose energy during all its path ($`0.008`$ GeV/($`\mathrm{g}/\mathrm{cm}^2`$)), while if such effects are considered, the muon energy loss rises up to $`0.3`$ GeV/($`\mathrm{g}/\mathrm{cm}^2`$) at $`10^{14}`$ eV. This is due to the fact that both MBR and MPP have the possibility to produce hard events, responsible for the more significant losses shown in figure 5 (b). On the other hand, when MBR and MPP are disabled, the muon energy loss comes from the emission of KNO electrons (soft and hard) and the muon decay that implies a total loss of less than 0.1 % of the primary energy, even in the case of horizontal showers.

However, the muon decay may affect the first stage of the average shower development and, in fact, it is the responsible of the initial ($`X_s<1500`$ $`\mathrm{g}/\mathrm{cm}^2`$) peak of figure 5 (a) in the MBR-MPP off case. To understand more clearly the origin of this effect, let us consider the probability, $`P_d`$, that the muon of energy $`E`$ decays before undergoing any process of knock-on electron emission (For simplicity we are not taking into account MBR and MPP in this analysis). A straightforward calculation yields,

$$P_d=1_0^{\mathrm{}}\mathrm{exp}\left\{\left[\frac{X}{\lambda _{KNO}}+\frac{\mathrm{}_X}{\mathrm{}_D}\right]\right\}\frac{dX}{\lambda _{KNO}}$$
(45)
where $`\lambda _{KNO}`$ is the knock-on mean free path in $`\mathrm{g}/\mathrm{cm}^2`$, $`\mathrm{}_D`$ is the decay mean free path in meters, $`X`$ is the matter path measured from the location of the particle and along its trajectory, and $`\mathrm{}_X`$ is the metric path corresponding to $`X`$. $`\mathrm{}_X`$ depends also on the location of the muon (represented by its depth $`X_s`$) and the atmospheric model used, and $`\mathrm{}_D`$ depends on the muon energy. Equation (45) can be conveniently evaluated numerically considering a realistic atmospheric model. The results are plotted in figure 6 where $`P_d`$ is represented as a function of $`X_s`$. As expected, $`P_d`$ is always small and diminishes as long as $`X_s`$ grows. At the top of the atmosphere $`P_d6\times 10^4`$. This means that in a batch of, say, $`10^6`$ showers, an average of 600 showers will be initiated by muon decay. Such showers will be characterized by an initial electron (or positron) carrying a significant fraction of the primary energy, and capable of generating a major electromagnetic shower. These electromagnetic showers are responsible for the initial peak that shows up in figure 5 (a). Notice that the maximum number of particles for $`10^{14}`$ eV electromagnetic shower is (roughly) $`10^5`$. When averaging, such showers contribute attenuated by a factor $`P_d`$. This gives $`N_{max}6\times 10^4\times 10^5=60`$, result that is in agreement with the corresponding plot at figure 5 (a).

### B Evolution of air showers

The influence of MBR and MPP in the global observables of air showers has been exhaustively studied using mainly the representative case of a proton primary.

Although the relative frequencies of MBR and MPP in all cases are small compared with other muonic events like KNO, in some conditions the influence of these processes in the development of the shower is not negligible (For example, as we have just mentioned in the last paragraph, the MBR and MPP may generate sub-showers). Figure 7 displays the percentage of muonic events for $`10^{19}`$ eV proton showers with zenith angles of $`45^{}`$ (solid lines) and $`85^{}`$ (dashed lines). The bars correspond, respectively, to KNO, MBR, MPP, and muon decay (MDY). As it is expected, the KNO processes always account for the largest frequency, 96.36 % (96.64 %) in the $`45^{}`$ ($`85^{}`$) case. The other processes are by far less frequent: 2.29 %, 1.32 % and 0.03 % for MDY, MPP and MBR respectively ($`45^{}`$ case). Comparing the percentages of muonic events in the $`85^{}`$ case against the $`45^{}`$ ones, it can be seen that both MPP and MBR rates are slightly increased (about 3 % and 0.05 % respectively), while the MDY relative frequency diminishes.

In figure 8 the energy distributions of the different muonic events are represented. The initial conditions of the shower are the same as in figure 7. In agreement with the MFP’s of figure 3, MBR and MPP occur, in average, at relative large energies. On the other hand, and as expected, MDY takes place at lower energies. Figure 8 shows that the energy spectrum of the muons that undergo the studied events moves slightly towards large energies if the zenith angle of the shower is increased.

We have also studied the frequency of the muonic events as a function of the primary energy. The results are shown in figure 9 where the percentages of muonic events are plotted versus the primary energy. The main characteristic of these plots are the following: For primary energies above $`10^8`$ GeV, all the percentages remain practically invariant. The KNO effect is always the one with maximum relative percentage. The MDY presents a noticeable dependence with the primary energy in the region below $`10^8`$ GeV (from 10 % at $`10^3`$ GeV down to 0.4 % at $`10^8`$ GeV). On the other hand, the MPP grows with the primary energy although the difference between extremes is less significant than in the MDY case (from 0.9 % at $`10^3`$ GeV up to 3.5 % at $`10^8`$ GeV). MBR behaves similarly than MPP, but this fraction is about two orders of magnitude smaller than the MPP one.

The particular behavior of these fractions can be explained considering the characteristics of the energy distribution of the different muonic events, plotted in figure 10 for several primary energies. All the spectra can, in principle, extend up to the primary energy of the shower. When the primary energy is less than $`10^{14}`$ eV, this cutoff is clearly visible in the plots of figure 10. In these cases, the muons generated during the shower have a non negligible decay probability. For primary energies above $`10^{14}`$ eV, the energy distribution of muons broadens, but the spectrum of decaying muons remains bounded in the region $`E10^{12}`$ eV, due to the fact that the decay probabilities become very small for energies above that limit. As a consequence, the total fraction of decaying muons diminishes progressively with the primary energy, as shown in figure 9.

When the primary energy is much larger than $`10^{14}`$ eV, the energy distribution of muons is concentrated in the region of energies lower than that indicative value, and only a small tail extends to higher energies. As a consequence, most significant part of the energy distributions for all the muonic events become almost independent of the primary energy, and so the fractions plotted in figure 9 do not present important changes at the highest primary energies.

The MPP relative fraction depends mainly on the number of high energy muons, which rises significantly with the primary energies for $`E_{prim}10^{15}`$ eV and stabilizes above that energy.

The very small variations in the fractions of figure 9 at the highest energies (increase of MDY and decrease of MPP fractions) can be regarded as secondary effects of the variations of the characteristics of the hadronic processes that take place at the beginning of the shower development. Cross sections and multiplicity of hadronic collisions rising with energy are some of the aspects that need to be taken into account in this sense. A detailed discussion on the characteristics of the inelastic hadronic collisions is beyond the scope of this work; the interested reader can consult reference .

#### 1 Longitudinal development

The following paragraphs contain a description of the modifications induced on different shower observables due to the MBR and MPP effects. We consider first the case of $`3\times 10^{20}`$ eV proton showers inclined $`85^{}`$.

One of the most evident modifications induced by the MBR and MPP effects is the increase of the size of the residual electromagnetic shower produced during the late stages of the shower development (well beyond the shower maximum).

This shows up clearly in figure 11 where the number of gammas (a1) and electrons and positrons (b1) are plotted against $`X_s`$. Notice that, accordingly to our calculations, there are no visible differences in the position of the maximum of the shower ($`X_{max}`$). In order to show the increase of the electromagnetic shower it is convenient to define the relative difference between the cases where the MBR and MPP are or are not taken into account, that is:

$$\mathrm{\Delta }=\frac{\mathrm{N}_{\mathrm{MBR}/\mathrm{MPP}\mathrm{On}}}{\mathrm{N}_{\mathrm{MBR}/\mathrm{MPP}\mathrm{Off}}}1.$$
(46)
$`\mathrm{\Delta }N_\gamma `$ and $`\mathrm{\Delta }N_e`$ have been plotted in figure 11 (a2) and (b2), respectively. For clarity, these plots include only the tail of the showers ($`X_s>2200`$ $`\mathrm{g}/\mathrm{cm}^2`$). It can be noticed that the relative increase of the number of gammas and electrons is about 20 % at the very late stages of the shower development.

Similar plots describing the development of the average energy of gammas and electrons and positrons are displayed in figure 12. The fact that the relative increments of the energies are similar to the corresponding relative increments of particle numbers indicates that the energies of the electromagnetic particles are not substantially modified by the inclusion of the MBR and MPP effects, as expected.

The influence of MBR and MPP is less significant on the muonic component: The number of muons during the development of the shower practically does not change if these effects are considered. The longitudinal development of muon energy appears in figure 13 (a). This plot shows that the energy of the muons diminishes about 3 % at the tail of the shower if the MBR and MPP effect are enabled. It is also observed that the sum of the energies of all muons divided by the average number of muons,

$$\xi _\mu =\frac{E_\mu }{N_\mu }$$
(47)
is not significantly modified when the effects are considered for the primary energy mentioned above (see figure 13 (b)).

We have also investigated whether or not the MBR and MPP generates modifications in the shower front arrival time profile for different particles of the shower (muons, electrons and gammas). We have not found any significant alteration when comparing the cases when the effects of MBR and MPP are or not taken into account.

We have studied the modifications introduced by the MBR and MPP effects for different primary energies. The influence of these processes in the electromagnetic component of the showers with smaller primary energy is similar to the case of $`3\times 10^{20}`$ eV, described above.

In figure 14, $`\xi _\mu `$ is plotted versus $`s=X_s/X_{max}`$ for different primary energies. All the curves show a similar behavior: (i) In the region $`0<s<1`$ $`\xi _\mu `$ decreases with $`s`$, direct consequence of the multiplicative processes that take place in this phase and increase the number of shower secondaries, thus reducing the average energy per secondary. (ii) For $`s>1`$, low energy muons decay progressively, and therefore the mean muon energy is shifted as long as $`s`$ grows. When comparing the curves corresponding to different primary energies, it is possible to see that $`\xi _\mu `$ increases monotonically as long as the primary energy decreases from $`E_{prim}=3\times 10^{20}`$ eV (a) to $`E_{prim}10^{14}`$ eV (c); and decreases when the energy continues decreasing below $`10^{14}`$ eV (curves (c) and (d) for $`3\times 10^{11}`$ eV; no intermediate cases were plotted for simplicity). These behavior is consistent with the characteristics of the energy distributions of figure 10, already explained at the beginning of section III B: The low energy range (curves (c) and (d)) are characterized by muon spectra bounded by the primary energy and thus significantly changing when it varies, and with a mean value increasing with the primary energy. On the other hand, the high energy range (curves (a) to (c)) is characterized by muon spectra weakly correlated with the primary energy, and enhanced fraction of low energy muons at the highest primary energies.

When comparing the curves corresponding to the cases MBR-MPP disabled (dashed lines) and enabled (solid lines), it is possible to notice that the differences between pairs of curves is always small, with a maximum of 5 % for curve (c) at $`s=10`$. In general $`\xi _\mu `$ decreases when MPP and MBR are switched on. However, a critical combination of event probabilities (see figure 9) determines that $`\xi _\mu `$ remains unchanged or is slightly increased for primary energies around $`10^{15}10^{17}`$ eV.

From figure 14 one can see that in the region around the shower maximum, that is, where $`s`$ ranges approximately between 0 and 2, $`\xi _\mu `$ practically does not present significant changes when comparing the cases where the MBR and MPP are enabled or disabled. On the other hand, progressively significant differences can appear for $`s`$ larger than 2. $`s`$ can be regarded as a measure of the stage of the shower development, ranging from 0 at the top of the atmosphere to a final value $`s_g`$ at ground which depends on the zenith angle. $`s_g`$ is a measure of the quantity of matter the shower has to pass through from its beginning until reaching the ground level. From the plot in figure 14, it can be inferred that when $`s_g`$ is less than 2, that is, for zenith angles less than $`45^{}`$ in the conditions of our simulations, there will be no noticeable modifications on the shower development due to MBR and MPP.

We have also analyzed the influence of MBR and MPP in the showers initiated by different primary particles like, for example, iron nuclei. We have observed that the modifications that MBR and MPP introduce in the late stages of the shower development have approximately the same characteristics of the ones introduced for a proton shower for the case of the electromagnetic component of the mentioned shower. For the case of muons observables the differences are less significant. For example, $`\mathrm{\Delta }_{E_\mu }`$ reaches a maximum of 1.5 % in the late stages of the shower development for the case of iron shower of $`3\times 10^{20}`$ eV and $`85^{}`$ of zenith angle while for proton showers, in the same initial conditions, $`\mathrm{\Delta }_{E_\mu }`$ is 3 %.

#### 2 Lateral distributions

We have also studied the distributions of the particles at ground level (lateral distributions). In the case of very inclined showers, the intersection with the ground plane occurs well beyond the shower maximum, and lateral distributions are somewhat different with respect to the “typical” distributions corresponding to showers with small zenith angles. Among other differences, we can mention: (i) Substantially smaller number of particles. (ii) The densities of the electromagnetic particles are of the same order of magnitude than the density of muons (In the case of quasi vertical showers, the muons account for only about 1 % of the ground particles).

In figure 15 the densities of $`\gamma `$, $`e^+e^{}`$, and $`\mu ^+\mu ^{}`$, are plotted versus the distance to the shower axis, for the case of $`10^{19}`$ eV proton showers inclined $`70^{}`$. The ground level is located at 875 $`\mathrm{g}/\mathrm{cm}^2`$ . The analysis of the data shows that the number of $`\gamma `$ and $`e^+e^{}`$ is slightly modified $``$when MBR and MPP are taken into account$``$ near the shower axis, while the lateral distribution of muons remains virtually unaltered when the MBR and MPP interactions are enabled. It is worthwhile to mention that in this case the geomagnetic field is taken into account in order to simulate a real situation. In fact, we have chosen the conditions corresponding to the site of El Nihuil (Mendoza, Argentina) with the aim of studying the characteristics of showers to be measured by the future Auger Observatory  that is currently being constructed at that site.

The measurement of the lateral distributions of particles at ground generated by ultra high energy cosmic rays is usually performed by means of water C̆erenkov detectors. Such devices are sensible to both electromagnetic particles and muons, and the signal they produce is the sum of both components. The signal produced by any particle hitting a water C̆erenkov detector must be estimated by a specific Monte Carlo simulation which takes into account the characteristics of the detectors. The detailed simulation of water C̆erenkov detectors is beyond the scope of our work; we have instead evaluated estimations of these signals using a direct conversion procedure that retrieves average signals . Such averages were evaluated on the basis of simulations performed with the AGASIM program .

We have plotted in figure 16 the ratio between the electromagnetic and the muonic component of the detector signal, as a function of the distance to the shower axis, for the cases when MBR and MPP are or not taken into account in the simulation of the showers.

The increase of the size of the residual electromagnetic shower that takes place when MBR and MPP are enabled, produces a larger signal close to the shower axis, as it clearly shows up from the plots of figure 16. For distances to the axis less than 30 m the electromagnetic to muon signal ratio increases slightly when the MBR and MPP are switched on. On the other hand, this increment is smaller for larger distances, and is virtually negligible beyond 200 m from the shower axis.

It is worthwhile to remark that the fact that the main modifications in the electromagnetic to muon signal ratios are concentrated in a narrow zone around the shower axis makes it unlikely that the incorporation of MBR and MPP in the air shower simulation engine will significantly affect the results obtained in references  where the measurements of inclined showers performed at the Haverah Park experiment  are analyzed with the help of air shower simulations using AIRES without taking into account MBR and MPP.

Notice also that the data plotted in figure 16 corresponding to the case when the MBR and MPP are switched off can be compared with the corresponding data presented in reference . It is easy to see that there is a good qualitative agreement between the two sets of data and that the small differences between the two works most probably come from differences in the ground altitude and/or specific parameters used to calculate the detector responses.

## IV Conclusions

We have studied the influence of the MBR and MPP in the development of air showers initiated by ultra high energy astroparticles. We have incorporated in the AIRES air shower simulation system the corresponding procedures to emulate these effects and have then performed simulations in a number of different initial conditions.

The analysis of the evolution of a single muon indicates that such particle can eventually generate secondary showers after undergoing hard MBR and MPP processes. This indicates clearly that these interactions cannot be simulated accurately as continuous energy loss processes.

For $`3\times 10^{20}`$ eV proton and iron primaries the main modifications introduced by MBR and MPP affect the electromagnetic component of the showers. The number and energy of gammas and electrons increase significantly in the late stages of the shower development (well beyond the shower maximum), but the mentioned effects do not generate visible changes in the position of $`X_{max}`$.

The changes generated by MBR and MPP for muon observables are less significant: The number of muons practically does not change and their energies diminish about 3 % (1.5 %) for the case of proton (iron) showers at the tail of the shower.

The shower front arrival time profile does not present modifications due to the MBR and MPP processes.

For primary energies below $`3\times 10^{20}`$ eV the modifications in the electromagnetic shower induced by MBR and MPP are qualitatively similar to the ones described in the preceding paragraphs. In the case of muon observables like $`\xi _\mu `$ we have found, in the entire range of primary energies that was considered, small variations due to MBR and MPP. Such small modifications correspond, in general, to decrements in the average muon energies when MBR and MPP are switched on. However, critical combinations of event probabilities determine that $`\xi _\mu `$ can remain unchanged or be slightly increased for primary energies around $`10^{15}10^{17}`$ eV.

The fact that the alterations in the electromagnetic showers are only significant in the late stages of the development of the showers, i.e., $`X_s>X_{max}`$, implies that in normal conditions there will be no visible changes in the electromagnetic shower at ground level, for showers with zenith angle less than $`45^{}`$.

For showers with zenith angles larger than $`60^{}`$, the MBR and MPP processes are responsible for an increment of the density of electromagnetic particles at ground, which is most important in a narrow region around the shower axis. In connection with this result, we have also found that the signal produced by C̆erenkov detectors will be also larger near the shower axis if the mentioned effects are taken into account.

## V Acknowledgments

This work was partially supported by Consejo Nacional de Investigaciones Científicas y Técnicas, Agencia Nacional de Programación Científica of Argentina, FOMEC program, and Fundación Antorchas. We wish to thank Professors H. Fanchiotti and C. A. García Canal (University of La Plata) for enlightening discussions. We are also indebted to Dr. C. Hojvat (Fermilab, USA) for helping us to access many of the references cited in this work.


format: text

text:
# The Physics of Blazar Optical Emission Regions II: Magnetic Field Orientation, Viewing Angle and Beaming

## 1. Introduction

Jet are very important for AGN unification. All radio loud objects have powerful jets. In the orientation Unification Scheme, blazars, with their extreme flux and polarization variability, have the jets beamed towards us. However viewing angle unification is not the complete story. Do all the jets have similar physical properties on average, so that orientation is able to explain most of the variations of observed properties among blazars?

Optical polarization of blazars provides magnetic field information very close to the jet formation region and is therefore a very important tool for probing jet physics. We have investigated a sample of 51 blazars with many optical polarization measurements made over the past 20-30 years. In paper I, we have investigated the relationship between the optical polarization position angle and the VLBI radio structure. The results are consistent with a shock model. In this paper, we focus on the variability of optical polarization position angles. Assuming a simple shock model, we link the angle variability with jet orientation and speed to explore the jet physical parameters among blazars.

## 2. Data and Analysis

### 2.1. Optical Polarization

Since a blazar’s optical polarization swings over a wide range of angles, there exists a 180 ambiguity. Linear statistics based on normal distributions is not strictly correct. Therefore we use circular statistics – a technique little used in astronomy but commonly used for example, to analyze directions of returning homing pigeons, or wind directions. Thus, from the polarization position angles $`\theta _i`$, we calculate the circular standard deviation $`c(\theta )`$ about the mean angle, to represent their scatter (Fisher 1993). Polarization data are axial, i.e. only cover the range $`0^{}`$ to $`180^{}`$. However, the $`2\theta _i`$ are truly circular data. Therefore, we first calculate $`c(2\theta )`$, then divide by 2 to get $`c(\theta )`$:

$$c(\theta )=\left[\mathrm{ln}\left(\frac{\sqrt{(_{i=1}^n\mathrm{cos}2\theta _i)^2+(_{i=1}^n\mathrm{sin}2\theta _i)^2}}{n}\right)\right]^{\frac{1}{2}}$$

where $`n`$ is the total number of measurements. If the sample has a large intrinsic scatter, the determination of the real scatter depends on the accurate shape of the wrapped distribution between $`\pi /2`$ and $`\pi /2`$. If the sampling of data is not enough to reveal the detailed shape of the distribution, the $`c(\theta )`$ calculation algorithm gives a maximum $`c(\theta )66^{}`$.

### 2.2. Radio Core Dominance

In the unification scheme of powerful radio sources, the core-dominated sources are the Doppler-boosted jet, viewed end-on. At larger angles the core becomes much weaker and collimated jets link the well resolved double lobes. Doppler boosting in core-dominant sources is supported by their very high brightness temperatures and by superluminal motions. Statistical analyses of radio source samples support the unified scheme (M. Chiaberge, this conference). The ratio of core-to-lobe luminosity, or the core dominance, R, is therefore a good indicator of viewing angle. However, in blazar samples, viewing angle is not the only parameter affecting Doppler boosting. From work of other authors (Orr and Browne 1982), we know that R depends on three physical parameters and can be expressed as

$$\mathrm{log}\text{R}=\mathrm{log}(\text{R}_T\times (1\beta \times \mathrm{cos}(\varphi ))^3)$$

where $`\varphi `$ is the viewing angle of the jet, $`\beta `$ is the bulk speed of the emitting electrons in the jet and $`R_T`$ is the tangential value of R measured for an average edge-on jet in radio galaxies. We used $`R_T=0.0024`$ at 5GHz rest frequency (e.g. Figure 1 of Hoekstra, Barthel & Hes 1997). While R is of course sensitive to $`\varphi `$, for smaller $`\varphi `$ (e.g. $`<20^{}`$), it is increasingly sensitive to $`\beta `$.

For each object we calculated R at 5 GHz rest frequency, using core and lobe flux densities from literature. In Section 4, we will calculate log R from models using the formula above and compare the models with observed values.

## 3. Correlation

Figure 1 (a) shows that there is a correlation between $`c(\theta )`$ and log R. The correlation is significant at the 99.9% level. However, we note that the correlation is not linear. For the log R$`<`$2 region, which includes most blazars, $`c(\theta )`$ increases with log R, but for log R$`>`$2, there is a lot of scatter. No object has $`c(\theta )`$ larger than 66, as expected from the algorithms for large scatter and small sample size.

## 4. Model

The model we use makes the following assumptions:

1. Optical emission comes from a transverse shock propagating along the jet.
2. The shock compresses the magnetic field so much that the magnetic component along the jet direction is negligible.
3. In the plane of the shock, there is a dominant magnetic field direction that changes randomly with time. This magnetic field direction corresponds to the emission region whose radiation happens to be boosted towards us. Different ejection directions observed for emerging blobs on VLBI scales and complex VLBI structure suggests that these dominant regions change rapidly.
4. The shock is optically thin.

When the shock is viewed face on, the magnetic field direction projected on the sky will vary randomly with time. As the viewing angle is increased, the projected magnetic field direction will lie increasingly perpendicular to the projected jet direction.

A simplification is afforded by the Lorentz invariance of the ratios of Stokes parameters, so our calculations of projected magnetic field direction in the co-moving frame (viewing angle $`\varphi ^{}`$) can be referred to the observer’s frame using the relation (Bjornsson, C.-I. 1982):

$$\mathrm{sin}(\varphi ^{})=\frac{\mathrm{sin}(\varphi )}{1\beta \mathrm{cos}(\varphi )}\sqrt{1\beta ^2}$$

For a given set of $`\beta `$ and $`\varphi `$, we generate 100 randomly oriented magnetic field directions to simulate the random changes over time inside the shock plane and project them onto the sky in the co-moving frame. The circular standard deviation of those angles can then be compared with real observations, using the relation between R and $`\varphi `$ given earlier. $`\text{R}_T`$ is fixed by observations, so the only free parameter in this model is $`\beta `$.

The results are shown in Figure 1 (b). Each curve represents a model with fixed $`\beta `$, with $`\varphi `$ decreasing from bottom left to top right. The values of $`\varphi `$ are marked along the curve. From left to right, $`\beta =0.93,0.95,0.97,0.99`$. The upper limit on $`c(\theta )`$ is the result of the computational limit. The lower limits are caused by the fact that we only have 100 simulated vectors to do statistics and any one vector that happens to be parallel to the jet direction contributes a lot to the scatter. In principle, if we had an increasing number of vectors, the lower limit for $`c(\theta )`$ would approach 0. On our model, the $`c(\theta )`$ vs. log R correlation arises because $`\beta 0.95`$ for most blazars, and so the viewing angle dependence dominates for log R$`<2`$. For log R$`>2`$, the scatter is the result of a tail to higher jet speeds.

Our simple model with only one free parameter, $`\beta `$, accounts for the correlation we find between the scatter in optical polarization angles, c($`\theta `$), and the core-dominance, log R. Figure 1 shows that the average jet speed differs significantly amongst the powerful blazars, with a range of values comparable with those derived from superluminal motion and from the statistics of core-dominance in radio-source surveys (Urry & Padovani 1995). In Unified Schemes it is important to take into account that log R is not simply a measure of orientation. Jet speed introduces considerable scatter, especially at log R $`>`$ 2.

## References

Bjornsson, C.-I. 1982, ApJ, 260, 855

Fisher, N I. 1993, Statistical Analysis of Circular Data, Cambridge University Press.

Hoekstra, H., Barthel, P. D. & Hes, R. 1997, A&A, 319, 757

Orr, M. J. L. & Browne, I. W. A. 1982, MNRAS, 200, 1067

Urry, C. M. & Padovani, P. 1995, PASP, 107, 803


text:
# DISK-JET CONNECTION IN THE MICROQUASAR GRS 1915+105 AND IR AND RADIO EMISSION

## 1. INTRODUCTION

Two Galactic X-ray transient sources GRS 1915+105 and GRO J1655-40 are known to produce relativistic radio jets (Mirabel & Rodriguez (1994); Tingay et al. (1995)). The combination of relativistic jets and a suspected central black hole has earned these two objects the name “microquasars” as they seem to be stellar mass analogs of the massive black hole systems in quasars and other active galactic nuclei (AGNs) (Morgan et al. (1997); Belloni et al. 1997b ; Orosz & Bailyn (1997)). Since microquasars are much smaller, closer and show faster variability than the extragalactic systems, they are potential “laboratories” for the study of black hole accretion/relativistic jet systems. GRS 1915+105 has shown spectacular X-ray variability since its discovery in 1992 (Castro-Tirado et al. (1994)). Recently, Belloni et al. (2000) have classified a large sample of RXTE/PCA observations in 12 separate classes on the basis of their light curves and the color-color diagrams. Out of these, the $`\beta `$ class is described as the most complex, has all the three basic states of the source (Belloni et al. (2000)), and is always accompanied by IR/radio flares. The presence of a “spike” in the X-ray light curve, which separates the dips with hard and soft spectra, clearly distinguishes this class from others. The “spike” coincides with the beginning of IR flares seen during simultaneous X-ray and IR observations, suggesting its role in initiating IR flares (Eikenberry et al. (1998)). These flares are termed as “baby jets” from energy consideration. Another simultaneous observations of GRS 1915+105 in the X-ray, IR, and radio wavelengths confirm that the IR and radio flares are associated with the X-ray dips (Mirabel et al. (1998)).

At present, the disk-jet connection is indirect at best for relativistic radio jets. Harmon et al. (1997) have shown a long-term correlation between hard X-ray flux and jet activity in GRS 1915+105. Fender et al. (1999) have observed four relativistic radio ejections in 1997 October/November and have presented high quality MERLIN radio images of ejecta in 400$``$5000 AU scales which are consistent with ballistic motion. Recently, Eikenberry et al. (2000) have reported detection of faint IR flares and have classified the ejection events associated with IR and radio flares into three classes: (A) the relativistic events producing bright superluminal radio jets of $``$ 1 Jy with decay times of several days (Mirabel & Rodriguez (1994) ; Fender et al. (1999)), (B) the “baby jets” associated with $``$ 100$``$200 mJy IR and radio flares with decay times of several minutes (Eikenberry et al. (1998); Mirabel et al. (1998)), and (C) faint IR flares with peak amplitude of $``$ 0.5 mJy and duration of 8-10 minutes (Eikenberry et al. (2000)). In this paper we present the results of our detailed analysis of RXTE/PCA observations spread over a year of $`\beta `$ and $`\lambda `$ classes which actually belong to the same class as shown here. We investigate the role of the “spike” in initiating the “baby jets” and study the accretion disk corresponding to IR/radio flares of varying peak amplitudes (Table 1). We present evidence of a direct disk $``$ jet connection for radio ejections of class A. We discuss the dissipation of accretion energy in the inner part of the disk in terms of a comptonising cloud and the escaping mass which produces synchrotron radiation.

## 2. OBSERVATIONS AND ANALYSIS

The observations discussed here are selected from the publicly available RXTE/PCA data for the X-ray transient source GRS 1915+105 (Jahoda et al. (1996)). In Table 1 we list details of these observations. The source was in a high/flaring state during these observations. A portion of 2$``$13 keV light curves for different days added for all available PCA units are shown in Figure 1. The “spike” is present in all the observations shown in Figure 1 (all of them belong to class $`\beta `$) and IR/Radio flares were observed during these observations. The peak amplitude of IR/radio flares are given in Table 1. The IR flares with a spacing of $``$30 minutes were observed during simultaneous X-ray/IR observation on 1997 August 14 (Eikenberry et al. (1998)). During another simultaneous observation in the X-ray, IR, and radio on 1997 September 9, similar flares were seen in both the IR and radio bands (Mirabel et al. (1998)). Fender et al. (1999) have observed four relativistic radio jets during 1997 October/November and have seen radio flares on 1997 October 30-31 in 15 GHz band (see their figure 7). The other observation of 1996 October 7 belongs to class $`\lambda `$ (shown in Figure 2) which is a “radio-quiet” state. Muno et al. (1999) have defined the “radio-quiet” state when the radio flux in the 8.3 GHz band (or the 15 GHz band) is less than 15 mJy. The general approach of our analysis is described elsewhere (Yadav & Rao (2000)). The two X-ray colors HR1 and HR2 are defined as the ratio of the flux in the 5$``$13 keV band to flux in the 2$``$5 keV band and the ratio of the flux in the 13$``$60 keV band to flux in the 2$``$13 keV band respectively. The sub-second rms variability is calculated from 0.1 s data.

To study the spectral properties during these observations, we have used simultaneous fits to the PCA data in the energy range of 5$``$60 keV and the HEXTE data in the energy range of 15$``$100 keV. We used “Standard 2” mode PCA data (16 s time resolution) with 128 channel spectra and 1% systematic error added. The 64 channel archive mode data of HEXTE (cluster 0) was rebinned to improve the statistics of the counts. We have used a sum of a multicolor disk black body (DISKBB) and a power law to approximate the energy spectrum (Yadav (2000)).

## 3. RESULTS AND INTERPRETATION

The light curves seen on 1997 October 31 (class $`\beta `$) and on 1996 October 7 (class $`\lambda `$) are shown in the top panels of Figure 2. The HR2 color is shown in the middle panels and sub-second mean variability is shown in the bottom panels of Figure 2. In the top panels, durations of the low hard state (quiescent phase), the high soft state (burst phase) and the soft dips (low soft state) are marked by ‘C’, ‘B’ and ‘A’, respectively. Although Belloni et al. (2000) have found short soft dips of state ‘A’ in the class $`\lambda `$ (in the beginning of burst phase and during the fast fluctuations), we have marked it whole as state ‘B’. This is done to emphasize the long dips of state ‘A’ found in class $`\beta `$. The HR2 is high (0.09$``$0.12) during the quiescent phase (state C) whereas it is low (0.03$``$0.06) during the burst phase (state B) and the soft dips (state A). The similarity of the two plots in the middle panels suggests that the “spike” in the light curves of class $`\beta `$ marks the beginning of the burst phase. This is shown by a vertical arrow in the top left panel. The mean variability plotted in bottom panels is consistent with this suggestion. The mean variability is 2$``$5 % during the burst phase, 10$``$15 % during the quiescent phase and 5$``$9 % during the soft dips. In the region of fast fluctuations during phase B of both the classes, the mean variability is mostly in the range of 5$`9`$% indicating the presence of short soft dips during these fast oscillations. These results suggest that classes $`\lambda `$ and $`\beta `$ of Belloni et al. (2000) are very similar except for (1) the presence of long soft dip in the beginning of the burst phase in case of class $`\beta `$ while it is short and shallow in case of class $`\lambda `$ (see figure 4 of Belloni et al. (2000)), and (2) the short but frequent soft dips during the fast oscillations are relatively deeper and for longer times in case of class $`\beta `$ than that in class $`\lambda `$. It may be noted here that the HR2 distinguishes the quiescent phase from other two phases (the burst phase and the soft dips) while the mean variability highlights the difference between the burst phase and the soft dips.

The class $`\beta `$ has been studied simultaneously in IR and radio bands. Eikenberry et al. (1998) have compared X-ray and IR profiles during 1997 August 14 observations and have found that the spikes in the X-ray band coincide with the beginning of IR flares (see their figure 3). The IR intensity reaches its peak flux shortly after the X-ray peak. The IR flares have decaying phases very similar in their smoothness and time scale to the rising phases. On the other hand, the X-ray flux begins the fast oscillations after reaching their peak level as mentioned earlier. An excess in IR is seen during the fast X-ray oscillations. This excess in IR is explained in terms of many faint IR flares superposed on one another (Eikenberry et al. (2000)). If we assume that each X-ray oscillation has an associated faint IR flare, an IR excess of 1.3 mJy is estimated, which is consistent with the observed excess of $``$ 1.0 mJy for 1997 August 14 observations. During another simultaneous observation of $`\beta `$ class in IR and radio on 1997 September 9, Mirabel et al. (1998) confirmed the role of the spike in initiating IR flares and have found that the radio flares follow the IR flares with a time delay consistent with broad band synchrotron emission. The profiles of IR and radio flares are quite similar. In contrast, the $`\lambda `$ class belongs to radio-quiet state. Pooley & Fender (1997) have observed GRS 1915+105 in radio during 1996 October but no radio flares were reported during observations of $`\lambda `$ class X-ray activity on 1996 October 7.

To investigate further the similarities and differences between $`\beta `$ and $`\lambda `$ classes, we have studied the mean X-ray color HR1, which is a measure of soft photons and hence the disk temperature. On the other hand, HR2 is a measure of hard photons and a indicator of appearance/disappearance of the advective/halo disk (Yadav et al. (1999)). The mean HR1 during the quiescent phase is plotted as a function of mean quiescent flux in Figure 3 for all the observations in Table 1. The mean HR1 increases with the quiescent flux and follows a single path. The HR2 during the quiescent phase also shows similar results except it decreases with the quiescent flux (in the inset of Figure 3). These results reinforce the above suggestion that classes $`\beta `$ and $`\lambda `$ are very similar and actually belong to a single class.

We have studied the spectral properties of both $`\beta `$ and $`\lambda `$ classes using simultaneous fits to PCA and HEXTE data in the energy range 5$``$100 keV as described in the previous section. The segments of the quiescent phase, the burst phase, and the soft dips are analysed separately. We have selected $`\beta `$ class data of 1997 October 30 (Figure 1c) for spectral analysis as the first relativistic jet observed by Fender et al. (1999) was close to this X-ray observation, and also the $`\lambda `$ class data of 1996 October 7 (Figure 2 right panels). The characteristic radius of the inner disk, Rin is derived from the best-fit spectral parameters and the known distance and inclination (the latter is taken as the same as the inclination of the radio jets to the line of sight). From Rin and Tin, we infer the mass accretion rate using $`\dot{m}=8\pi R_{in}^3\sigma T_{in}^4/3GM`$ where $`\sigma `$ is Boltzmann constant, G is the gravitational constant and M is the mass of the black hole. The value of Rin, which is smaller during the burst phase, cannot be less than the innermost stable orbit around a black hole. For a Schwarzschild black hole ($`R_{in}6GM/c^2`$), an Rin of $``$ 20 km during the burst phase puts M $``$ 2.4 M (Belloni et al. 1997a ). The best-fit parameters are given in Table 2. Note that the derived values for the inner disk radius can be an underestimate due to scattering effects and the approximations made in fitting the Comptonised part of the spectrum as a power-law (see Shrader & Titarchuk 1998). We use these numbers here only for a qualitative description.

The results in Table 2 are consistent with the general picture: the spectrum is hard during the quiescent phase ($`\mathrm{\Gamma }`$ 2.5) and it is soft during the burst phase ($`\mathrm{\Gamma }`$ 3.6). The soft X-ray luminosity (5$``$15 keV) increases over a factor of ten during the burst phase while the hard X-ray luminosity (20$``$100 keV) reduced to half than that during the quiescent phase for 1996 October 7 observations. The total X-ray luminosity (5$``$100 keV) drops by a factor of $``$ 5.5 during the quiescent phase than that in the burst phase. Yadav et al. (1999) have suggested that the change of states during the flaring of the source (like during class $`\lambda `$) is due to the appearance/disappearance of an advective/halo disk. If we assume the above scenario, most of the energy ($``$ 80%) disappears during the quiescent phase through advection into black hole (Narayan et al. (1997)). During the quiescent phase, the total X-ray flux in the 5$``$100 keV band is increased from a value of 7.1 $`\times `$ 10-9 erg s-1 cm-2 on 1996 October 7 to a value of 16.4 $`\times `$ 10-9 erg s-1 cm-2 on 1997 October 30 while the B/C ratio for 5$``$15 keV X-ray is reduced from 10.16 to a value of 3.01 (see Table 2). The simplest explanation of these results may be that the disk disrupts due to mass ejection on 1997 October 30 before X-ray luminosity could reach its peak value during the burst phase. During the burst phase the 5$``$100 keV X-ray flux is $``$ 3.4 $`\times `$ 10-8 erg s-1 cm-2 on 1997 October 30 whereas it is $``$ 3.9 $`\times `$ 10-8 erg s-1 cm-2 on 1996 October 7 when short, soft dips were observed by Belloni et al. (2000) as discussed earlier which suggest the beginning of mass ejection (radio flux is still low, see Table 1). These results indicate a luminosity threshold to start the disk evacuation which produces the long, soft dip after a spike in the X-ray light curve.

The B/A luminosity ratio is $``$ 2 in the whole range of 5$``$100 keV on 1997 October 30 (see Table 2) while $`\mathrm{\Gamma }`$ increases to 3.77 during the soft dip from a value of 3.57 during the burst phase suggesting a portion of the inner disk and the halo are blown away/disappears during the soft dips. Feroci et al. (1999) have also proposed of a disappearance of the inner disk during soft dips from BeppoSAX observations of GRS 1915+105. Note that the viscous time scales of the inner disk and that of the halo have almost same value of about 1 s which agrees with the fall time scale of the soft dips (Belloni et al. (2000); Yadav et al. (1999)). As the halo disappears more disk is observable reducing the Rin during the soft dips (see Table 2). Note here that the “spike” marks the beginning of the burst phase as described earlier. The accretion rate is reduced from $``$ 6.94 $`\times `$ 1018 g s-1 during the burst phase to $``$ 1.89 $`\times `$ 1018 g s-1 during the soft dip phase suggesting $``$ 5 $`\times `$ 1018 g s-1 mass is being blown away during these soft dips. In the above scenario, the missing energy of X-rays in the 5$``$100 keV range during the major soft dips is $``$ 1.4 $`\times `$ 1039 erg s-1 using the known distance of the source. Fender et al. (1999) have estimated a minimum power of $``$ 2 $`\times `$ 1039 erg s-1 is required for 12 hours of rising time of the relativistic radio jets observed on 1997 October/November on the basis of equipartition of energy. Including one proton for each electron approximately doubles this energy and requires a mass flow rate of $``$ 1018 g s-1. The energy requirement is broadly in agreement with our estimate of missing X-ray energy (5$``$100 keV) of $``$ 1.4 $`\times `$ 1039 erg s-1 as X-rays below 5 keV should significantly increase this estimate. During the soft dip, the accretion rate drops by $``$ 5 $`\times `$ 1018 g s-1 from that during the burst phase which is in agreement with the required mass outflow rate.

The average radio flux is given in Table 1 for all X-ray observations discussed here. As HR1 increases the radio flux also increases. The peak amplitude of IR/radio flares also increases with HR1 during class $`\beta `$ observations (Table 1). Our spectral analysis has suggested that X-ray luminosity increases with HR1. In Figure 4 we show the radio flux at 2.25 GHz from the NSF-NRAO-NASA Green Bank Interferometer for 1997 October 30-31 (data for $``$ 48 hours) on the y-axis (left). The X-ray average color HR1 during the quiescent phase is plotted on the y-axis (right) during the same time (data for $``$ 24 hours). The timing of our X-ray data falls between the first and second relativistic radio jets observed by Fender et al. (1999). Clearly, results shown in Figure 4 suggest a correlation between the radio and X-ray data. A least-square fit to the HR1 and radio data results a relation $`radio_{Jy}=0.615\times HR_10.322`$ with correlation coefficient close to one. This linear relation is based on a limited set of data and is valid in a narrow range of HR1 from 0.94 to 1.2 during the relativistic radio jets. The issues like extrapolation of this relation beyond the valid range of HR1, its applicability during the baby jets and fitting a non-linear relation (instead of a linear relation) will be discussed separately (Yadav (2000)). In Figure 5, we show the decay profile of this radio flare observed during 1997 October 30 - November 6 (offset by 50748.2 MJD). The profile is consistent with an exponential decay having a time constant of 3.157 day (dotted line). We have also plotted the decay profile of another radio flare of class A observed on 1999 June 8-15 (offset by 51334.5 MJD) in Figure 5. This flare has been observed in the X-ray band by IXAE/PPC (Naik et al. (2000)) as well as partially by RXTE/PCA. The decay profiles of these two flares of class A ($``$ two years apart) are identical except some deviation in the middle (on days 6 through 8). The average X-ray color HR1 during the low hard state is also plotted in Figure 5 for both the radio flares using the linear relation between the HR1 and radio flux density discribed above. These points are in good agreement with the decay profiles over the duration of observations when the radio flux drops by a factor of $``$ 2.

Let us summarise the physical picture emerging here. When the source changes state from a quiescent phase to the burst phase, the soft X-ray luminosity increases by a large factor ($``$ 10). If the total X-ray luminosity is above a certain threshold, it promptly triggers a major ejection event of synchrotron - emitting plasma producing a large soft dip after a spike in the X-ray light curve. Subsequent short but frequent dips and the undetectable continuous ejection during the burst phase, and probably also during the quiescent phase (if HR1 is high), produce an enhanced level of quasi - steady emission. This conclusion agrees with the observations of faint IR flares during the hard state which start even before the X-ray oscillations (Eikenberry et al. (2000)). On 1997 October 30, the total 5$``$100 keV X-ray flux during the quiescent phase was similar to that found during a “plateau” state ($``$ 1.7 $`\times `$ 10-8 erg s-1 cm-2) with the enhanced quasi-steady radio emission of $``$ 50 mJy (15 GHz) lasting for 20 days (Fender et al. (1999)) which would suggest mass ejection during the quiescent phase. It is also consistent with the value of k = 2 (suggesting continuous jets) derived from the flux densities of approaching and receding components of the relativistic radio jets observed on 1997 October/November (Fender et al. (1999)).

The radio emission during “baby jets” is consistent with synchrotron emission from an adiabatically expanding small cloud with a flat radio spectrum (Fender et al. (1997); Eikenberry et al. (1998)). A magnetic field of 8-16 G and total energy $``$ 1040 erg for each half hour flare are estimated assuming equipartition. The dominant decay mechanism is adiabatic expansion losses. For a pure adiabatic expansion, the synchrotron radio luminosity Iν(r) $``$ r-2p, where r is the distance (eq. (19.34) of Longair (1994)). This changes to Iν(r) $``$ r-2 for a flat radio spectrum. For a cloud of size $``$ 10 AU (minor axis at 15 GHz) and a relativistic speed of $``$ 10 AU hr-1 (Dhawan et al. (2000)), a unit peak radio luminosity will decay to a value of $``$ 0.44 in 15 minutes time which is in good agreement with the observations (Eikenberry et al. (1998); Mirabel et al. (1998); Pooley & Fender (1997)).

The radio emission during the relativistic jets is consistent with synchrotron emission from a extended radio cloud with ballistic motion on the scale of 300-5000 AU and steepening of the spectrum with time ($`\alpha `$ changes from 0.5 to 1.0 in a few days). Atoyan & Aharonian (1999) have shown that a single population of relativistic particles accelerated at the time of ejection cannot explain the radio emission from the relativistic jets and it requires continuous replenishment of energetic particles with energy dependent losses. Recently, Kaiser et al. (2000) have developed an internal shock model for the origin of relativistic radio jets in microquasars assuming quasi-continuous jet ejection. In this model, much of the energy needed to produce the radio emission is stored in the material of the continuous jet which was ejected by the central source prior to the formation of the shock fronts. The shock is caused by the collision of shells of jet material moving at different velocities. The formation of the internal shock ‘lights up’ the relativistic jets by accelerating particles which emit the observed synchrotron radiation. The relativistic jet mode in which strong internal shocks are produced, is thought to have strong variations in the jet speed. This is consistent with the observed X-ray color HR1 $``$ 1.1 in the initial phase of both the relativistic radio jets (Figure 5), which drops to a value of 0.94 in less than a day’s time (Figure 4). On the other hand, “baby jets” are supposed to represent a relatively stable mode where weak internal shocks (if any) do not play any significant role in the synchrotron emission. This is in agreement with observed values of X-ray color HR1 which remains in the range of 0.9$``$1.0 for both the data sets of 1997 August 14 and September 9 (Figure 3).

## 4. Conclusions

We present evidence of a direct accretion disk $``$ jet connection for relativistic radio jets. We have provided here an explanation for the most complex light curves with a “spike” observed to date in the Galactic microquasar GRS 1915+105. We suggest that the “spike” which separates the dips with hard and soft spectra marks the beginning of the burst phase, when the luminosity of the soft X-rays increases by a large factor ($``$ 10). This produces a major ejection episode of the synchrotron - emitting plasma termed as “baby jets” with the half-hour spacing widely reported in the literature. Subsequent short but frequent soft dips produce overlapping faint flares which result in an enhanced level of quasi-steady emission. The radio emission during “baby jets” is consistent with synchrotron emission from an adiabatically expanding small cloud (within a distance of few tens of AU of the inner accretion disk) with a flat radio spectrum and a decay time of several minutes. On the other hand, the radio emission during the relativistic radio jets is consistent with synchrotron emission from the extended radio cloud with ballistic motion on a scale of a few hundreds to few thousands AU and a decay time of a few days. These jets require continuous replenishment of energetic particles. Our results support the relativistic jet model with quasi-continuous mass ejection and strong variations in the jet speeds required for the formation of strong internal shocks to continuously accelerate the particles.

Author thanks the referee for his very useful comments which has greatly improved the contents and the presentation of the paper. He thanks A. R. Rao and K. P. Singh for their valuable comments on the manuscript. He also thanks RXTE and NSF-NRAO-NASA Green Bank Interferometer teams for making their data publicly available. The Green Bank Interferometer is a facility of the National Science Foundation operated by the NRAO in support of NASA High Energy Astrophysics programs.


text:
# References

Quantum information processing, science of111Article by E. H. Knill and M. A. Nielsen accepted for Supplement III, Encyclopaedia of Mathematics (publication expected Summer 2001). See also “http://www.wkap.nl/series.htm/ENM”. E. H. Knill is with the Los Alamos National Laboratory, MS B265, Los Alamos NM 87545, USA, and M. A. Nielsen is with the Center for Quantum Computer Technology, Department of Physics, University of Queensland 4072, Australia.  \- The theoretical, experimental and technological areas covering the use of quantum mechanics for communication and computation. Quantum information processing includes investigations in quantum information theory, quantum communication, quantum computation, quantum algorithms and their complexity, and quantum control. The science of quantum information processing is a highly interdisciplinary field. In the context of mathematics it is stimulating research in pure mathematics (e.g. coding theory, $``$-algebras, quantum topology) as well as requiring and providing many opportunities for applied mathematics.

The science of quantum information processing emerged from the recognition that usable notions of information need to be physically implementable. In the 1960s and 1970s researchers such as R. Landauer, C. Bennett, C. Helstrom and A. Holevo realized that the laws of physics give rise to fundamental constraints on the ability to implement and manipulate information. Landauer repeatedly stated that “information is physical”, providing impetus to the idea that it should be possible to found theories of information on the laws of physics. This is in contrast to the introspective approach which led to the basic definitions of computer science and information theory as formulated by A. Church, A. Turing, C. Shannon and others in the first half of the 20th century.

Early work in studying the physical foundations of information focused on the effects of energy limitations and the need for dissipating heat in computation and communication. Beginning with S. Wiesner’s work on applications of quantum mechanics to cryptography in the late 1960s, it was realized that there may be intrinsic advantages to using quantum physics in information processing. Quantum cryptography and quantum communication in general were soon established as interesting and non-trivial extensions of classical communication based on bits. That quantum mechanics may be used to improve the efficiency of algorithms was first realized when attempts at simulating quantum mechanical systems resulted in exponentially complex algorithms compared to the physical resources associated with the system simulated. In the 1980s, P. Benioff and R. Feynman introduced the idea of a quantum computer for efficiently implementing quantum physics simulations. Models of quantum computers were developed by D. Deutsch, leading to the formulation of artificial problems that could be solved more efficiently by quantum than by classical computers. The advantages of quantum computers became widely recognized when P. Shor (1994) discovered that they can be used to efficiently factor large numbers — a problem believed to be hard for classical deterministic or probabilistic computation and whose difficulty underlies the security of widely used public key encryption methods. Subsequent work established principles of quantum error-correction to ensure that quantum information processing was robustly implementable. See  for introductions to quantum information processing and a quantum mechanics tutorial.

In the context of quantum information theory, information in the sense of C. Shannon is referred to as *classical* information. The fundamental unit of classical information is the *bit*, which can be understood as an ideal system in one of two states or configurations, usually denoted by $`0`$ and $`1`$. The fundamental units of quantum information are qubits (short for “quantum bits”), whose states are identified with all “unit superpositions” of the classical states. It is common practice to use the bra-ket conventions for denoting states. In these conventions, the classical configurations are denoted by $`|0`$ and $`|1`$, and superpositions are formal sums $`\alpha |0+\beta |1`$, where $`\alpha `$ and $`\beta `$ are complex numbers satisfying $`|\alpha |^2+|\beta |^2=1`$. The states $`|0`$ and $`|1`$ represent a standard orthonormal basis of a two-dimensional Hilbert space. Their superpositions are unit vectors in this space. The state space associated with $`n>1`$ qubits is formally the tensor product of the Hilbert spaces of each qubit. This state space can also be obtained as an extension of the state space of $`n`$ classical bits by identifying the classical configurations with a standard orthonormal basis of a $`2^n`$ dimensional Hilbert space.

Access to qubit states is based on the postulates of quantum mechanics with the additional restriction that they are *local* in the sense that elementary operations apply to one or two qubits at a time. Most operations can be expressed in terms of standard measurements of a qubit and two-qubit *quantum gates*. The standard qubit measurement has the effect of randomly projecting the state of the qubit onto one of its classical states; this state is an output of the measurement (accessible for use in a classical computer if desired). For example, using the tensor product representation of the state space of several qubits, a measurement of the first qubit is associated with the two projection operators $`P_0^{(1)}=P_0I\mathrm{}`$ and $`1P_0^{(1)}`$, where $`P_0|0=|0`$ and $`P_0|1=0`$. If $`\psi `$ is the initial state of the qubits, then the measurement outcome is $`0`$ with probability $`p_0=P_0\psi ^2`$, in which case the new state is $`P_0\psi /p_0`$, and the outcome is $`1`$ with probability $`1p_0=P_1\psi ^2`$ with new state $`P_1\psi /(1p_0)`$. This is a special case of a *von Neumann measurement*. A general two-qubit quantum gate is associated with a unitary operator $`U`$ acting on the state space of two qubits. Thus $`U`$ may be represented by a $`4\times 4`$ unitary matrix in the standard basis of two qubits. The quantum gate may be *applied* to any two chosen qubits. For example, if the state of $`n`$ qubits is $`\psi `$ and the gate is applied to the first two qubits, then the new state is given by $`(UI\mathrm{})\psi `$. Another important operation of quantum information processing is preparation of the $`|0`$ state of a qubit, which can be implemented in terms of a measurement and subsequent applications of a gate depending on the outcome.

Most problems of theoretical quantum information processing can be cast in terms of the elementary operations above, restrictions on how they can be used and an accounting of the physical *resources* or *cost* associated with implementing the operations. Since classical information processing may be viewed as a special case of quantum information processing, problems of classical information theory and computation are generalized and greatly enriched by the availability of quantum superpositions. The two main problem areas of theoretical quantum information processing are quantum computation and quantum communication.

In studies of quantum computation (cf. quantum computation) one investigates how the availability of qubits can be used to improve the efficiency of algorithmic problem solving. Resources counted include the number of quantum gates applied and the number of qubits accessed. This can be done by defining and investigating various types of quantum automata, most prominently quantum Turing machines, and studying their behavior using approaches borrowed from the classical theory of automata and languages. It is convenient to combine classical and quantum automata, for example by allowing a classical computer access to qubits as defined above, and then investigating the complexity of algorithms by counting both classical and quantum resources, thus obtaining trade-offs between the two.

Most of the complexity classes for classical computation have analogues for quantum computation, and an important research area is concerned with establishing relationships between these complexity classes. Corresponding to the classical class $`𝐏`$ of polynomially decidable languages is the class of languages decidable in bounded error quantum polynomial time, $`\mathrm{𝐁𝐐𝐏}`$. While it is believed that $`𝐏`$ is properly contained in $`\mathrm{𝐁𝐐𝐏}`$, whether this is so is at present an open problem. $`\mathrm{𝐁𝐐𝐏}`$ is known to be contained in the class $`𝐏^{\text{\#}𝐏}`$ (languages decidable in classical polynomial time given access to an oracle for computing the permanent of $`0`$-$`1`$ matrices), but the relationship of $`\mathrm{𝐁𝐐𝐏}`$ to the important class of nondeterministic polynomial time languages $`\mathrm{𝐍𝐏}`$ is not known.

In quantum communication one considers the situation where two or more entities with access to local qubits can make use of both classical and quantum (communication) channels for exchanging information. The basic operations now include the ability to send classical bits and the ability to send quantum bits. There are two main areas of investigation in quantum communication. The first aims at determining the advantages of quantum communication for solving classically posed communication problems with applications to cryptography and to distributed computation. The second is concerned with establishing relationships between different types of communication resources, particularly with respect to noisy quantum channels, thus generalizing classical communication theory.

Early investigations of quantum channels focused on using them for transmitting classical information by encoding a source of information (cf. information, source of) with uses of a quantum channel (cf. quantum communication channel). The central result of these investigations is A. Holevo’s bound (1973) on the amount of classical information that can be conveyed through a quantum channel. Asymptotic achievability of the bound (using block coding of the information source) was shown in the closing years of the twentieth century. With some technical caveats, the bound and its achievability form a quantum information-theoretic analogue of Shannon’s capacity theorem for classical communication channels.

Quantum cryptography, distributed quantum computation and quantum memory require transmitting (or storing) quantum states. As a result it is of great interest to understand how one can communicate quantum information through quantum channels. In this case, the source of information is replaced by a source of quantum states, which are to be transmitted through the channel with high *fidelity*. As in the classical case, the state is encoded before transmission and decoded afterwards. There are many measures of fidelity which may be used to evaluate the quality of the transmission protocol. They are chosen so that a good fidelity value implies that with high probability, quantum information processing tasks behave the same using the original or the transmitted states. A commonly used fidelity measure is the Bures-Uhlmann fidelity, which is an extension of the Hilbert space norm to probability distributions of states (represented by *density operators*). In most cases, asymptotic properties of quantum channels do not depend on the details of the fidelity measure adopted.

To improve the reliability of transmission over a noisy quantum channel, one uses *quantum error-correcting codes* to encode a state generated by the quantum information source with multiple uses of the channel. The theory of quantum codes can be viewed as an extension of classical coding theory. Concepts such as minimum distance and its relationship to error-correction generalize to quantum codes. Many results from the classical theory, including some linear programming upper bounds and the Gilbert-Varshamov lower bounds on the achievable rates of classical codes have their analogues for quantum codes. In the classical theory, linear codes are particularly useful and play a special role. In the quantum theory, this role is played by the *stabilizer* or *additive* quantum codes, which are in one-to-one correspondence with self-dual (with respect to a specific symplectic inner product) classical $`\mathrm{GF}_2`$-linear codes over $`\mathrm{GF}_4`$ (cf. finite fields).

The capacity of a quantum channel with respect to encoding with quantum codes is not as well understood as the capacity for transmission of classical information. The exact capacity is known only for a few special classes of quantum channels. Although there are information theoretic upper bounds, they depend on the number of channel instances, and whether or not they can be achieved is an open problem. A further complication is that the capacity of quantum channels depends on whether one-way or two-way classical communication may be used to restore the transmitted quantum information .

The above examples illustrate the fact that there are many different types of information utilized in quantum information theory, making it a richer subject than classical information theory. Another physical resource whose properties appear to be best described by information-theoretic means is quantum entanglement. A quantum state of more than one quantum system (e.g. two qubits) is said to be entangled if the state can not be factorized as a product of states of the individual quantum systems. Entanglement is believed to play a crucial role in quantum information processing, as demonstrated by its enabling role in effects such as quantum key distribution, superdense coding, quantum teleportation, and quantum error-correction. Beginning in 1995 an enormous amount of effort has been devoted to understanding the principles governing the behavior of entanglement. This has resulted in the discovery of connections between quantum entanglement and classical information theory, the theory of positive maps  and majorization .

The investigation of quantum channel capacity, entanglement, and many other areas of quantum information processing involves various quantum generalizations of the notion of entropy, most notably the von Neumann entropy. The von Neumann entropy is defined as $`H(\rho )=\text{tr}\rho \mathrm{log}_2(\rho )`$ for density operators $`\rho `$ ($`\rho `$ is positive Hermitian and of trace $`1`$). It has many (but not all) of the properties of the classical information function $`H()`$ (cf. information, amount of). Understanding these properties has been crucial to the development of quantum information processing (see  for reviews). Probably the most powerful known result about the von Neumann entropy is the strong subadditivity inequality. Many of the bounds on quantum communication follow as easy corollaries of strong subadditivity. Whether still more powerful entropic inequalities exist is not known.

An important property of both classical and quantum information is that although it is intended to be physically realizable, it is abstractly defined and therefore independent of the details of a physical realization. It is generally believed that qubits encapsulate everything that is finitely realizable using accessible physics. This belief implies that any information processing implemented by available physical systems using resources appropriate for those systems can be implemented as efficiently (with at most polynomial overhead) using qubits. It is noteworthy that there is presently no proof that information processing based on quantum field theory (cf. quantum field theory) is not more efficient than information processing with qubits. Furthermore, the as-yet unresolved problem of combining quantum mechanics with general relativity in a theory of quantum gravity prevents a fully satisfactory analysis of the information processing power afforded by fundamental physical laws.

Much effort in the science of quantum information processing is being expended on developing and testing the technology required for implementing it. An important task in this direction is to establish that quantum information processing can be implemented robustly in the presence of noise. At first it was believed that this was not possible. Arguments against the robustness of quantum information were based on the apparent relationship to analogue computation (due to the continuity of the amplitudes in the superpositions of configurations) and the fact that it seemed difficult to observe quantum superpositions in nature (due to the rapid loss of phase relationships called *decoherence*). However, the work on quantum error-correcting codes rapidly led to the realization that provided the physical noise behaves locally and is not too large, it is at least in principle possible to process quantum information fault tolerantly. Research in how to process quantum information reliably continues; the main problems is improving the estimates on the maximum amount of tolerable noise for general models of quantum noise and for the types of noise expected in specific physical systems. Other issues include the need to take into consideration restrictions imposed by possible architectures and interconnection networks.

There are many physical systems that can potentially be used for quantum information processing . An active area of investigation involves determining the general mathematical features of quantum mechanics required for implementing quantum information. More closely tied to existing experimental techniques are studies of specific physical systems. In the context of communication, optical systems are likely to play an important role, while for computation there are proposals for using electrons or nuclei in solid state, ions or atoms in electromagnetic traps, excitations of superconductive devices etc. In all of these, important theoretical issues arise. These issues include how to optimally use the available means for controlling the quantum systems (*quantum control*), how to best realize quantum information (possibly indirectly), what architectures can be implemented, how to translate abstract sequences of quantum gates to physical control actions, how to interface the system with optics for communication, refining the theoretical models for how the system is affected by noise and thermodynamic effects, and how to reduce the effects of noise.

E. H. Knill  M. A. Nielsen

AMS 2000 Subject Classification: 81P68, 68Q05