text:
# Defibrillation via the Elimination of Spiral Turbulence in a Model for Ventricular Fibrillation

## Abstract

Ventricular fibrillation (VF), the major reason behind sudden cardiac death, is turbulent cardiac electrical activity in which rapid, irregular disturbances in the spatiotemporal electrical activation of the heart makes it incapable of any concerted pumping action. Methods of controlling VF include electrical defibrillation as well as injected medication. Both methods yield results that are subject to chance. Electrical defibrillation, though widely used, involves subjecting the whole heart to massive, and often counterproductive, electrical shocks. We propose a defibrillation method that uses a very-low-amplitude shock (of order mV) applied for a brief duration (of order 100 ms) and over a coarse mesh of lines on our model ventricle. Our proposal is based on a detailed numerical study of a model for ventricular fibrillation.

Ventricular fibrillation (VF), the leading cause of sudden cardiac death, is responsible for about one out of every six deaths in the USA. In the absence of any attempt at defibrillation, VF leads to death in a few minutes. The importance of this problem can hardly be overemphasized, so it is not surprising that many studies of VF have been carried out in various mammalian hearts and mathematical models. Recent experiments have shown that, though VF is turbulent cardiac activity, it is associated with the formation, and subsequent breakup, of electrophysiological structures that emit spiral waves; these are referred to as rotors or simply spirals. Such structures and their breakup have also been obtained in a set of partial differential equations of the Fitzhugh-Nagumo type; this set has been proposed as a model for VF. We show that the spiral turbulence associated with VF in this model arises because of spatiotemporal chaos; this is similar to the spatiotemporal chaos in a related model for the catalysis of CO on Pt(110) which we have studied elsewhere. Our defibrillation scheme for the VF model of Ref. relies on the understanding we have developed for spatiotemporally chaotic states in the CO catalysis model.

Electrical defibrillation entails the application of electrical jolts to the fibrillating heart to make it start beating normally again. Such defibrillation works only about two-thirds of the time and often damages the heart in the process of reviving it. In external defibrillation, electrical shocks ($`5`$ kV) are applied across the patient’s chest; they depolarise all heart cells simultaneously and essentially reset the pacemaking nodes of the heart. Slightly lower voltages ($`600`$ V) suffice in open-heart conditions, i.e., when the shock is applied directly on the heart’s surface. In our study of the VF model we achieve defibrillation by using very-low-amplitude pulses (of order mV) applied for a brief duration (of order 100 ms) and over a coarse mesh of lines on our model ventricle. Thus, if our defibrillation scheme can be realised in an internal defibrillator, it would constitute a significant advance in these devices.

Before proceeding to our quantitative results it is useful to define the model for ventricular fibrillation. For simplicity we use the model for isotropic cardiac tissue; in this case the equations governing the excitability $`e`$ and recovery $`g`$ variables are

$$\begin{array}{ccc}e/t\hfill & =\hfill & ^2ef(e)g,\hfill \\ g/t\hfill & =\hfill & ϵ(e,g)(keg).\hfill \end{array}$$
(1)
The function $`f(e)`$, which specifies fast processes (e.g., the initiation of the action potential) is piecewise linear: $`f(e)=C_1e`$, for $`e<e_1`$, $`f(e)=C_2e+a`$, for $`e_1ee_2`$, and $`f(e)=C_3(e1)`$, for $`e>e_2`$. We use the physically appropriate parameter values given in Ref., namely, $`e_1=0.0026`$, $`e_2=0.837`$, $`C_1=20`$, $`C_2=3`$, $`C_3=15`$, $`a=0.06`$ and $`k=3`$. The function $`ϵ(e,g)`$, which determines the dynamics of the recovery variable, is $`ϵ(e,g)=ϵ_1`$ for $`e<e_2`$, $`ϵ(e,g)=ϵ_2`$ for $`e>e_2`$, and $`ϵ(e,g)=ϵ_3`$ for $`e<e_1`$ and $`g<g_1`$ with $`g_1=1.8`$, $`ϵ_1=1/75`$, $`ϵ_2=1.0`$, and $`ϵ_3=0.3`$ (which lies in the range $`2>ϵ_3>0.1`$ suggested in Ref.). We solve model (1) by using a forward-Euler integration scheme. We discretise our system on a grid of points in space with spacing $`\delta x=0.5`$ dimensionless units and use the standard five- and seven-point difference stencils for the Laplacians in spatial dimensions $`d=2`$ and 3, respectively. Our spatial grid consists of a square lattice with $`L\times L`$ points or cubic lattices with $`L\times L\times L_z`$ points; in our studies we have used values of $`L`$ ranging from $`128`$ to $`512`$ and $`2L_z16`$. Our time step is $`\delta t=0.022`$ dimensionless units. As in Ref., we define dimensioned time $`T`$ to be $`5`$ ms times dimensionless time and $`1`$ spatial unit to be $`1`$ mm, such that the period and wavelength of a spiral wave are approximately 120ms and 32.5 mm, respectively. The dimensioned value of the conductivity constant in model (1) is 2 cm2/s . Our initial condition is a broken plane wave which we allow to evolve into a state displaying spiral turbulence; we control this eventually by our defibrillation scheme. On the edges of our simulation region we use no-flux (Neumann) boundary conditions since the ventricles are electrically insulated from the atria. We impose these boundary conditions numerically by adding an extra layer of grid points on each side of our simulation grid and requiring the values of $`e`$ be equal pointwise to their values on the points one layer within the boundary.

Before constructing an efficient defibrillation scheme for VF in model (1) it is important to appreciate the following points: (a) Ventricular fibrillation in this model arises because the system evolves to a state in which large spirals break down . (b) This state is a long-lived transient whose lifetime $`\tau _L`$ increases rapidly with $`L`$, the linear size of our system (e.g., for $`d=2`$, $`\tau _L850`$ ms for $`L=100`$ whereas $`\tau _L3200`$ ms for $`L=128`$). This property of model (1) is in qualitative accord with the experimental finding that the hearts of small mammals are less prone to fibrillation than those of large mammals . For time $`t\tau _L`$, a quiescent state with $`e=g=0`$ is obtained. (c) In systems with $`L128`$, $`\tau _L`$ is sufficiently long that a nonequilibrium statistical steady state is established. This state displays spatiotemporal chaos. For example, we find that there are several positive Lyapunov exponents $`\lambda _i`$ (averages for $`\lambda _i`$ are performed for $`\tau _0<t<\tau _L`$, where $`\tau _0`$ is the time of decay for initial transients); the number of positive $`\lambda _i`$ increases with $`L`$ (e.g., the Kaplan-Yorke dimension $`D_{KY}`$  increases from 7 to $`\stackrel{>}{}`$ 35 as $`L`$ increases from 128 to 256). (d) Model (1) is akin to a model for the catalysis of CO on Pt(110) in so far as both show spiral breakup and spatiotemporal chaos . Our recent studies of the CO catalysis model  have shown that Neumann or no-flux boundary conditions tend to absorb spiral defects and, indeed, the spirals do not last for appreciable periods of time on small systems.

Given that the long transient which leads to VF in model (1) is spatiotemporally chaotic, we might guess that the fields $`e`$ and $`g`$ have to be controlled globally to achieve defibrillation. In fact, some earlier studies  of spiral breakup in models for ventricular fibrillation have used global control. Here we show that a judicious choice of control points (on a mesh specified below) leads to an efficient defibrillation scheme for model (1). For $`d=2`$, we divide our simulation domain (of size $`L\times L`$) into $`K^2`$ smaller blocks and choose the mesh size such that it effectively suppresses the formation of spirals. For $`d=3`$, we use the same control mesh but only on one of the square faces of the $`L\times L\times L_z`$ simulation box. In our defibrillation scheme we apply a pulse to the $`e`$ field on a mesh composed of lines of width $`3\delta x`$. A network of such lines is used to divide the region of simulation into square blocks whose length in each direction is fixed at a constant value $`L/K`$ for the duration of control. (The blocks adjacent to the boundaries can turn be out to be rectangular). The essential point here is that, if a pulse is applied to the $`e`$ field at all points along the mesh boundaries for a time $`\tau _c`$, then it effectively simulates Neumann boundary conditions (for the block bounded by the mesh) in so far as it absorbs spirals formed inside this block (just like Neumann boundary conditions absorb spiral defects in the CO catalysis model ). Note that $`\tau _c`$ is not large at all since the individual blocks into which the mesh divides our system are of a linear size $`L/K`$ which is so small that it does not sustain long, spatiotemporally chaotic transients. Nor does $`K`$, which is related to the mesh density, have to be very large since the transient lifetime, $`\tau _L`$, decreases rapidly with decreasing $`L`$. We find that, for $`d=2`$, $`L=128`$, $`K=2`$ and $`\tau _c=41.2`$ ms is required for defibrillation. In Fig. 2 we illustrate such defibrillation with $`\tau _c=44`$ ms. For $`d=2`$, $`L=512`$, $`K=8`$ and $`\tau _c=704`$ ms suffices. Finally we show that a slight modification of our defibrillation scheme also works for $`d=3`$ (Fig. 3).

The efficiency of our defibrillation scheme is fairly insensitive to the height of the pulse we apply to the $`e`$ field along our control mesh so long as this height is above a threshold. To obtain the value of $`e`$ in mV units we have scaled the peak amplitude of a spike in the $`e`$ field (which has an amplitude of 0.9 in dimensionless units) to be equal to 110 mV. The latter is a representative average value of the peak voltage of an electrical wave in the heart . With this voltage scaling, dimensioned excitability is computed as $`110/0.9122.22`$ mV times dimensionless excitability. We find, e.g., that, for $`L=128`$, the smallest pulse which yields defibrillation is 57.3 mV/ms for the parameter values we use; however, we have checked that even stronger pulses (e.g., 278.3 mV/ms) also lead to defibrillation. We use a capacitance density of 1 $`\mu `$F/cm2 , which then yields a current density of 57.3 $`\mu `$A/cm2. Note that this threshold is of the order of the smallest potential ($`22.18`$ $`\mu `$A/cm2) required to trigger an action-potential spike in model (1) with $`g=0`$; the nullclines for model (1), in the absence of the Laplacian, are such that this smallest potential increases with increasing $`g`$. This is physically consistent with the increase in the refractory nature of the heart tissue with increasing $`g`$.

We have checked that (a) small, local deformations of our control mesh or (b) the angle of the mesh axes with the boundaries do not affect the efficiency of our defibrillation scheme. Furthermore the application of our control pulse on the control mesh does not lead to an instability of the quiescent state; thus it cannot inadvertently promote VF. We have checked specifically that, with $`e=0`$, $`g=0`$ as the initial condition at all spatial points, a wave of activation travels across the system when the control pulse is initiated; this travels quickly ($``$ 200 ms for $`L`$ = 128) to the boundary where it is absorbed and quiescence is restored. Our defibrillation method also works if $`g`$ is stimulated instead of $`e`$. This can be implemented by pharmaceutical means in an actual heart.

Our two-dimensional defibrillation scheme above applies without any change to thin slices of cardiac tissue. However, it is important to investigate whether it can be extended to three dimensions which is clearly required for real ventricles. A naive extension of our mesh into a cubic array of sheets will, of course, succeed in achieving defibrillation. However, such an array of control sheets cannot be easily implanted in a ventricle. We have tried to see, therefore, if we can control the turbulence in a three-dimensional version of model (1) on a $`L\times L\times L_z`$ domain but with the control mesh present only on one $`L\times L`$ face. For the open faces we use open boundary conditions and for the other faces we use no-flux Neumann boundary conditions. Our defibrillation scheme works if $`L_z4`$ (we have checked explicitly for $`L=220`$) but not for $`L_z>4`$. A slight modification of this scheme is effective even for $`L_z>4`$: Instead of applying a pulse for a duration $`\tau _c`$, we apply a sequence of $`n`$ pulses separated by a time $`\tau _{ip}`$ and each of duration $`\tau _w`$. We find that, if $`L=256`$ and $`L_z=8`$, defibrillation occurs in $`2002`$ ms with $`\tau _{ip}=22`$ ms, $`\tau _w=0.11`$ ms, $`n=30`$ and a control pulse amplitude of 57.3 $`\mu `$A/cm2; if $`L=128`$ and $`L_z=16`$, defibrillation occurs in $`1760`$ ms with $`\tau _{ip}=22`$ ms, $`\tau _w=0.11`$ ms, $`n=15`$ and a control pulse amplitude of 57.3 $`\mu `$A/cm2. We further find that optimal defibrillation is obtained in our model if $`\tau _{ip}`$ is close to the absolute refractory period for model (1) without the Laplacian term. The efficacy of our control scheme in the three-dimensional case can be understood heuristically as follows: Control via a steady pulse does not work for $`L_z>4`$ since the propagation of this pulse in the $`z`$ direction (normal to our control mesh) is blocked once the medium in the interior of our simulation domain becomes refractory. However, if we use a sequence of short pulses separated by a time $`\tau _{ip}`$, then, provided $`\tau _{ip}`$ is long enough for the medium to recover its excitability, the control-pulse waves can propagate in the $`z`$ direction and lead to successful defibrillation.

Typical electrical defibrillation schemes use much higher voltages than in our study. Recent studies  have explored low-amplitude defibrillation methods in model systems. They constitute an advance over conventional methods, but lack some of the appealing features of our defibrillation scheme. For example, the scheme of Ref. works only when the slow variable (the analog of our $`g`$) is controlled; though this can be done, in principle, by pharmaceutical means, it is clearly less direct than control via electrical means. Reference  uses electrical defibrillation, but has been demonstrated to prevent only one spiral from breaking up, as opposed to the suppression of a spatiotemporally chaotic state with broken spirals by our defibrillation scheme. We have checked explicitly that, for the spatiotemporally chaotic state of model (1), a straightforward implementation of the defibrillation scheme of Ref. is ineffective. \[This scheme applies pulses to the fast variable ($`e`$ in model (1)) on a two-dimensional, discrete lattice of points.\] The control current density in Ref.  is comparable to our $`57.3\mu `$A/cm2, which is much lower than 139$`\mu `$A/cm2, the maximum value of the ionic current during depolarization in the Beeler-Reuter model.

We have checked that our defibrillation scheme is not sensitively model dependent. For example, we have used the same scheme to eliminate spatiotemporal chaos associated with spiral breakup in the model for the catalysis of CO on Pt(110) mentioned above. We have found recently that our defibrillation scheme also works for the biologically realistic Beeler-Reuter model  for VF.

In conclusion, then, we have developed an efficient method for defibrillation by the elimination of spiral turbulence in model (1). Our method has the attractive feature that it uses very-low-amplitude pulses, applied only for a short duration on a coarse control mesh of lines. And, to the best of our knowledge, our study is the only one which shows how to attain defibrillation by the control of spatiotemporal chaos and spiral turbulence in a model for VF in both two and three dimensions. We hope our work will stimulate experimental tests of the efficacy of our defibrillation method.

We thank CSIR (India) for support, SERC (IISc, Bangalore) for computational facilities, and A. Pumir and N.I. Subramanya for discussions.

Figure Legends

Figure 1. Pseudograyscale plots showing the evolution of spatiotemporal chaos through spiral breakup in model (1) for physical times $`T`$ between 880 ms and 1430 ms. The panels on the left show the excitability $`e`$ and those on the right the recovery $`g`$ for all points ($`x`$,$`y`$) on a two-dimensional $`L\times L`$ spatial grid with $`L=128`$. The initial condition used is described in the text.

Figure 2. Pseudograyscale plots of the $`e`$ field in $`d`$ = 2 for $`L`$ = 128 (top) and $`L`$ = 512 (bottom) illustrating defibrillation by our control of spiral breakup in model (1). The control mesh divides the domain into four equal squares for $`L`$ = 128 and, for $`L`$ = 512, into 49 equal squares, 28 equal rectangles along the edges and 4 small equal squares on the corners. For $`L`$ = 128 we apply 57.3 $`\mu `$A/cm2 from $`T`$ = 891 ms to $`T`$ = 935 ms and by $`T`$ $`\stackrel{>}{}`$ 1500 ms spatiotemporal chaos is all but eliminated ($`|e|,|g|10^{13}`$ at all grid points). For $`L`$ = 512, we apply 250 $`\mu `$A/cm2 from $`T`$ = 55 ms to $`T`$ = 759 ms to the $`e`$ field. By $`T`$ = 1650 ms, spatiotemporal chaos is all but eliminated ($`|e|,|g|10^7`$ at all grid points) (see http://theory1.physics.iisc.ernet.in/rahul/images/heart.mpg for an animated figure).

Figure 3. Isosurface ($`e=0.6`$) plots of initial states (left panels) and pseudograyscale plots of the $`e`$ field on the top square face (right panels) illustrating defibrillation by the control of spiral breakup in model (1) for $`d=3`$: (top) $`L=128`$ and $`L_z=16`$, with the control mesh dividing the top face into 4 equal squares; and (bottom) $`L=256`$ and $`L_z=8`$, with the control mesh dividing the top face into 64 equal squares. In both cases we apply pulses on the control mesh (top face only) of 57.3 $`\mu `$A/cm2 with $`\tau _{ip}=22`$ ms and $`\tau _w=0.11`$ ms (see text). In the former case, with $`n=15`$, spatiotemporal chaos is all but eliminated by $`T`$ = 1760 ms ($`|e|,|g|10^4`$ at all grid points); in the latter, with $`n=30`$, it is all but eliminated by $`T`$ = 2002 ms ($`|e|,|g|10^3`$ at all grid points).


format: text

text:
# Decoupled oscillations and resonances of three neutrinos in matter

## I Introduction

Recently new experimental results for neutrino appear rapidly . More and more theoretical articles are published attempting at fixing their mass and mixing angles . Solar neutrino and some other experiments relate to the mixing and oscillations of the neutrino in matter . For studying these data, theoretical methods for solving the propagation and oscillation equations of neutrino in matter are required. Many important results have been obtained for two flavors. An important fact, discovered more than ten years ago, is that there is resonance of the mixing angle in matter due to the interaction of charged current . Because of the increase of mixing parameters the situation is more complex in three flavors than two. An important advance was gained in the research of MSW oscillation in three flavors in 1986. It was found that the three flavor oscillations are approximately decoupled into two independent oscillations and resonances between two flavors . Many important results from the researches of two flavors oscillation can be used in the three flavor case . Many works have been done and papers are published. They have been used in analyzing the experimental results, analytically and numerically . However a new theoretical framework and analytical method may be helpful to this task. We have made a series of efforts in this respect. In a previous paper, the author gave a new theoretical framework to discuss three neutrino mixing in matter . Results are given analytically in explicit, symmetric, and very simple form. Now we apply them to the approximate solutions and discuss three flavor oscillations in matter.

In section II, the results obtained in the first paper  and will be used in the paper, are quoted.

In section III, using and only using the condition of $`\mathrm{\Delta }m_1^2\mathrm{\Delta }m_2^2,`$ a strict approximate mathematical method is developed to deal with the eigenvalue equation of mass in matter for three neutrinos. By the method, a cubic algebra equation is reduced into two of quadratic, in two regions around $`A\mathrm{\Delta }m_1^2`$ and $`A\mathrm{\Delta }m_2^2`$ respectively. Thus, we get four solutions at each region, but only three of them satisfies our condition imposed on eigenvalues used for reducing the equation. In this way, decoupled oscillations and resonances can be read out directly. In addition, the solution which does not take part in the resonance can be calculated to higher precision which will be very useful in derivating propagating equation in a followed paper.

In section IV, a consistent rigorous mathematical derivation is presented using our theoretical framework  to expose the behaviors of the mixing parameters in matter around $`A\mathrm{\Delta }m_1^2`$ and $`A\mathrm{\Delta }m_2^2`$ when $`\mathrm{\Delta }m_1^2\mathrm{\Delta }m_2^2.`$ Approximately decoupled oscillation and resonance behavior are obtained between $`U_{e1}^m`$ and $`U_{e2}^m`$ around the $`A\mathrm{\Delta }m_1^2,`$ and between $`U_{e2}^m`$ and $`U_{e3}^m`$ around $`A\mathrm{\Delta }m_2^2`$ respectively. This phenomena is the partner of the oscillation and resonance between cos$`\theta ^m`$ and sin$`\theta ^m`$ in two flavors . In three flavor, this had been discovered as the decoupled resonances of $`\theta _1^m`$ around $`\mathrm{\Delta }m_1^2`$ and $`\theta _3^m`$ around $`\mathrm{\Delta }m_2^2`$ by T.K. Kuo and his collaborator previously . However, there are some phenomena which are not discussed or proved explicitly before. In the lower resonance, the amplitude, $`\sqrt{U_{e1}^{m2}+U_{e2}^{m2}}\sqrt{\eta _1^2+\eta _2^2}<1.`$ It is determined by the vacuum parameters. However, in the higher resonance, $`\sqrt{U_{e2}^{m2}+U_{e3}^{m2}}1`$ which is independent on the vacuum mixing parameters. There is notable difference between two flavor and three flavor oscillation and resonance in matter, even if in the case of decoupled resonances in three flavors. When the vacuum mixing angles are not small, the new observation is meaningful.

In the end, in section V, the results are discussed and the conclusion is given.

## II Theoretical framework

We begin from the general equation of the eigenvalve problem and some of the rigorous analytical solutions obtained in a previous paper . All symbols and definitions used in that paper will be used in here without explanation.. The equation dominated the eigenvalues of three neutrinos in matter is.

$$\lambda ^3+a\lambda ^2+b\lambda +c=0$$
(1)
where

$$a=(A+\mathrm{\Delta }m_2^2)$$
(2)
$$b=(\mathrm{\Delta }m_1^2)^2+\left[\mathrm{\Delta }m_1^2\left(\eta _{11}\eta _{22}\right)+\mathrm{\Delta }m_2^2\left(\eta _{11}+\eta _{22}\right)\right]A$$
(3)
and

$$c=(\mathrm{\Delta }m_1^2)^2\mathrm{\Delta }m_2^2\mathrm{\Delta }m_1^2\left[\mathrm{\Delta }m_2^2\left(\eta _{11}\eta _{22}\right)\mathrm{\Delta }m_1^2\eta _{33}\right]A$$
(4)
The solutions for $`\nu _e`$ mixing matrix elements in matter are

$$U_{eu}^m=N^u\underset{i=1}{\overset{3}{}}\frac{\eta _i^2}{\lambda _u\mathrm{\Delta }_i}u=1,2,3$$
(5)
where $`N^u`$ is the normalizing constant.

$$N^u=\left[\underset{i=1}{\overset{3}{}}\left(\frac{\eta _i}{\lambda _u\mathrm{\Delta }_i}\right)^2\right]^{\frac{1}{2}}u=1,2,3$$
(6)
and $`\mathrm{\Delta }_is`$ are

$$\begin{array}{c}\mathrm{\Delta }_1=\frac{1}{2}\left(m_2^2m_1^2\right)=\mathrm{\Delta }m_1^2\hfill \\ \mathrm{\Delta }_2=\frac{1}{2}\left(m_2^2m_1^2\right)=\mathrm{\Delta }m_1^2\hfill \\ \mathrm{\Delta }_3=m_3^2\frac{1}{2}\left(m_2^2+m_1^2\right)=\mathrm{\Delta }m_2^2\hfill \end{array}$$
(7)
The mass eigenvalues in matter expressed by $`\lambda _u`$ are

$$M_u^2=\lambda _u+\frac{1}{2}\left(m_1^2+m_2^2\right)$$
(8)
Remembering that we have proved

$$\eta _1^2+\eta _2^2+\eta _3^2=1$$
(9)
it is easy to see

$$U_{e1}^{m2}+U_{e2}^{m2}+U_{e3}^{m2}=1$$
(10)
This relation is helpful to our task below. Now we make use of these formulae and the condition $`\mathrm{\Delta }m_1\mathrm{\Delta }m_2`$ to derive approximate solutions in which decoupled resonances emerged explicitly.

## III Approximate solutions of eigenvalues

### A Dimensionless Equation

Before we begin work, it is convenient to make the equation dimensionless in favor of estimating the quantitative order for each term in the Eq.(1)-Eq.(4). Taking $`\lambda ,\mathrm{\Delta }m_1^2`$ and $`\mathrm{\Delta }m_2^2`$ in the unit of $`A`$, we introduce the following variable transformation.

$$\overline{\lambda }=\frac{\lambda }{A},\overline{\mathrm{\Delta }m}_1^2=\frac{\mathrm{\Delta }m_1^2}{A}and\overline{\mathrm{\Delta }m}_2^2=\frac{\mathrm{\Delta }m_2^2}{A}$$
(11)
We get

$$\overline{\lambda }^3+\overline{a}\overline{\lambda }^2+\overline{b}\overline{\lambda }+\overline{c}=0$$
(12)
where

$$\overline{a}=(1+\overline{\mathrm{\Delta }m}_2^2)$$
(13)
$$\overline{b}=(\overline{\mathrm{\Delta }m}_1^2)^2+\left[\overline{\mathrm{\Delta }m}_1^2\left(\eta _{11}\eta _{22}\right)+\overline{\mathrm{\Delta }m}_2^2\left(\eta _{11}+\eta _{22}\right)\right]$$
(14)
and

$$\overline{c}=(\overline{\mathrm{\Delta }m}_1^2)^2\overline{\mathrm{\Delta }m}_2^2\overline{\mathrm{\Delta }m}_1^2\left[\overline{\mathrm{\Delta }m}_2^2\left(\eta _{11}\eta _{22}\right)\overline{\mathrm{\Delta }m}_1^2\eta _{33}\right]$$
(15)
Three flavor neutrino oscillations can be decoupled approximately into two independent oscillations between two flavors and resonances, if the neutrino masses have the hierarchical characteristic, or more exactly if $`\mathrm{\Delta }m_1^2\mathrm{\Delta }m_2^2,`$ where $`\mathrm{\Delta }m_1^2=\left(m_2^2m_1^2\right)/2`$ and $`\mathrm{\Delta }m_2^2=m_3^2\left(m_2^2+m_1^2\right)/2.`$ This fact have been observed by T.K.Kuo and his collaborator previously. Now, we derive it in this new theoretical framework.

### B Solutions around $`A\mathrm{\Delta }m_1^2`$

First, we search for the approximate solutions around the region $`A\mathrm{\Delta }m_1^2`$. In present case,

$$\overline{\mathrm{\Delta }m}_1^21and\overline{\mathrm{\Delta }m}_2^21$$
(16)
#### 1 Solution to the small eigenvalues $`\lambda \mathrm{\Delta }m_1^2.`$

We try to search for the solutions which satisfy $`\lambda \mathrm{\Delta }m_1^2`$. Then $`\overline{\lambda }\overline{\mathrm{\Delta }m}_1^21.`$ Neglecting all the terms which is of order 1 in Eq.(12)-Eq.(15), we obtain the following approximate equation,

$$\overline{\lambda }^2\left(\eta _{11}+\eta _{22}\right)\overline{\lambda }\left[(\overline{\mathrm{\Delta }m}_1^2)^2\overline{\mathrm{\Delta }m}_1^2\left(\eta _{11}\eta _{22}\right)\right]=0$$
(17)
There are two solutions for it

$$\overline{\lambda }_{}=\frac{\left(\eta _{11}+\eta _{22}\right)\sqrt{\left(\eta _{11}+\eta _{22}\right)^2+4\left[(\overline{\mathrm{\Delta }m}_1^2)^2\overline{\mathrm{\Delta }m}_1^2\left(\eta _{11}\eta _{22}\right)\right]}}{2}$$
(18)
Both of them satisfy the condition $`\overline{\lambda }\overline{\mathrm{\Delta }m}_1^2`$ $`1`$. Thus we obtain two solutions

$$\lambda _{1,2}=\frac{A\left(\eta _{11}+\eta _{22}\right)\sqrt{A^2\left(\eta _{11}+\eta _{22}\right)^2+4\left[(\mathrm{\Delta }m_1^2)^2A\mathrm{\Delta }m_1^2\left(\eta _{11}\eta _{22}\right)\right]}}{2}$$
(19)
#### 2 Solution to the large eigenvalue $`\lambda \mathrm{\Delta }m_2^2.`$

We try to search for the third solution as $`\lambda \mathrm{\Delta }m_2^2.`$ That is the solution satisfying $`\overline{\lambda }\overline{\mathrm{\Delta }m}_2^21.`$ We can neglect the terms whose orders are equal or smaller than $`\overline{\mathrm{\Delta }m}_2^2`$ in Eq.(12)-Eq.(15). We get an approximate equation as

$$\overline{\lambda }^2\left(1+\overline{\mathrm{\Delta }m}_2^2\right)\overline{\lambda }+\overline{\mathrm{\Delta }m}_2^2\left(\eta _{11}+\eta _{22}\right)=0$$
(20)
There are two solutions

$$\overline{\lambda }_\pm =\frac{1}{2}\left[1+\overline{\mathrm{\Delta }m}_2^2\pm \sqrt{\left(1+\overline{\mathrm{\Delta }m}_2^2\right)^24\overline{\mathrm{\Delta }m}_2^2\left(\eta _{11}+\eta _{22}\right)}\right]$$
(21)
but, only the solution $`\overline{\lambda }_+`$ satisfies $`\overline{\lambda }\overline{\mathrm{\Delta }m}_2^2`$ which is used for obtaining the reduced equation. Therefore we get a maximum solution

$$\lambda _3=\frac{1}{2}\left[A+\mathrm{\Delta }m_2^2+\left(A+\mathrm{\Delta }m_2^2\right)\sqrt{1\frac{4A\mathrm{\Delta }m_2^2\left(\eta _{11}+\eta _{22}\right)}{\left(A+\mathrm{\Delta }m_2^2\right)^2}}\right]$$
(22)
Remaining to the order $`A,`$ we get

$$\lambda _3\mathrm{\Delta }m_2^2+A\eta _{33}$$
(23)
### C Solutions around $`A\mathrm{\Delta }m_2^2.`$

When $`A\mathrm{\Delta }m_2^2`$

$$\overline{\mathrm{\Delta }m}_2^21and\overline{\mathrm{\Delta }m}_1^21$$
(24)
#### 1 Larger solution $`\lambda \mathrm{\Delta }m_2^2.`$

We try to search for the solutions it satisfy $`\overline{\lambda }\overline{\mathrm{\Delta }m}_2^2.`$ Neglecting all the terms having the order smaller than $`\left(\overline{\mathrm{\Delta }m}_2^2\right)^2`$, from Eq.(12)-Eq.(15), we obtain an approximate equation

$$\overline{\lambda }^2(1+\overline{\mathrm{\Delta }m}_2^2)\overline{\lambda }+\overline{\mathrm{\Delta }m}_2^2\left(\eta _{11}+\eta _{22}\right)\lambda =0$$
(25)
There are two solutions for it

$$\overline{\lambda }_{}=\frac{1}{2}(1+\overline{\mathrm{\Delta }m}_2^2)\frac{1}{2}\sqrt{(1+\overline{\mathrm{\Delta }m}_2^2)^24\overline{\mathrm{\Delta }m}_2^2\left(\eta _{11}+\eta _{22}\right)}$$
(26)
Both of them satisfy the condition $`\lambda \mathrm{\Delta }m_2^2`$. Thus we obtain

$$\lambda _{2,3}=\frac{1}{2}(A+\mathrm{\Delta }m_2^2)\frac{1}{2}\sqrt{(A+\mathrm{\Delta }m_2^2)^24A\mathrm{\Delta }m_2^2\left(\eta _{11}+\eta _{22}\right)}$$
(27)
#### 2 Smaller solution $`\lambda \mathrm{\Delta }m_1^2`$

We try to search for the third solution in $`\overline{\lambda }\overline{\mathrm{\Delta }m}_1^21.`$ In Eq.(12)-Eq.(15), neglecting the terms whose order are equal or smaller $`\left(\overline{\mathrm{\Delta }m}_1^2\right)^2,`$ we obtain an approximate equation

$$\left(\eta _{11}+\eta _{22}\right)\overline{\lambda }\overline{\mathrm{\Delta }m}_1^2\left(\eta _{11}\eta _{22}\right)=0$$
(28)
Its solution is

$$\overline{\lambda }=\overline{\mathrm{\Delta }m}_1^2\frac{\eta _{11}\eta _{22}}{\eta _{11}+\eta _{22}}$$
(29)
That is

$$\lambda _1=\mathrm{\Delta }m_1^2\frac{\eta _{11}\eta _{22}}{\eta _{11}+\eta _{22}}$$
(30)
In fact, in our present method, higher precision can be reached for this solution by neglecting only the terms in order $`\left(\overline{\mathrm{\Delta }m}_1^2\right)^3.`$ Now we have a twice algebra equation

$$\begin{array}{cc}(1+\overline{\mathrm{\Delta }m}_2^2)\overline{\lambda }^2\left[\overline{\mathrm{\Delta }m}_1^2\left(\eta _{11}\eta _{22}\right)+\overline{\mathrm{\Delta }m}_2^2\left(\eta _{11}+\eta _{22}\right)\right]\overline{\lambda }\hfill & \\ \text{ }\left[(\overline{\mathrm{\Delta }m}_1^2)^2\overline{\mathrm{\Delta }m}_1^2\left(\eta _{11}\eta _{22}\right)\right]\overline{\mathrm{\Delta }m}_2^2+\left(\overline{\mathrm{\Delta }m}_1^2\right)^2\eta _{33}\hfill & =0\hfill \end{array}$$
(31)
There are two solutions for this equation, but only one of them satisfies $`\overline{\lambda }\overline{\mathrm{\Delta }m}_1^2`$ which is used to reduce the equation. It is the smaller one. The solution is a minimum solution for $`A\mathrm{\Delta }m_2^2`$. When we discuss propagating equation in a followed paper, this higher precision solution is very useful, but we are not necessary write it here explicitly.

## IV Approximately decoupled resonance behaviors of $`U_{e1}^m`$, $`U_{e2}^m`$, $`U_{e3}^m`$

In section II, we have introduced a set of solutions, Eq.(5) for the $`U_{e,u}^m`$ from a previous paper , because we are interested only in the mixing associated with $`\nu _e`$ in matter. They are expressed as functions of the vacuum parameters and potential $`A.`$ Now we use the equation to discuss the approximately decoupled resonant behaviors of $`U_{e1}^m,`$ $`U_{e2}^m,`$ and $`U_{e3}^m.`$

### A Resonant form for $`\lambda _1`$ and $`\lambda _2`$ around $`A\mathrm{\Delta }m_1^2`$.

When $`A\mathrm{\Delta }m_1^2,`$ we have get three approximate solutions. They are expressed in Eq.(19) and Eq.(22). There is a possible resonance between eigenstates of $`\lambda _1`$ and $`\lambda _2`$ when the quantity in the square root in Eq.(19) takes the minimum value. It can be reached when the following condition of $`A=A_l`$ is satisfied:

$$A_l=2\mathrm{\Delta }m_1^2\frac{\eta _{11}\eta _{22}}{\left(\eta _{11}+\eta _{22}\right)^2}$$
(32)
Introducing

$$\rho _l=A\left(\eta _{11}+\eta _{22}\right)2\mathrm{\Delta }m_1^2\frac{\eta _{11}\eta _{22}}{\eta _{11}+\eta _{22}}and\delta _l=\mathrm{\Delta }m_1^2\frac{2\eta _1\eta _2}{\eta _{11}+\eta _{22}}$$
(33)
we can rewrite the solutions Eq.(19) as

$$\lambda _{1,2}=\frac{1}{2}\left[A\left(\eta _{11}+\eta _{22}\right)\sqrt{\rho _l^2+4\delta _l^2}\right]A\mathrm{\Delta }m_1^2$$
(34)
It is clear that $`4\delta _l^2>0`$ when $`\eta _1\eta _20.`$ There is a possible resonance between eigenstates of $`\lambda _1`$ and $`\lambda _2`$ when the condition $`A=A_l`$ is satisfied. Because $`\lambda _3\mathrm{\Delta }m_2^2\lambda _{1,2}\mathrm{\Delta }m_1^2,`$ this resonance is decoupled approximately from $`\lambda _3.`$

Taking the traditional representation and let $`U=U_2U_3U_1`$ , we have

$$\frac{\eta _{11}\eta _{22}}{\left(\eta _{11}+\eta _{22}\right)^2}=\frac{c_1^2c_3^2s_1^2c_3^2}{\left(c_1^2c_3^2+s_1^2c_3^2\right)^2}=\frac{\mathrm{cos}2\theta _1}{c_3^2}$$
(35)
Thus, the lower resonance condition, Eq.$`\left(32\right),`$ becomes

$$A=A_l=2\mathrm{\Delta }m_1^2\frac{\mathrm{cos}2\theta _1}{c_3^2}$$
(36)
It is the same as the result obtained by Kuo et al.

### B The behaviors of $`U_{e1}^m`$, $`U_{e2}^m`$ and $`U_{e3}^m`$ in the neighborhood of $`A\mathrm{\Delta }m_1^2.`$

We can use the equations of $`U_{e1}^m`$ and $`U_{e2}^m`$ and take the approximation as follows

$$U_{e1,2}^mN^{1,2}\left(\frac{\eta _1^2}{\lambda _{1,2}+\mathrm{\Delta }m_1^2}+\frac{\eta _2^2}{\lambda _{1,2}\mathrm{\Delta }m_1^2}\right)$$
(37)
Correspondingly, we have

$$N^{1,2}\left[\left(\frac{\eta _1}{\lambda _{1,2}+\mathrm{\Delta }m_1^2}\right)^2+\left(\frac{\eta _2}{\lambda _{1,2}\mathrm{\Delta }m_1^2}\right)^2\right]^{\frac{1}{2}}$$
(38)
Substituting the Eq.(38) into the Eq.(37), we obtain

$$U_{e1,2}^m\frac{\lambda _{1,2}\mathrm{\Delta }m_1^2}{\left|\lambda _{1,2}\mathrm{\Delta }m_1^2\right|}\frac{\alpha _{1,2}\sqrt{\eta _1^2+\eta _2^2}}{\sqrt{\alpha _{1,2}^2+\left(\eta _1^2+\eta _2^2\right)\delta _l^2}}$$
(39)
Where

$$\alpha _{1,2}=\sqrt{\eta _1^2+\eta _2^2}\left(\lambda _{1,2}\frac{\eta _1^2\eta _2^2}{\eta _1^2+\eta _2^2}\mathrm{\Delta }m_1^2\right)$$
(40)
Using Eq.(34) of $`\lambda _{1,2}`$, we get

$$\alpha _{1,2}=\sqrt{\eta _1^2+\eta _2^2}\left(\frac{1}{2}\rho _l\frac{1}{2}\sqrt{\rho _l^2+4\delta _l^2}\right)$$
(41)
Then

$$U_{e1,2}^m\frac{\lambda _{1,2}\mathrm{\Delta }m_1^2}{\left|\lambda _{1,2}\mathrm{\Delta }m_1^2\right|}\frac{\alpha _{1,2}}{\left|\alpha _{1,2}\right|}\sqrt{\frac{\eta _1^2+\eta _2^2}{2}\left(1\frac{\rho _l}{\sqrt{\rho _l^2+4\delta _l^2}}\right)}$$
(42)
The $`\eta _3^m`$ is not independent. It is determined by $`U_{e1}^m+U_{e2}^m+U_{e3}^m=1.`$ from it, we obtain

$$U_{e3}^m\eta _3$$
(43)
When the resonance condition $`\rho _l=0`$ is satisfied, we have

$$U_{e1}^m\frac{\lambda _{1,2}\mathrm{\Delta }m_1^2}{\left|\lambda _{1,2}\mathrm{\Delta }m_1^2\right|}\frac{\alpha _{1,2}}{\left|\alpha _{1,2}\right|}\sqrt{\frac{\eta _1^2+\eta _2^2}{2}}$$
(44)
### C Resonant form for $`\lambda _2`$ and $`\lambda _3`$ around $`A\mathrm{\Delta }m_2^2`$

We have obtained three solutions when $`A\mathrm{\Delta }m_2^2`$. They are expressed in Eq.(27) and Eq.(30). There is a possible resonance between eigenstates of $`\lambda _2`$ and $`\lambda _3`$ when the quantity in the square root takes a minimum value. It is easy to show that the minimum can be obtained when the following condition of $`A=A_h`$ is satisfied:

$$A_h=\mathrm{\Delta }m_2^2\left(\eta _{11}+\eta _{22}\eta _{33}\right)$$
(45)
Let

$$\rho _h=A\mathrm{\Delta }m_2^2\left(\eta _{11}+\eta _{22}\eta _{33}\right)and\delta _h=\mathrm{\Delta }m_2^2\sqrt{\eta _{11}+\eta _{22}}\eta _3$$
(46)
the $`\lambda _{2,3}`$ can be rewrite as

$$\lambda _{2,3}=\frac{1}{2}\left(A+\mathrm{\Delta }m_2^2\sqrt{\rho _h^2+4\delta _h^2}\right)$$
(47)
It is clear that $`4\delta _h^2>0`$ when $`\left(\eta _{11}+\eta _{22}\right)\eta _{33}0.`$ We obtain a possible resonance. Because $`\lambda _1\mathrm{\Delta }m_2^2\lambda _{2,3}`$, this resonance between eigenstates of $`\lambda _2`$ and $`\lambda _3`$ is decoupled approximately from $`\lambda _1.`$

When we take the traditional representation and let $`U=U_2U_3U_1,`$ we have

$$\eta _1=c_1c_3,\eta _2=s_1c_3and\eta _3=s_3$$
(48)
Then,

$$\eta _{11}+\eta _{22}\eta _{33}=c_1^2c_3^2+s_1^2c_3^2s_3^2=\mathrm{cos}2\theta _3$$
(49)
The higher resonance condition, Eq.(45), becomes

$$A_h=\mathrm{\Delta }m_2^2\mathrm{cos}2\theta _3+O\left(\mathrm{\Delta }m_1^2\right)$$
(50)
It is the same as the result obtained by Kuo et al. at large $`A`$ resonance.

### D The behaviors of $`U_{e1}^m`$$`U_{e2}^m`$ and $`U_{e3}^m`$ in the neighborhood of $`A\mathrm{\Delta }m_2^2`$.

When $`A\mathrm{\Delta }m_2^2,`$ the higher resonance, $`\lambda _1\lambda _{2,3}`$. We can use the equations of $`U_{e2}^m`$ and $`U_{e3}^m`$ and make the approximation as follows

$$U_{e2,3}^mN^{2,3}\left[\frac{\eta _1^2}{\lambda _{2,3}}+\frac{\eta _2^2}{\lambda _{2,3}}+\frac{\eta _3^2}{\lambda _{2,3}\mathrm{\Delta }m_2^2}\right]$$
(51)
Correspondingly, we have

$$N^{2,3}\left[\left(\frac{\eta _1}{\lambda _{2,3}}\right)^2+\left(\frac{\eta _2}{\lambda _{2,3}}\right)^2+\left(\frac{\eta _3}{\lambda _{2,3}\mathrm{\Delta }m_2^2}\right)^2\right]^{\frac{1}{2}}$$
(52)
Substituting the Eq.(52) into the Eq.(51), we obtain

$$U_{e2,3}^m\frac{\lambda _{2,3}\left(\lambda _{2,3}\mathrm{\Delta }m_2^2\right)}{\left|\lambda _{2,3}\left(\lambda _{2,3}\mathrm{\Delta }m_2^2\right)\right|}\frac{\alpha _{2,3}}{\sqrt{\alpha _{2,3}^2+\delta _h^2}}$$
(53)
Where

$$\alpha _{2,3}=\lambda _{2,3}\mathrm{\Delta }m_2^2\left(\eta _1^2+\eta _2^2\right)$$
(54)
Using Eq.(47) of $`\lambda _{2,3},`$ we get

$$\alpha _{2,3}=\frac{1}{2}\rho _h\frac{1}{2}\sqrt{\rho _h^2+4\delta _h^2}$$
(55)
Then

$$U_{e2,3}^m\frac{\lambda _{2,3}\left(\lambda _{2,3}\mathrm{\Delta }m_2^2\right)}{\left|\lambda _{2,3}\left(\lambda _{2,3}\mathrm{\Delta }m_2^2\right)\right|}\frac{\alpha _{2,3}}{\left|\alpha _{2,3}\right|}\sqrt{\frac{1}{2}\left(1\frac{\rho _h}{\sqrt{\rho _h^2+4\delta _h}}\right)}$$
(56)
The $`U_{e1}^m`$ is not independent. It is determined by $`U_{e1}^m+U_{e2,3}^m+U_{e3}^m=1.`$ From it, we obtain

$$U_{e1}^m0$$
(57)
When the resonance condition $`\rho _h=0`$ is satisfied, we have

$$U_{e2,3}^m\frac{\lambda _{2,3}\left(\lambda _{2,3}\mathrm{\Delta }m_2^2\right)}{\left|\lambda _{2,3}\left(\lambda _{2,3}\mathrm{\Delta }m_2^2\right)\right|}\frac{\alpha _{2,3}}{\left|\alpha _{2,3}\right|}\sqrt{\frac{1}{2}}$$
(58)
## V Discussion and conclusion

In situation of three neutrinos in matter, the mixing is generally different to two neutrinos. They do not occur in any plane spanned by two of the three eigenvectors $`\nu ^{m,u}`$ $`(u=1,2,3)`$. We can not use a simple mixing angle to describe their mixing as done in two neutrino situation. However in case of $`\mathrm{\Delta }m_1^2\mathrm{\Delta }m_2^2,`$ we have two approximately decoupled resonances. In each case the oscillation do occurs between two mass eigenvectors, but there is important difference from two neutrinos situation.

Around $`A\mathrm{\Delta }m_1^2`$, matter effect make $`U_{e1}^m`$ and $`U_{e2}^m`$ variation when A change. The variation has oscillation characteristic and resonance. Oscillation occurs between $`\nu ^{m,1}`$ and $`\nu ^{m,2}`$ but it has an amplitude $`\sqrt{U_{e1}^m+U_{e2}^m}=\sqrt{\eta _1^2+\eta _2^2}<1.`$ This amplitude is determined by the vacuum mixing parameters. In plane spanned by $`\nu ^{m,1}`$ and $`\nu ^{m,2},`$ there exist a mixing angle $`\beta _l`$ which varies with $`A`$. The $`U_{e1}^m=a\mathrm{cos}\beta _l`$ and $`U_{e2}^m=a\mathrm{sin}\beta _l`$ with $`a=\sqrt{\eta _1^2+\eta _2^2}<1`$. In this meaning, we have oscillation and resonance in a plane spanned by $`\nu ^{m,1}`$ and $`\nu ^{m,2}`$. It constitutes only a part of $`\nu _e.`$ The $`\nu ^{m,3}`$ contribution is nonzero and has a approximately constant mixing coefficient $`U_{e3}^m\nu ^{m,3}\eta _3\nu ^{m,3}.`$ Noticeably $`\nu ^{m,3}`$ enter only due to vacuum mixing but not matter effect. It is not oscillation with the potential $`A`$ change. When $`U_{e3}=\eta _3`$ is not small, the nature is important. It is different with two neutrino oscillation and resonance in matter.

Around $`A\mathrm{\Delta }m_2^2`$, the case is same. However, because of the small mass, the mixing contribution of $`\nu ^{m,1}`$ is negligible to $`U_{e2}^m`$ and $`U_{e3}^m`$ oscillation. There exist a angle $`\beta _h`$ which varies with $`A`$ in a plane spanned by the eigenvector $`\nu ^{m,2}`$ and $`\nu ^{m,3}.`$ The $`U_{e2}^m=\mathrm{cos}\beta _h`$ and $`U_{e3}^m=\mathrm{sin}\beta _h`$ with a amplitude 1. The oscillation and resonance occur when potential A change. In this case, $`U_{e1}^m0`$ whatever the values of $`\eta _i`$ is taken. Then the oscillation and resonance is more as one of two neutrinos.

We have used a theoretical framework developed in a previous paper which is convenient for dealing with the three flavor neutrino oscillations in matter. In condition of $`\mathrm{\Delta }m_1^2\mathrm{\Delta }m_2^2,`$ approximately decoupled resonance behaviors is studied and discussed in details. New interesting physical results are obtained. There are important differences in resonant phenomena between two flavors and three flavors even in the case of decoupled resonances, in particular when $`\eta _3`$ is large. In traditional mixing angle description, there is mixing order problem. In some case, it may mask the physical characteristic and lead to misunderstanding.

###### Acknowledgements.

I would like to express my sincere thanks to Prof. A. S. Hirsch, Prof. T. K. Kuo and the members of High Energy Theory Group, Physics Department, Purdue University where part of the work was completed, for the kind hospitality extended to me during my visit from January to April, 2000. Many fruitful discussions with Prof. T. K. Kuo are very helpful to this work. I would like to think Mr. Hai-Jun Pan and Mr. Tao Tu for their help in doing this work. The author is supported in part by the National Science Foundation, China.


text:
# FLIERs as stagnation knots from post-AGB winds with polar momentum deficiency

## 1 Introduction and model description

FLIERs in planetary nebulae were originally identified with the structures previously known as ansae in elliptical planetary nebulae (Aller, 1941; Balick et al. 1993). Their peculiar characteristics are now recognized in a much wider variety of PNe (e.g. Guerrero, Vázquez & López 1999; Goncalves et al. 2000). Their nature has resisted a consistent explanation to date. Initially, FLIERs were considered as pairs of knots located symmetrically with respect to the PN nucleus and characterised by outflow radial velocities of the order of 30-50 $`\mathrm{km}\mathrm{s}^1`$. Ionization gradients decrease outwards from the nebular core and in often ‘head-tail’ morphologies are observed (see Balick et al. 1998).

However, and as pointed out by López (2000), the concept of FLIERs has been used in recent times to encompass nearly any \[N II\] bright knot in the periphery of PN-shells found to be traveling either with nearly null or very high radial velocity. A single model can hardly account for all the properties observed. Symmetric ejecta, shocks, recombination fronts or simply dynamical instabilities may form in the nebular rim developing dense, low ionization knots and expanding with the rim. These effects have not always been distinguished.

Regarding the observed position and velocities, projection effects have to be taken into account as well. A number of symmetric FLIERs, like the outer pair in NGC 7009, appear to be nearly aligned with the plane of the sky. Although they show only a few $`\mathrm{km}\mathrm{s}^1`$ radial velocity, the deprojected speed is estimated by Reay & Atherton (1985) as 160 (D/2) km/s, with D the distance in kpc. In some sources with extended ansae or strings of symmetric knots a linear increase of speed with distance has been found with values rising up to $`630`$ $`\mathrm{km}\mathrm{s}^1`$(Bryce et al. 1997, Corradi et al. 1999, O’Connor et al. 2000). Obviously, FLIERs are only part of a much richer variety of dynamic phenomena in PNe.

Models have proliferated with varying degrees of success in explaining the observed properties of these structures. For example, ionization fronts (IF) in localized dense knots and collimated flows which produce fast knots ramming through the shell of the PN have been discussed by Balick et al. (1998). Dopita (1997) has discussed FLIERs in terms of shocks immersed in the strongly radiative PN environment.

FLIERs do not usually show a physical connection between themselves and the nebular core. This fact has to be accounted for by any collimation model. In this regard, Frank, Balick & Livio (1996) discussed the case of a cooling wind that avoids fast reexpansion via a collimation mechanism similar to the one described by Cantó et al. (1988). In this model the wind slides along the aspherical outer shock towards the axis forming a dense, narrow jet. Alternatively, Redman & Dyson (1999) have presented a model in which FLIERs represent recombination fronts (RF) in mass-loaded jets.

On a different approach, García-Segura & López (2000) have recently developed 3D magneto- hydrodynamic models which successfully reproduce ansae-type structures by considering a shocked stellar wind carrying a toroidal magnetic field (see Ró$`\dot{\mathrm{z}}`$yczka & Franco, 1996) with mass-loss rates $``$ 10-7 M yr-1. Above this limit, jet structures develop in their model. In the MHD model, and the others mentioned above, the material cooling to form FLIERs comes from the stellar wind or a mixture with material from the ambient (in the case of mass-loading).

In this paper we elaborate on an alternative hydrodynamical model for the formation of axial, high-velocity ansae, introduced by Steffen & López (2000). In this model, FLIERs are formed from external gas swept-up by a fast, low-density wind plowing into the ambient medium of the PN, the latter formed from the slow dense wind of the central star during its AGB phase. It is worth noting that the group of elliptical nebulae with typical FLIERS show very extended halos, such as NGC 3242 (Meaburn, López & Noriega-Crespo, 2000) and NGC 7009 (Moreno-Corral, de la Fuente & Gutierrez, 1998)

### 1.1 The stagnation knot model

The idea of a ”stagnation knot” was used for the first time to reproduce the large-scale structure of the giant envelope of the PN KjPn8 (Steffen & López, 1998). We have generalized the model in such a way that it is able to produce symmetric high-density knots of jets from an uncollimated stellar wind or even from precessing jets. In our model, we postulate that there is a fast low-density outflow from the central object of a planetary nebula with a deficiency of momentum along the symmetry axis of the nebula (as compared to the momentum flowing along higher latitudes, see Figure 1a). The deficiency of momentum flux along the axis may arise from an interacting binary system, as described by Soker & Rappaport (2000). Alternatively, precessing collimated outflows or jets produce conical outflows (Peter & Eichler, 1995, and Lim, 2001) which also have the required flow geometry.

The reduced momentum flux along the axis causes the bow-shock to advance more slowly near the axis and hence the shock becomes concave instead of convex (as seen from outside; see Figure 2). The ambient medium passing through the oblique regions of this section is then refracted towards the axis, instead of away from it (see the velocity vectors on the right of Figure 1; these have been drawn in the reference frame of the stagnation region of the bow-shock). This converging motion is contrary to what occurs in conventional bow-shocks where the shocked external gas has diverging streamlines. Note that the ”concave” shape refers to a usually spherical symmetry of the stellar wind, i.e. a somewhat flat-top or ”boxy” geometry may be enough to produce ”polar caps” similar to those seen in the ”Cat’s Eye” nebula (NGC 6543), even if it is not strictly ”concave”.

The reverse shock going into the fast stellar wind, has a similar shape, but causes the fast wind to diverge (see Figure 1). Consequently, the accumulated external material in the stagnation region is surrounded and confined by the high pressure of the surrounding hot gas from the diverging shocked stellar wind. The confinement allows sufficient time for the collected external gas to cool and be compressed to a dense knot or long jet-like feature. These are the structures we refer to as stagnation knots. The fast and dilute stellar wind diverges from the axis and is too hot to cool during the required timescales of order $`10^3`$ to $`10^4`$ years. This is an important difference between this model and previous models.

In many cases, the symmetric ansae are seen well outside and detached from the main envelope. In our model this occurs once the fast stellar wind stops or reduces its power, i.e. in the case of a non-steady outflow (see Figure 1b). Variable PN nuclear winds have been clearly shown to operate in cases like the Stingray Nebula (Bobrowsky et al., 1998), LMC-N66 (Peña et al., 1997) and Lo 4 (Werner et al. 1992). These few well documented cases indicate that variable winds in evolving PNe may be a common condition, rather than the exception. In our model the wind is considered to cease or decrease temporarily a few hundred or thousand years after the formation of the PN and the expansion of the envelope slows down. The stagnation knot has by now accumulated sufficient mass and, given its higher density and associated specific momentum, it slows down at a lower rate than the envelope and consequently may escape from it and propagates into the ambient medium. The evolution of such a dense knot propagating through a thin ambient medium has been studied in some detail by Jones, Kang & Tregillis (1994) using hydrodynamic simulations and has been also discussed by Soker & Regev (1998) in the context of FLIERs.

In this paper we shall concentrate on the dynamical and kinematical properties of the formation and initial evolution of the stagnation knots. The study of the ionization structure around the stagnation knot, taking into account the complex shock structures and the photo-ionization from the central star, is out of the scope of the present work and will be addressed in a future paper. We show here that at least three different types of structure can be produced: first, slow ansae which remain inside or near the rim of the main nebula; second, fast ansae which move a relatively large distance out of the rim; third, very long, jet-like strings with an internal linear velocity increase as a function of distance from the source.

## 2 Simulations

The 2D-hydrodynamical simulations in axisymmetry have been calculated using the Coral-code with a 5-level binary adaptive grid (Raga et al. 1995). This code solves the equations of mass, momentum, and energy conservation using a flux-vector-splitting scheme (van Leer 1982). The non-equilibrium cooling as described in Biro, Raga, & Cantó (1995) has been used. For low temperatures, energy loss from the collisional excitation of \[O I\] and \[O II\] has been taken into account. In order to simulate the lower limits on temperature imposed by photo-ionization from the central star (without explicit calculation of the photo-ionization) a lower limit of 5000 K has been kept for the temperature. The grid size at full resolution is $`513\times 257`$ grid cells. The full domain has a physical size of $`1\times 10^{18}`$ cm by $`5\times 10^{17}`$ cm. The outflow was initialized on a sphere with a radius of $`5\times 10^{16}`$ cm. We have used a Courant-number of 0.1 times the smallest signal crossing-times combined with the cooling time scale of all active cells on the adaptive grid.

The three-dimensional simulation has been produced with the Reefa-code described in Lim & Steffen (2001). This code uses a similar scheme as Coral (Raga et al. 1995) to solve the hydrodynamic equations on a binary adaptive 3D-grid. The grid size is variable and increases as the expanding dynamical region nears the boundary (at the time shown, the simulation in Figure 3d has dimensions of $`345\times 345\times 289`$ cells). Again we use the flux-vector splitting method due to Van Leer (1982), with the mass density conservation equation written separately for each of the ionic and neutral components. The radiative energy loss is computed with the prescription described by Biro et al. (1995). The collisional ionization of hydrogen by electron impact is included with the rate coefficients of Cox (1970) and recombination with the interpolation formula of Seaton (1960).

### 2.1 Initial and boundary conditions

The main observed dynamical and structural properties which served as a guideline for our simulations are the following: radius of the brightest sections of the nebula between 0.25 and 0.5 parsec, rim densities of order 1,000 cm-3 with expansion speeds near 50$`\mathrm{km}\mathrm{s}^1`$ at this size. The properties of the axial FLIERs should include a density of the order of $`10^4`$ cm-3 and speeds between that of the rim and around 200$`\mathrm{km}\mathrm{s}^1`$, which may rise linearly with distance up to velocities over 600$`\mathrm{km}\mathrm{s}^1`$ , as observed in the case of MyCn 18 (Bryce et al. 1997, O’Connor et al. 2000).

The initial density distribution $`n(\theta )`$ has been assumed to be static varying with polar angle $`\theta `$ according to the equation

$$n(\theta ,r)=n_0\left[(1q)\mathrm{sin}^\delta \theta +q\right](r/r_0)^\alpha .$$
(1)
This is equivalent to the one used by Kahn & West (1985), except for the free power-law parameter $`\alpha `$. It takes into account possible deviations from the inverse square law, which might arise from time-variations of the AGB-wind. Here $`n_0`$ is the density at the fiducial distance $`r_0`$ from the star, $`q`$ is the pole to equator density ratio, $`\delta `$ controls how fast the density rises from pole to equator. This density distribution in the environment fits calculations of stellar outflows taking into account slow stellar rotation (Reimers, Dorfi, & Höfner, 2000). In the 3D-simulation the density of environment has been assumed to be uniform.

The latitudinal momentum distribution in the fast wind will depend on the model adopted for its origin. Precessing jets are one possibility to produce the required momentum distribution for the formation of stagnation knots. On the other hand, rather uncollimated interacting winds from double stars are another possible mechanism to produce a deficiency of momentum perpendicular to the plane of the binary (Soker & Rappaport, 2000). However, at this time it is not clear whether this sort of axial reduction of momentum can hold up to the distances required for the formation of stagnation knots. Currently single stellar wind scenarios do not seem to be able to produce polar momentum deficiencies. Considerable simplifications in the theory of winds from rapidly rotating stars do, however, leave some open ground for winds with polar momentum deficiencies (Petrenz & Puls, 2000; S. Cranmer private communication). Hence, the details of the velocity distribution used in this paper are to be regarded only as one possibility of many. Simulations with other distributions show essentially the same qualitative results.

Since there is no simple analytic description of the wind momentum as a function of polar angle, for this initial qualitative study, we therefore adopt the velocity distribution $`v(\theta )`$ of the fast wind (equation 2) to follow the same function of polar angle as the density of the slow wind, with the difference that the maximum speed $`v_0`$ is reached at some half opening angle $`\theta _0`$. At higher polar angles than $`\theta _0`$ the wind has constant velocity. The numerical parameters $`\kappa `$ and $`ϵ`$ may of course have different values to $`q`$ and $`\delta `$ in equation 1.

$$\frac{v(\theta )}{v_0}=(1ϵ)\mathrm{sin}^\kappa \left(\frac{\pi \theta }{2\theta _0}\right)+ϵ$$
(2)
For the simulations we have used the parameters shown in Table 1.

## 3 Results and Discussion

We have performed a series of simulations varying the following parameters: the pole-equator ratios of the velocity of the fast flow, the densities, the opening-angle of the zone of momentum depletion, and the pole-equator density ratio in the ambient medium. Figure 1 shows representative results.

In the simulations we identify the three main structures which are observed in many PNe. First, a high density envelope or rim of shocked external gas which propagates at a few tens of kilometers per second (see also Figure 1). Second, small knots formed by instabilities in the envelope, which might correspond to low-ionization knots which are not aligned with the symmetry axis of the nebula (Dwarkadas & Balick 1998). Third, fast dense knots near the axis which we associate with the ansae and which are the subject of this paper.

As long as the central outflow is on, driving the expanding shock, the stagnation knot moves roughly at the same speed as the rest of the bow-shock and remains at the bright rim formed by shocked ambient gas. As soon as the outflow ceases, however, the cooling envelope slows down rather quickly, whereas the dense knot continues to move at its original speed. Later, it slows down as it expands and continues to add mass from the ambient medium. The deceleration of the expanding envelope is ever so faster than that of the stagnation knot, since the ambient density in the equatorial region is higher than that in the polar region.

Instabilities in the thin cold envelope (representing the main PN-shell) may develop. In these simulations the development of instabilities is susceptible to the details of the numerical procedure, e.g. the Courant number, resolution and diffusion of the code. The kinematic results on them are therefore not very reliable (see Dwarkadas and Balick, 1998, for a discussion of the instabilities). The stagnation knot is not noticeably sensitive to the numerical details, because it does not form on the scale of the grid resolution.

Figures 3a and 3b show stagnation knots of the slow (RUN 1) and the fast (RUN 2) types, respectively (see Table 1 for the parameters). The slow knot starts to expand significantly before it separates much from the rim of the nebula. Small deviations from the transverse viewing angle will make the knot appear to be inside the PN-shell. Its speed at the time of the density image shown in Figure 1a is around 140$`\mathrm{km}\mathrm{s}^1`$. The fast knot in Figure 3b is able to reach a greater distance before it is slowed down as a consequence of its expansion. The expansion causes an increase of the rate at which mass is swept up, and hence greater deceleration. The speed of the fast stagnation knot reaches 240$`\mathrm{km}\mathrm{s}^1`$ at the tip.

Inspection of the structure of the axial FLIERs in the simulations reveal different possible morphologies. We find head-tail structures directed both ways, towards and away from the source, as those observed in the HST images of NGC 3242, 6826 and 7009, presented by Balick et al. (1998). Tails directed towards the central star are, however, short-lived (see below). RUN 1 (Figure 3a) shows an extended blob, with some head-tail structure, with the head towards the source. RUN 2 is similar to the structures found in the Saturn-nebula (NGC 7009), where we find a crossing shock feature on the axis with the highest density on the near side of the star. This feature lies within the main rim. Outside the rim, there is the main stagnation knot, which has developed a bow-shock structure, with highest density near the tip. Lower density wings extend back towards the source.

The main differences in the conditions from which me obtain one or another type of stagnation knot is as follows. The slow knot is produced in the case where there is a strong reduction of momentum along the axis (by a factor of 0.2), a relatively high ambient density along the axis and a smaller wind speed. The reverse is true for the fast knot. Here the reduction of momentum is only moderate (by a factor of 0.5). Densities in both stagnation knots vary roughly between $`10^3`$ and $`10^4`$ cm-3.

Figure 3c (table 1 RUN 3) shows a very long stagnation knot or jet. Its nature is very different from that of jets from young stellar objects or active galaxies, which are collimated very near the central source and are refueled continuously or in pulses. In the stagnation knots, we have swept-up material which is not replaced from the central source. The most important parameters for the production of such a long jet are a small opening angle of the region of reduced momentum flow and a large pole-equator density ratio in the environment.

The axial jet-like structure initially has its highest density on the far-side as seen from the star, since this is where catastrophic cooling starts first. At that time, this feature would be seen as a head-tail structure with the head and lowest-ionization on the far side. This situation is, however, transient (about 100 years in this simulation). After this time, the reverse is true, i.e. the denser region with presumably lower ionization faces the star. The reason for this is that the far side collapses first and reexpands first as well, such that later it has lower density as compared to the region closer to the centre.

We find that in the stagnation ”jets”, the velocity increases linearly with distance (see inset of Figure 3c). This velocity structure is very similar to that seen in the magneto-hydrodynamic models of ansae developed by García-Segura et al. (1999). Both models, the stagnation knots and the MHD-model reproduce this important observational result. This problem deserves a separate study.

The three-dimensional calculations (Figure 3d, table 1 RUN 3D) show that the formation of stagnation knots is not a consequence of the symmetry imposed in the axisymmetric calculations. In 3D, the knots form and develop in a very similar fashion to 2D simulations.

From our simulations we can discuss the ionization state of the knots only as suggested by the density structure, assuming ionization balance with respect to incoming photons. Thus to first order we take that high density in the knots implies low ionization and vice-versa. A detailed ionization study will have to include not only photo-ionization from the central star, but also collisional and possibly local photo-ionization from strong and fast shocks, since these are ubiquituous in our model. These local sources of ionization may complicate the ionization structure considerably. Such a study is outside the scope of the present paper.

As for the geometry of the polar region predicted by the stagnation knot model, a comparison with observations has to focus on young PNe, since the expected concavity (in the extended definition proposed here) is only short-lived (of the order of a few hundred years). Hence, we expect to find the concavity only in a rather small fraction of objects in which, in addition, the axial FLIERs have not yet fully developed. The young PN IRAS 17150-3224 is such an example. It shows a pair of elongated, truncated lobes with concave ends (Kwok, Su & Hrivnak, 1998). An additional clear example of a PN showing a concave bow-shock structure is K 3-24 (Manchado et al. 1996). This type of structure is best obtained with a conical momentum distribution, i.e. with a clear maximum at some polar angle (Steffen & López, 2000). Physical mechanism for producing such a momentum distribution can be, e.g., a rapidly precessing jet (Peter & Eichler, 1995; Lim 2000) or a close binary system where one of the components blows a collimated fast wind (Soker & Rappaport, 2000).

## 4 SUMMARY

We have shown that the high-velocity axial FLIERs found in some bipolar/elliptical planetary nebulae could be due to a polar deficiency of the momentum flow of the fast stellar wind which interacts with the ambient medium. This wind configuration accumulates and compresses shocked ambient gas in the stagnation region of the concave (locally inverted) bow-shock. We find basically three different types of stagnation features. First, those which remain on or within the rim of the PN. Second, knots which escape the rim and, third, long jet-like features which may disintegrate to form a series of knots. The stagnation knots can be produced with speeds of up to several hundred kilometers per second. The stagnation jets show a constant positive velocity gradient with distance from the source. Both the magnitude and the gradients of the velocities are consistent with observations.

W.S. is grateful for travel support from project DGAPA-UNAM IN114199 and University College London visiting researchers programme. J.A.L. acknowledges support from project DGAPA-UNAM IN114199. A.J.L. acknowledges a PPARC research associateship. We thank Guillermo García-Segura and Alejandro Raga for useful discussions. We are also greatful to the constructive suggestions by the anonymous referee.


text:
# Gamma Ray and Hadron generated Čerenkov Photon Spectra at Various Observation Altitudes

## 1 Introduction

In the highest energy end of the electromagnetic spectrum, the energy range between 50 keV and 30 GeV has been successfully covered by satellite based missions. At higher the energies steeply falling spectrum makes it impossible to observe the $`\gamma `$ray sky by satellite based detectors. The ground based atmospheric Čerenkov technique has proved to be the most sensitive technique in exploring the celestial $`\gamma `$rays in the energy range 300 GeV - 10 TeV. The basic detection technique in this case is to deploy one or several parabolic mirrors at the observation level fitted with fast photo tubes at the focus. Very high energy celestial $`\gamma `$rays initiate an electromagnetic cascade in the atmosphere as they enter the atmosphere. The electrons and positrons in the cascade, being relativistic, emit Čerenkov light as they propagate down the atmosphere resulting in a faint ($`10^4`$ times the brightness of the star light background) flash of light lasting a few nanoseconds. This fast Čerenkov flash is detected electronically by coincidence technique. However the main drawback of the technique is the presence of much more abundant cosmic ray background which severely limit the sensitivity of this technique. Various ingenious techniques are employed to distinguish between the Čerenkov light produced by cosmic $`\gamma `$rays from that by the cosmic rays. Among these, the imaging technique, wavefront sampling technique, angular and spectral separation techniques are currently employed by different experiments .

In spite of the necessity of carrying out these observations during clear, moon-less nights which scales to only about 10-15% duty cycle, the technique has been successful mainly due to two reasons. Firstly, a large collection area which is $`10^410^5`$ times larger than the satellite detectors provides a high event rate. Secondly, with high angular resolution and some recently developed $`\gamma `$ hadron separation techniques  provide high sensitivity to the method. Large ground based telescopes based on this technique are now being built to cover the intermediate energy range of 30 - 300 GeV .

In the absence of standard mono-energetic beams of cosmic rays or $`\gamma `$rays, one has to depend on the simulation techniques to understand and optimize the detector response to these radiations. For this purpose we have carried out detailed simulations with VHE $`\gamma `$rays and cosmic rays of various primary energies initiating extensive air showers in the atmosphere. In order to ensure the reliability of the conclusions drawn from these studies it is imperative to verify the correctness of the basic results derived from these simulations. There are a few simulation packages painstakingly developed by high energy physicists and we chose one of them called CORSIKA  developed by the KASCADE collaboration.

CORSIKA (version 562),  has been used here to simulate Čerenkov light emission in the earth’s atmosphere by the secondaries of the extensive air showers generated by cosmic ray primaries or $`\gamma `$rays. This program simulates interactions of nuclei, hadrons, muons, electrons and photons as well as decays of unstable secondaries in the atmosphere. It uses EGS4 code  for the electromagnetic component of the air shower simulation and GHEISHA for the simulation of hadronic interactions at TeV energies. The Čerenkov radiation produced within the specified bandwidth (300-550 $`nm`$) by the charged secondaries is propagated to the ground. The US standard atmosphere parameterized by Linsley  has been used. The position, angle, time (with respect to the first interaction) and production height of each Čerenkov photon hitting the detector at the observation level are recorded.

We have mainly used Pachmarhi (longitude: 78 26 E, latitude: 22 28$`{}_{}{}^{}N`$ and altitude: 1075 $`m`$) as the observer’s location where an array of 25 Čerenkov detectors each of area 4.35 $`m^2`$ is deployed in the form of a rectangular array. We have assumed 17 detectors in the E-W direction with a separation of 25 $`m`$ and 21 detectors in the N-S direction with a separation of 20 $`m`$. This configuration, similar to the Pachmarhi Array of Čerenkov Telescopes (PACT)  but much larger, is chosen so that one can study the core distance dependence of various observable parameters. Simulations were also carried out at sea level and a location at an altitude of 2 $`km`$ above mean sea level. Mono-energetic primaries consisting of $`\gamma `$ rays and protons incident vertically on the top of the atmosphere with their cores at the centre of the array have been simulated in the present study.

The longitudinal shower development profiles of Čerenkov photons in the atmosphere have been generated using CORSIKA for $`\gamma `$ray and proton primaries of various energies. Tracks of charged particles are followed through in steps $`20gcm^2`$ .

In the present work the Čerenkov photon growth curves have been obtained from the simulations for proton and $`\gamma `$ray primaries of various energies and for the three observation levels mentioned above. These have been compared with those derived from independent empirical relations. The Čerenkov photon spectra as seen at the observation levels are then derived for various primary energies of $`\gamma `$rays and protons after taking into account of the wavelength dependent attenuation in the atmosphere. Then we estimate the relative fraction of UV photons in these spectra as a function of primary energy and observation level.

## 2 Calculation of longitudinal profiles

### 2.1 $`\gamma `$ray primaries

For showers initiated by photons of energy $`E_\gamma `$ the longitudinal shower development curve is derived from the equation for cascade curve under approximation A. The average number of electrons of all energies $`N(E_\gamma ,x)`$ as a function of atmospheric depth $`x`$, is given approximately, in the region where the number of particles is large, by :

$$N(E_\gamma ,x)=\frac{0.31}{\sqrt{\mathrm{ln}\left(\frac{E_\gamma }{E_t}\right)}}\mathrm{exp}(t(11.5\mathrm{ln}s))$$

where $`t`$ is the depth $`x`$ measured in units of radiation length in the atmosphere (37.2 $`gcm^2`$) and $`E_t`$ is the the electron critical energy at the depth $`t`$. This is the energy below which the electron multiplication stops. It is approximately the energy at which radiation losses and collisional losses are equal and has a value 84.2 MeV in air. The shower age parameter, $`s`$, which is a measure of the stage of the longitudinal shower development, is given by:

$$s=\frac{3t}{(t+2\mathrm{ln}\frac{E_\gamma }{E_t})}$$

$`s`$ = 1 at shower maximum while $`s<1`$ above the shower maximum and and $`s>1`$ below. The relation between the atmospheric height $`h`$, (as measured from sea level in $`m`$) and the depth, $`x`$ ($`gcm^2`$ as measured from the top of atmosphere) is taken from Rao:

$$h=(6740+2.5x)\mathrm{ln}\left(\frac{1030}{x}\right)m$$

where the scale height is dependent on the atmospheric depth. The refractive index $`n`$ at a given height $`h`$ is assumed to be:

$$n=1+n_0$$

where $`n_0`$ is given by:

$$n_0=2.9\times 10^4e^{\frac{h}{7100}}$$

For calculating the electron growth curve the atmosphere is divided into slabs of thickness 333 $`m`$, after the height of first interaction of the primary. The electron threshold energy for the production of Čerenkov radiation at a depth $`t`$ is given by

$$E_{th}=\left(\frac{1}{\sqrt{1(1+n_0)^2}}1\right)m_e$$

where $`m_e`$ is the electron rest mass (0.51 MeV). The fraction of electrons with energy above the Čerenkov threshold is given by :

$$f_{\stackrel{ˇ}{C}}=\left(\frac{0.89E_01.2}{E_0+E_{th}}\right)^s(1+0.0001sE_{th})^2$$

where $`E_0`$ = 26 *MeV* when $`s0.4`$
= $`4417(s1.46)^2MeV`$ for $`s>0.4`$

The number of electrons above the Čerenkov threshold is then computed in each of the slabs. For each straight section of the electron path, the number of Čerenkov photons produced per unit path-length is given by,

$$\frac{dN}{dl}=6.28\alpha \left(\frac{1}{\lambda _1}\frac{1}{\lambda _2}\right)\left(1\frac{1}{\beta ^2n^2}\right)$$

in the wavelength band bounded by $`\lambda _1`$ and $`\lambda _2`$ and $`\beta `$ is the electron velocity and $`\alpha `$ is the fine structure constant. Assuming all the electrons in the cascade to be close to the shower axis and they travel with the velocity close to that of light in vacuum, the number of Čerenkov photons produced in a path-length $`dl`$ (= 333 $`m`$ in the present case) and in the wavelength range 300-550 $`nm`$ at a depth $`xgcm^2`$ is then given by,

$$N_{ph}=0.402N\left(\frac{x}{1030}\right)dl$$

### 2.2 Proton primaries

Because of the differences in the kinematics of shower development in the case of hadronic primaries the longitudinal development profiles too are different as compared to that of $`\gamma `$ray primaries. Using the scaling model for nuclear interactions, Gaisser & Hillas  find that the average number of particles, $`N(E_p,x)`$ at an atmospheric depth of $`xgcm^2`$, measured from the first point of interaction, in a shower initiated by a proton of energy $`E_p`$ (GeV) can be adequately represented by the empirical relationship

$$N(E_p,x)=S_0\frac{E_p}{ϵ}e^{t_m}\left(\frac{t}{t_m}\right)^{t_m}e^t$$

where

$$t_m(E_p)=\frac{x_0}{\lambda }\mathrm{ln}\left(\frac{E_p}{ϵ}\right)1$$

where $`S_0`$=0.045, $`ϵ`$=0.074, $`x_0`$= 37.2 $`gcm^2`$, $`t=\frac{x}{\lambda }`$ and $`\lambda `$ =proton interaction mean free path in air (70 $`gcm^2`$). We chose the electron energy spectrum given by Zatsepin and Chudakov which is independent of the stage of cascade development .

$$F(E,x)dE=0.75N(x)\frac{d\delta }{1+\delta }$$

where $`F(E,x)`$ and $`N(x)`$ are the number of electrons of energy $`E`$ and all energies respectively at a depth $`x`$, $`\delta `$ = 2.3$`E/\kappa `$ and $`\kappa `$ = 72 MeV with $`E`$ in MeV. Even though this approximates well at shower maximum in the energy range 20-300 $`MeV`$, the minor departure at other depths affects only the shape of the lateral distribution of Čerenkov photons and does not affect the present estimates significantly. The integral electron energy spectrum is then derived:

$$F(E,x)=0.75N(x)\mathrm{ln}\left(1+\frac{2.3}{\kappa }E\right)$$

Substituting $`E_{th}`$ for $`E`$, the number of electrons above the Čerenkov threshold at a particular depth are calculated using the above equation.

The longitudinal Čerenkov photon profiles are then derived for $`\gamma `$rays of primary energy 50, 100, 250, 500 and 1000 GeV and protons of energy 100, 200, 500, 1000 and 2000 GeV. The growth curves are simulated using CORSIKA for two primary energies for $`\gamma `$rays (500 GeV & 1000 GeV) and protons (1000 & 2000 GeV). Figure 1 shows a comparison of the simulated longitudinal profiles for $`\gamma `$ rays (a & c) as well as protons (b & d) with calculated profiles. The simulated profiles (histograms) are averaged over 30 showers for $`\gamma `$ rays and 50 showers for protons. The agreement between the simulated and calculated longitudinal shower development profiles in terms of Čerenkov photons is good for $`\gamma `$ray as well as proton primaries of energies considered here. Table 1 summarizes the total number of Čerenkov photons in these showers. The total number of Čerenkov photons agree with that from simulations within 1 & 2% for $`\gamma `$ray and proton primaries respectively.

The depth of shower maximum is expected to increase logarithmically with primary energy both for electromagnetic as well as hadronic showers as more cascade generations are required to degrade the secondary particle energies . Figure 2 shows a plot of the depth of shower maximum as a function of primary energy, both for $`\gamma `$ray and proton primaries in order to illustrate this point as well as to confirm the correctness of the present calculations. The energy dependence of the depth of shower maximum is parameterized as follows:

$$x_{max}=30.3\mathrm{ln}E_\gamma +118.8gcm^2$$

for $`\gamma `$rays and,

$$x_{max}=32.6\mathrm{ln}E_p+118.6gcm^2$$

for protons.

## 3 Photon attenuation in the atmosphere

To study the attenuation of Cerenkov light in the atmosphere Elterman’s atmospheric attenuation model  is used, which provides the attenuation coefficients for the Rayleigh and aerosol scattering as well as ozone absorption in an altitude dependent form for the wavelength range 270-1260 $`nm`$ .

Atmospheric attenuation model computes optical parameters spectrally and with altitude as follows: (1) pure air attenuation parameters are determined by utilizing Rayleigh scattering cross sections with molecular number densities from standard atmosphere; (2) ozone absorption parameters are derived from coefficients applied to a representative atmospheric ozone distribution; (3) seven sets of aerosol measurements are compared and a profile of aerosol attenuation coefficients as a function of altitude is developed. Tabulation permits calculations for vertical path transmission at one kilometer intervals up to an altitude of 50 $`km`$, individually for each attenuating component or for overall atmospheric extinction (molecular + ozone + aerosol). Above 700 $`nm`$ light of night sky (LONS) increases rapidly because of emission lines of OH and water bands in upper atmosphere while the intensity of Čerenkov light drops considerably( inversely proportional to the square of wavelength). At lower wavelengths ($`<300nm`$) light undergoes strong absorption by ozone molecules in the upper atmosphere.

The number of photons, integrated over the bandwidth of 300-550 $`nm`$ (dictated by the photo-tube band-width) and transmitted through a slab (1 $`km`$ thick) using the corresponding extinction coefficient is given by:

$$N=N_1e^B$$

where N is the number of photons transmitted through the slab with absorption coefficient B and $`N_1`$ is the mean number of photons entering the slab. For each altitude, the photon transmission spectrum is convolved with the Cerenkov photon emission spectrum to get the transmitted photon spectrum. Showers pass through several such slabs with different transmission coefficients, resulting in a modified longitudinal development profile for the Cerenkov photons that reach the observation level. The ratio of the total number of photons within the band-width in a shower received at a particular observation level to the total photons produced is defined as the transmission coefficient. The average transmission coefficients ($`T_c`$) are listed in tables 2 and 3 for $`\gamma `$ray and proton primaries respectively for Pachmarhi altitude. $`T_c`$ can be expressed as a power law in primary energy, E as

$$T_c=aE^b$$
(1)
The values of $`a`$ and $`b`$ are listed in table 4 both for $`\gamma `$ray and proton primaries at three different observation levels. The variation in the average transmission coefficients for Čerenkov photons in the atmosphere is shown in figure 3 for (a) $`\gamma `$ray and (b) protons as a function of primary energy at three different observation levels, illustrating the power law behavior.

The increase in $`T_c`$ with primary energy implies that higher energy primaries penetrate deeper in the atmosphere and hence pass through lesser air mass. The proton primaries reach the shower maximum lower down in the atmosphere compared to a $`\gamma `$ray primary of the same energy since the interaction mean free path in air for the former is nearly twice the radiation length. As a result the protons have marginally larger ($`1.2\%`$ at 100 GeV to $`0.6\%`$ at 1 TeV) average transmission coefficient compared to a $`\gamma `$ray primary of the same energy.

The attenuation of optical photons due to Rayleigh and aerosol scattering is more significant at lower altitudes. As a result the Čerenkov photon transmission is expected to depend on the observation level as also seen in figure 3. The average photon transmission coefficient increases almost linearly with decreasing atmospheric depth and the rate is $`0.1\%`$ for every 100 $`gcm^2`$ of air.

Figure 4 shows the longitudinal development profiles of Čerenkov photons at production in the atmosphere and that for the photons detected at the observation level ($`i.e.`$ “without” and “with” atmospheric attenuation correction) for the same wavelength band. The differences in the two profiles is due to the wavelength dependent absorption of Čerenkov photons in the atmosphere. From the plots it is clear that the shape of the longitudinal profile remains largely unchanged (except for the total number of photons). The range of atmospheric heights from which the detected photons originate also remains unchanged as well, irrespective of the primary energy or species. Thus the height of shower maximum remains unchanged due to atmospheric attenuation.

The average path-length of Čerenkov photons reaching different core distances differ resulting in different attenuation. Figure 5 shows such a radial variation of the transmission coefficient for Čerenkov photons generated by $`\gamma `$rays and protons of two different energies. Also shown are the radial variation of mean production heights. It can be seen that the radial variation of transmission coefficients is rather small but varies inversely as the mean production height.

## 4 Čerenkov photon spectrum

The Čerenkov photon spectrum at the observation level is an important input for the design of an atmospheric Čerenkov experiment. The fraction of the photon spectrum bracketed by the photo-multiplier bandwidth has an important bearing on the sensitivity of a TeV $`\gamma `$ray telescope. Hence we computed the photon spectrum at different observation levels for $`\gamma `$ray and proton primaries of various energies. The bandwidth considered in our calculation is 270-550 $`nm`$. Figure 6 shows photon spectra, corresponding to Pachmarhi level, generated by (a) 50 GeV $`\gamma `$rays & 100 GeV protons (\* & diamond respectively) as well as (b) 1 TeV $`\gamma `$rays & 2 TeV protons. It can be seen from the figure that the wavelength at peak intensity is a function of the primary energy in both the cases. It shifts from around 350 $`nm`$ at 50 GeV to around 330 $`nm`$ at 1 TeV. This is because the higher energy primaries reach the shower maximum lower down in the atmosphere compared to lower energy primaries. As a result, absorption of shorter wavelength photons which is primarily due to atmospheric Ozone is comparatively larger for lower energy primaries.

### 4.1 Fraction of UV component at various altitudes

The UV filters were used in atmospheric Čerenkov experiments  whose purpose was two-fold: firstly, it helps reduce night sky background while minimizing the loss of Čerenkov light in the blue and near UV range. This would provide better stability at higher photo-tube gains which in turn will improve the sensitivity of an atmospheric Čerenkov telescope. Secondly, it could serve to identify showers with a larger UV fraction in the Čerenkov light. Since the Čerenkov light generated by proton primaries traverse lesser air mass as compared to $`\gamma `$ray primaries it is expected to have a larger UV content. This property could be exploited to discriminate against hadronic showers .

Hence the UV fraction in Čerenkov light at observation level is a useful signature. So we estimated this fraction for primaries of various energies. We divided the Čerenkov photon spectrum as seen at the observation level into two groups $`viz.`$, the UV range comprising wavelengths 270-300 $`nm`$ and the visible range comprising the wavelength band 300-550 $`nm`$. The ratio of the number of photons in the UV range to that in the visible range is defined as $`R_p`$ for proton primaries and $`R_\gamma `$ for $`\gamma `$ray primaries. Table 5 summarizes the $`R_\gamma `$ and $`R_p`$ values for different primary energies as seen at three different observation altitudes.

Figure 7 shows a plot of $`R_p`$ and $`R_\gamma `$ as a function of primary energy for three different observation levels. Also shown in figure 7 are the relative excess of UV content in proton showers *vis-a-vis* $`\gamma `$ray showers, as a function of primary energy at three different altitudes.

It can be readily seen from the figure that sensitivity to UV excess is better at higher altitudes as compared to sea level.

Figure 8 shows radial variation of the ratio of UV to visible photons for $`\gamma `$ray and proton primaries at an altitude of 1 *km*. The relative UV fraction is marginally higher around the hump region since the photons in this region are contributed by higher energy electrons which have shortest path-length in air. The radial dependence of the UV fraction is more pronounced for proton primaries for reasons already mentioned before.

## 5 Discussion and Conclusions

In spite of the complexity of the interaction kinematics of high energy cosmic rays in the atmosphere, the average behavior as derived from detailed simulation studies agrees reasonably well with analytical calculations. This demonstrates that the simulation package does take into account almost all the interaction characteristics giving credence to the conclusions drawn from simulation studies. The position of shower maximum obtained by us agrees well with those of Miller & Westerhoff .

Armed with the above result we proceed to make analytical calculations for the production of Čerenkov light produced by both $`\gamma `$rays and protons of various energies at three different observation altitudes. This is ideal to study the average shower properties as this is much faster than the detailed simulations. We then apply a detailed wavelength dependent corrections due to the attenuation of Čerenkov photons during their propagation in the atmosphere. This has been done for $`\gamma `$ray and proton primaries of various energies and for three observation levels. It has been found that the fraction of the transmitted Čerenkov photons at an observation level of 2 $`km`$ above sea level is about 36% more than that at sea level. This is primarily because the average distance to shower maximum is reduced resulting in lesser absorption. Even though the effective collection area is reduced at a higher altitude, the light intensities are larger near the shower core. In addition, increased UV content in the Čerenkov light makes higher observation levels better suited for atmospheric Čerenkov work. Simulations carried out for Mt. Hopkins observatory site (altitude of 2.3 $`km`$) suggest an advantage factor in the range 1.2-1.4 compared to sea level , is consistent with the present result. We also derived the observed Čerenkov photon spectrum for a $`\gamma `$ray of energy 3 TeV at an observation altitude of 2 $`km`$. This spectrum agrees well with that reported by Mirzoyan, et al. for HEGRA observation level.

The use of a UV filter to preferentially accept Čerenkov light has been suggested even though its efficiency in improving the signal to noise ratio of an atmospheric Čerenkov telescope is in doubt . On the other hand pushing the sensitivity to near UV wavelengths does improve the system sensitivity . Also using suitable filters to suppress the photo-tube sensitivity in the visual band, existing imaging telescopes have been modified so that observations could be made during moderately moonlit nights. This technique has been shown to be successful in detecting $`TeV\gamma `$ray signal from Crab and Mkn 421 despite higher energy threshold and reduced sensitivity .

It has been observed by Zyskin,et al.,  that the relative signal strengths in the UV range (200-300 $`nm`$) to that in the visible range (300-600 $`nm`$) increases with the angular size of the image in the visible. It increases up to $`(8\pm 1.8)\%`$ for large images($`\sqrt{a.b}0.5^{}`$ where $`a`$ and $`b`$ are the semi-major and semi-minor axis of the image in the optical). This could be due to the fact that those showers whose maximum is closer to the observation level produce larger Čerenkov images. This is equivalent to raising observation level which effectively reduces the attenuation of UV photons as discussed in §4.1 From the present calculations the relative signal strength at a primary proton energy of 4 TeV (energy threshold of Zyskin et al., ) is estimated to be $`10\%`$ which is quite consistent with their measurement. It may be noticed that the wavelength ranges used by Zyskin et al.,  is broader in UV and while the visible band in our estimate was modified to match theirs. However the observed Čerenkov photon spectrum falls steeply below $`300nm`$, as could be seen in figure 6. Extending the wavelength range below 270 $`nm`$ does not change the UV content significantly.

We find that the relative strengths of UV photons to the visible is higher in the case of proton primaries by about 16% at 50 GeV and decreases to 12% at 1 TeV. Hence, the hadron discrimination efficiency based on the UV content in Čerenkov light is relatively better at lower primary energies and at near core distances. Thus, measurement of relative UV content of a shower could be a good parameter in order to discriminate against hadrons especially for large ground based arrays with low energy thresholds ($`2050GeV`$) and employing the wavefront sampling technique . However one has to keep in mind that for this technique to be successful one has to make a very accurate (better than $`1\%`$) estimate of the UV and visible light contents of the shower.