text:
# Orbital dynamics of three-dimensional bars: I. The backbone of 3D bars. A fiducial case

## 1 Introduction

A thorough understanding of the orbital structure in a barred galaxy potential can provide useful insight to the stellar dynamics of barred galaxies and therefore to the dynamical evolution of these objects, as reviewed e.g. by Athanassoula (1984), Contopoulos & Grosbøl (1989), Sellwood & Wilkinson (1993) and Pfenniger (1996). Stable periodic orbits trap around them regular orbits and thus constitute the backbone of galaxy structure (Athanassoula, Bienayme, Martinet et al. 1983). Thus the appearance of a given morphological feature can often be associated with the properties of one of the main families of periodic orbits. In the ’90s, starting with Athanassoula (1992a; 1992b), many papers have pointed out that the gaseous response to steady barred potentials is, to a large degree, determined by the morphology of the periodic orbits in the corresponding stellar case. Thus, orbital and gaseous dynamics are linked. This has provided added incentive for studies of the morphology and the stability of periodic orbits in Hamiltonian systems representing disc galaxies.

Orbital theory has often provided useful information on the structure of galactic bars. Thus it is now understood that a bar is basically due to regular orbits trapped around the so called ‘x1’ periodic orbits, which are elongated along the bar major axis \[Contopoulos & Grosbøl 1989\]. Such orbits do not extend beyond the corotation resonance, and in many cases no suitable elongated orbits can be found beyond the 4:1 resonance. This led orbital theory to predict that bars should end at or before corotation \[Contopoulos 1980\]. Orbital theory was also able to predict - at the right distance from the center - the loops of the near-infrared isophotes (see the case of NGC 4314 in Patsis, Athanassoula & Quillen 1997). Yet not all important morphological features have been so far explained with the help of periodic orbits. Thus orbital theory has difficulties to explain the boxy isophotes surrounding the bars of, mainly, early-type barred galaxies (Athanassoula 1996, Patsis et al. 1997). Another point still under discussion is the morphology of the peanut-shaped bulges observed in edge-on disc galaxies. They are considered by many authors as revealing the presence of a bar, and to be associated with the 2:1 vertical resonance. It is not clear, however, which families can make this vertical structure. Could a bar without a vertical 2:1 resonance be boxy or peanut-shaped when viewed edge-on? Could we have stellar rings out of the equatorial plane at the ILR region? Furthermore, the detailed dynamics of the corotation region and the differences in the vertical structure between fast and slow bars remain open issues.

In this series of papers we use orbital theory to address the above questions. This is a first step towards understanding both the orbital behavior in $`N`$-body models and the responses of gaseous discs to potentials derived from near-infrared observations. The differences between our model and the well studied corresponding 2D case of the Ferrers bar \[Athanassoula 1992a\] reflect the changes due to the inclusion of the third dimension. In separate papers we address the question of the morphology of the peanut shaped-bulges (Patsis, Skokos & Athanassoula, 2002a - paper III) and of the boxy isophotes of bars seen face-on (Patsis, Skokos & Athanassoula, 2002b - paper IV).

Our first goal is to make a thorough study of the orbital structure in 3D barred potentials, to classify the important families, and to follow their morphological evolution as a function of the Jacobi integral. We start with a fiducial case. Many of the families we find in this model have been previously mentioned (e.g Heisler, Merritt & Schwarzschild 1982; Pfenniger 1984, 1985b; Martinet & de Zeeuw 1988; Hasan, Pfenniger & Norman 1993). However, other, equally important families, have not yet been studied.

Studying the orbital stability in a Hamiltonian system approximating the dynamics of a barred galaxy we get the periodic orbits that could be used as building blocks for a density model. The general rule is to look for stable periodic orbits, since they trap around them the regular orbits. Not all of them, however, are equally important. The isodensities of the model we use show us the topological limits within which we should look for the significant orbits. Stable representatives of families of periodic orbits which do not support the imposed morphology, i.e. that of a bar embedded in an axisymmetric disc with a central bulge, should be considered as less important. As an example, in 2D models, let us mention the case of the retrograde family x4, which is stable over a very large interval of energies (Athanassoula et al. 1983). This family, however, should get a minimum weight when one tries to construct a self consistent model using a Schwarzschild (1979) or a Contopoulos & Grosbøl (1988) method. In a model of a 3D disc galaxy, besides the counter-rotating x4 family on the equatorial plane, one has to filter out also stable orbits with large $`\overline{|z|}`$, i.e. orbits with large mean vertical deviations, since these orbits do not contribute much to the density of the barred galaxy.

This paper is organized as follows: In section 2 we review briefly the parts of orbital theory that are necessary for understanding this paper. In particular we explain the use of characteristic and stability diagrams in following the dynamical evolution of a family of periodic orbits. We describe also the various types of instabilities encountered in 3D Hamiltonian systems and we introduce the nomenclature of the main families. The latter is necessary since a number of the families presented here have not been previously discussed and thus need to be incorporated in a unique nomenclature scheme. In section 3 we introduce our 3D model and the orbital structure in a 2D counterpart. In section 4 we present the main families x1, x2 and x3 and their bifurcations. In section 5 we describe the orbits around L4 (or L5) and around L1 (or L2), as well as families outside corotation. We conclude in section 6.

## 2 A short introduction to periodic orbits in the present context

### 2.1 Periodic orbits and their stability

In this section we will briefly review some parts of orbital theory which are necessary for the understanding of this paper. A clear, easily readable introduction to the subject has been given by Sellwood & Wilkinson (1993). We also refer the reader to the pioneering works of Pfenniger (1984, 1985b) and Contopoulos & Magnenat (1985).

We study the stability of simple-periodic orbits in a barred potential in cartesian coordinates. The 3D bar is rotating around its short $`z`$ axis. The $`x`$ axis is the intermediate and the $`y`$ axis the long one. The system is rotating with an angular speed $`\mathrm{\Omega }_b`$ and the Hamiltonian governing the motion of a test-particle can be written in the form:

$$H=\frac{1}{2}(p_x^2+p_y^2+p_z^2)+V(x,y,z)\mathrm{\Omega }_b(xp_yyp_x),$$
(1)
where $`p_x,p_y,`$ and $`p_z`$ are the canonically conjugate momenta. We will hereafter denote the numerical value of the Hamiltonian by $`E_j`$ and refer to it as the Jacobi constant or, more loosely, as the ‘energy’. The corresponding equations of motion are:

$`\dot{x}=p_x+\mathrm{\Omega }_by,`$ $`\dot{y}=p_y\mathrm{\Omega }_bx,`$ $`\dot{z}=p_z`$
$`\dot{p_x}={\displaystyle \frac{V}{x}}+\mathrm{\Omega }_bp_y,`$ $`\dot{p_y}={\displaystyle \frac{V}{y}}\mathrm{\Omega }_bp_x,`$ $`\dot{p_z}={\displaystyle \frac{V}{z}}`$ (2)
The space of section in the case of a 3D system is 4D. The equations of motion are solved for a given value of the Hamiltonian, starting with initial conditions $`(x_0,\dot{x}_0,z_0,\dot{z}_0)`$ in the plane $`y`$=0, for $`\dot{y}>0`$. The next intersection with the $`y`$=0 plane with $`\dot{y}>0`$ is found and the exact initial conditions for the periodic orbit are calculated using a Newton iterative method. A periodic orbit is found when the initial and final coordinates coincide with an accuracy at least 10-10. The integration scheme used was a fourth order Runge-Kutta scheme.

The estimation of the linear stability of a periodic orbit is based on the theory of variational equations. We first consider small deviations from its initial conditions, and then integrate the orbit again to the next upward intersection. In this way a transformation $`T:𝐑^4𝐑^4`$ is established, which relates the initial with the final point. The relation of the final deviations of this neighboring orbit from the periodic one, with the initially introduced deviations can be written in vector form as: $`\stackrel{}{\xi }=M\stackrel{}{\xi _0}`$. Here $`\stackrel{}{\xi }`$ is the final deviation, $`\stackrel{}{\xi _0}`$ is the initial deviation and $`M`$ is a $`4\times 4`$ matrix, called the monodromy matrix. It can be shown that the characteristic equation is written in the form $`\lambda ^4+\alpha \lambda ^3+\beta \lambda ^2+\alpha \lambda +1=0`$. Its solutions $`(\lambda _i,i=1,2,3,4)`$ obey the relations $`\lambda _1\lambda _2=1`$ and $`\lambda _3\lambda _4=1`$ and for each pair we can write:

$$\lambda _i,1/\lambda _i=\frac{1}{2}[b_i\pm (b_i^24)^{\frac{1}{2}}],$$
(3)
where $`b_i=1/2(\alpha \pm \mathrm{\Delta }^{1/2})`$ and $`\mathrm{\Delta }=\alpha ^24(\beta 2)`$.

The quantities $`b_1`$ and $`b_2`$ are called the stability indices. If $`\mathrm{\Delta }>0`$, $`|b_1|<2`$ and $`|b_2|<2`$, the 4 eigenvalues are on the unit circle and the periodic orbit is called ‘stable’. If $`\mathrm{\Delta }>0`$, and $`|b_1|>2`$, $`|b_2|<2`$, or $`|b_2|>2`$, $`|b_1|<2`$, two eigenvalues are on the real axis and two on the unit circle, and the periodic orbit is called ‘simple unstable’. If $`\mathrm{\Delta }>0`$, $`|b_1|>2`$, and $`|b_2|>2`$, all four eigenvalues are on the real axis, and the periodic orbit is called ‘double unstable’. Finally, $`\mathrm{\Delta }<0`$ means that all four eigenvalues are complex numbers but off the unit circle. The orbit is characterized then as ‘complex unstable’ (Contopoulos & Magnenat 1985, Heggie 1985, Pfenniger 1985a,b). We use the symbols $`S,U,D,\mathrm{\Delta }`$ to refer to $`stable`$, $`simpleunstable`$, $`doubleunstable`$ and $`complexunstable`$ periodic orbits respectively. For a general discussion of the kinds of instability encountered in Hamiltonian systems of N degrees of freedom the reader may refer to Skokos (2001).

The method described above has been initially presented by Broucke (1969) and Hadjidemetriou (1975), and has been used in studies of the stability of periodic orbits in systems of three degrees of freedom. The reader is referred to Pfenniger (1984) and Contopoulos & Magnenat (1985) for an extended description.

A diagram that describes the stability of a family of periodic orbits in a given potential when one of the parameters of the system varies (e.g. the numerical value of the Hamiltonian $`E_j`$), while all other parameters remain constant, is called a ‘stability diagram’ (Contopoulos & Barbanis 1985, Contopoulos & Magnenat 1985). With the help of such a diagram one is able to follow the evolution of the stability indices $`b_1`$ and $`b_2`$, and the transitions from stability to instability or from one to another kind of instability. We will loosely refer to the $`b=2`$ and $`b=2`$ lines on a stability diagram as the $`b=2`$ and $`b=2`$ axes. The S$``$U transitions, when one of the stability indices has an intersection with the $`b=2`$ axis, or tangencies of the stability curves with the $`b=2`$ axis, are of special importance for the dynamics of a system. In this case a new stable family is generated by bifurcation of the initial one and has the same multiplicity as the parent family. That means that the periodic orbits of the bifurcating family have, before closing, as many intersections with the plane $`y`$=0, for $`\dot{y}>0`$, as the orbits of the parent family. The new family may play an important role in the dynamics of the system. S $``$U transitions after an intersection of a stability curve with the $`b=2`$ axis, or tangencies of a stability curve with the $`b=2`$ axis, also generate a stable family but are accompanied by period doubling. This means that the number of intersections with the plane $`y`$=0 (always with $`\dot{y}>0`$), needed for the periodic orbits to close, is double the corresponding number of the parent family. Since the most important families we examine here are simple-periodic, i.e. of multiplicity 1, intersections or tangencies of their stability indices with the $`b=2`$ axis introduce in the system families of orbits with multiplicity 2. U$``$D and D$``$$`\mathrm{\Delta }`$ transitions do not bring new stable families in the system and thus in principle are only of theoretical interest. As we will see, however, the evolution of a family which is found to be initially unstable may be very important for the dynamics of our model. The family could simply become stable in another energy interval, or it may play a major role in a collision of bifurcations, an inverse bifurcation or other dynamical phenomena \[Contopoulos 1986\]. Finally in the case S$``$$`\mathrm{\Delta }`$ we have in general no bifurcating families of periodic orbits.

Another very useful diagram is the ‘characteristic’ diagram \[Contopoulos & Mertzanides 1977\]. It gives the $`x`$ coordinate of the initial conditions of the periodic orbits of a family as a function of their Jacobi constant $`E_j`$. In the case of orbits lying on the equatorial plane and starting perpendicular to the $`x`$ axis, we need only one initial condition, $`x`$, in order to specify a periodic orbit on the characteristic diagram. Thus, for such orbits this diagram gives the complete information about the interrelations of the initial conditions in a tree of families of periodic orbits and their bifurcations. However, even for orbits completely on the equatorial plane, but not starting perpendicular to the $`x`$ axis we need to give initial conditions as position–velocity pairs $`(x,\dot{x})`$ and the characteristic diagram is three-dimensional $`(E_j,x,\dot{x})`$. In the general case of orbits in a 3D system, one has a set of four initial conditions and the characteristic diagram is five-dimensional. The representation of such a diagram is difficult, but when necessary we will give just the $`(E_j,x)`$ projection. $`(E_j,x)`$ diagrams that can be compared with the corresponding 2D models will always be given. In all characteristic diagrams the region to which the orbits are confined is bounded by a curve known as the zero velocity curve (ZVC), since the velocity on it becomes zero.

### 2.2 The nomenclature of the main families

Our orbital study is more extended than previous ones and thus brings in new families of orbits which have not been studied so far. We were thus brought to introduce a new nomenclature system, extension of the Contopoulos & Grosbøl (1989) system, which covers all the new types of orbits.

For the main 2D families of simple periodic orbits the nomenclature in the present paper follows the standard notation of Contopoulos & Grosbøl (1989). We thus have the x1 family, where orbits are elongated along the bar and which is the main family in the case of barred potentials, families x2 and x3, whose orbits are elongated perpendicular to the bar, and the retrograde family x4. 2D families bifurcated from x1 at the 3:1 resonance region on the equatorial plane are denoted by t1, t2, …, for consistency with the names used in Patsis et al. (1997). 2D families bifurcated at the 4:1 resonance region on the equatorial plane are called q1, q2, q3,…. Planar orbits related with the 1:1 radial resonance will be called o1, o2…. They are encountered only in some models. The fiducial case presented in the present paper is not one of them.

Further planar families appear beyond the x1 family, at the gaps of the even resonances 4:1, 6:1, 8:1 etc. They are given the names ‘f’, ‘s’, ‘e’… respectively. These families, not directly related to the morphological problems we address in this series, will be discussed elsewhere.

We name the 3D families bifurcated from the basic family x1 at the vertical resonances as x1v$`n`$, where $`n`$ denotes the order of their appearance in our fiducial model A (see below section 4). This is a convenient model to be used for our nomenclature, since there are families of 3D orbits associated with all basic vertical resonances. So x1v1 is the one bifurcated at the first S$``$U transition, which happens at the vertical 2:1 resonance region, x1v2 is the one bifurcated at the U$``$S transition (second stability transition of the model also at the vertical 2:1 resonance region), x1v3 is the one bifurcated at the S$``$U transition at the vertical 3:1 resonance and so on. Further bifurcations of these x1v$`n`$ families are indicated with an ‘.$`n`$’ (for the $`n`$-th bifurcation) attached to the name of the parent family; i.e. the first bifurcation of x1v1 will be x1v1.1, the second x1v1.2, etc. Further bifurcations of these families will be indicated by further ‘.$`n`$’ attached to the name of the parent family. Thus x1v1.1.1 is the first bifurcation of x1v1.1. The naming system is thus extendable at will.

In general at each vertical resonance we have two bifurcating families introduced in the system. The number of oscillations along the rotation axis $`z`$ corresponds to the vertical resonance at which the family is born. E.g. families x1v1 and x1v2, which are bifurcated at the vertical 2:1 resonance region, have orbits with two oscillations along the $`z`$ axis. This determines only partially their morphology, since the bifurcating family can be introduced either in the $`z`$ or the $`\dot{z}`$ coordinate of the initial conditions. If we know the number of oscillations of a family along each axis and also whether it is a bifurcation in $`z`$ or $`\dot{z}`$, then we know its morphology. Families with similar morphology are similar in their corresponding $`(x,y)`$, $`(x,z)`$ and $`(y,z)`$ projections. In the fiducial case, where each vertical resonance is associated with two bifurcating families, the families x1v(2n-3) and x1v(2n-2) are born at the n:1 resonance.

We note, however, that the first vertical bifurcation is not in every model the x1v1 family, as in the fiducial case. In other models (Skokos, Patsis & Athanassoula, 2002 - paper II) it can happen that the first 3D bifurcation of x1 is not related with the 2:1 vertical resonance, but with a different one. In such a model the first vertical bifurcation of x1 will have the same name as the family of the fiducial model which has similar morphology. Equivalently, it will have the same name as the family of the fiducial model which is introduced in the same n:1 resonance and in the same ($`z`$ or $`\dot{z}`$) coordinate. This way we make sure that families with similar morphologies share the same name in the various models. In addition, if for some reason we have more than one vertical bifurcation of x1 associated with a vertical resonance, we introduce appropriate primes in our nomenclature. E.g. in a model with two vertical 4:1 resonances we will have the pairs of bifurcating families x1v5, x1v6 and x1v5, x1v6. By keeping the basic name of the family similar for all families associated with the same vertical resonance, we underline again the dependence of the name on the encountered morphology. Nevertheless, the basic names are given in the fiducial model, which thus becomes a reference case for all our work.

We use the same nomenclature not only for the bifurcations of the basic family x1, but in general for the vertical bifurcations of every 2D family. Their name consists of the name of the parent family, followed by ‘v$`n`$’, where $`n`$ indicates its $`n`$-th vertical bifurcation. Also the names of the bifurcations of the bifurcating families are characterized by the addition of ‘.1’, ‘.2’$`\mathrm{}`$ etc. at the end of the name of the 3D family, as described above for the corresponding families associated with x1.

We will use the same system in order to name also radially bifurcating families. In general a radial bifurcation will be named as ‘wr$`n`$’, where ‘w’ the name of the parent family. E.g. the $`n`$-th radial bifurcation of family f will be ‘fr$`n`$’ (fr1, fr2…etc.).

Let us now turn to orbits related with the axis of rotation. The family on the axis of rotation is called ‘z-axis’ family \[Martinet & de Zeeuw 1988\]. Its two first bifurcations are introduced at the first S$``$U transition and the first U$``$S transition respectively and they are the ‘sao’ and ‘uao’ orbits of Heisler, Merritt & Schwarzschild (1982). This nomenclature, however, does not lend itself to extension which can include what Poincaré (1899) called the ‘deuxième genre’ families (cf. Polymilis, Servizi & Skokos, 1997) which can play an important role in some Ferrers bars (paper II), so we will not adopt it here for other families related with the z-axis orbits. In practice ‘deuxième genre’ orbits are found on the stability diagrams as bifurcations of the parent family when this family is considered as being of higher multiplicity, i.e. if its orbits are repeated many times. Thus, the z-axis family, when its orbits are repeated twice, is called z2. Bifurcations of the z2 family are called z2.1, z2.2 etc. The same rule applies for the bifurcations of z3, i.e. for the bifurcations of z-axis if this is described three times. We then have z3.1, z3.2 and so on. These bifurcating families always come in pairs. A further index (s or u) is attached to their names and is related with their stability.

Around the Lagrangian points $`L_{4,5}`$ we have the long period banana-like orbits, which form a tree of families, and the short period orbits. For the latter we keep the Contopoulos & Grosbøl (1989) notation (spo). For the banana-like orbits we use the notation ban1, ban2,…ban$`n`$ in the 2D cases. Their 2D bifurcations are the families ban$`n`$.1, ban$`n`$.2,… etc. and their 3D bifurcations the families ban$`n`$v1, ban$`n`$v2,… etc. 3D banana-like orbits not related with a 2D one are named banv$`n`$.

A 2D family found around the unstable Lagrangian points $`L_{2,3}`$ is called $`\mathrm{}_1`$.

Throughout the papers we give also the names used by other authors for families that have been previously studied. However, since our study is more extended, there are several families mentioned for the first time.

## 3 The model

### 3.1 The 3D potential

For our calculations we used a 3D potential, which consists of a Miyamoto disc, a Plummer bulge and a Ferrers bar. Pfenniger and collaborators have made extensive use of this potential for orbital calculations (Pfenniger 1984, Pfenniger 1985a;b, Martinet & Pfenniger 1987, Pfenniger 1987, Pfenniger 1990, Hasan, et al. 1993, Olle & Pfenniger 1998). Our work is, in many ways, more extended. We make a much more extensive search for periodic families and we furthermore follow their stability. The latter allows us to find a number of ‘new’ families, which show interesting morphological characteristics. Furthermore, we vary the parameters of the model so that we are able to make comparisons between fast and slow rotating bars as well as between strong and weak bars (paper II). Finally, we focus our work more on tracing the orbital behaviour that could support observed morphological features and less on studying in depth qualitatively the dynamical phenomena that take place in this kind of Hamiltonian systems.

Our general model consists of 3 components. The disc is represented by a Miyamoto disc (Miyamoto & Nagai 1975), the potential of which reads:

$$\mathrm{\Phi }_D=\frac{GM_D}{\sqrt{x^2+y^2+(A+\sqrt{B^2+z^2})^2}},$$
(4)
where $`M_D`$ is the total mass of the disc, $`A`$ and $`B`$ are the horizontal and vertical scale lengths, and $`G`$ is the gravitational constant. The bulge is modeled by a Plummer sphere with potential:

$$\mathrm{\Phi }_S=\frac{GM_S}{\sqrt{x^2+y^2+z^2+ϵ_s^2}},$$
(5)
where $`ϵ_s`$ is the scale length of the bulge and $`M_S`$ is its total mass. The third component of the potential is a triaxial Ferrers bar, whose density $`\rho (x)`$ is:

$$\rho (x)=\{\begin{array}{ccc}\frac{105M_B}{32\pi abc}(1m^2)^2\hfill & \text{f}or& m1\\ & & \\ 0\hfill & \text{f}or& m>1\end{array},$$
(6)
where

$$m^2=\frac{y^2}{a^2}+\frac{x^2}{b^2}+\frac{z^2}{c^2},a>b>c,$$
(7)
$`a`$, $`b`$, $`c`$ are the semi-axes and $`M_B`$ is the mass of the bar component. The corresponding potential $`\mathrm{\Phi }_B`$ and the forces are given in Pfenniger (1984)111We made use of the offer of Olle & Pfenniger (1998) for free access to the electronic version of the potential and forces routines.. They are in a closed form, well suited for numerical treatment. For the Miyamoto disc we use A=3 and B=1, and for the axes of the Ferrers bar we set $`a:b:c=6:1.5:0.6`$, as in Pfenniger (1984). We note that these axial ratios are near the standard values given by Kormendy (1982). The masses of the three components satisfy $`G(M_D+M_S+M_B)=1`$. The length unit is taken as 1 kpc, the time unit as 1 Myr and the mass unit as $`2\times 10^{11}M_{}`$.

In Table 1 we give the parameters of our model. We give it the name A1, and it will be one of the models to be used in our comparative study in paper II.

### 3.2 The 2D Ferrers bar

The general orbital structure in potentials including a 2D Ferrers bar can be found in Athanassoula (1992a). The dynamics are dominated by the presence of the x1 family, which is in general stable. It is characterized by the presence of a narrow instability zone at the 3:1 resonance and a gap at the 4:1 region, which is generally of type 2 \[Contopoulos & Grosbøl 1989\]. The S$``$U$``$S transition at the 3:1 region introduces in the system a couple of simple periodic families of orbits, the importance of which remains local. Beyond the type 2 gap and above the local maximum of the characteristic of x1 at the 4:1 resonance (Fig. 2 in Contopoulos & Grosbøl 1989) one can find a large number of families squeezed close to the zero velocity curve. Finally the families x2 and x3 generally exist for a large energy range and their characteristics form a single bubble. As it is known, x2 is generally stable and x3 unstable.

In the next sections we describe the orbital behaviour in a 3D case where both radial and vertical 2:1 resonances exist. We will thus find the differences introduced in the morphology and stability of the families of periodic orbits by the inclusion of the third dimension. We will also examine how the 3D families of periodic orbits support the bar.

## 4 The x1-family and its bifurcations

### 4.1 A general description

Contrary to the 2D models, where a single family, the x1 family, provides the building blocks for the bar, in 3D models we have a tree of families consisting of 2D and 3D families related to the planar x1 orbits. In Table 2 we summarize the properties of these families. We list their name, the value of the energy at which they are born ($`E_j^{}`$), the $`E_j`$ intervals where they are stable and we indicate whether they are two-dimensional (2D) or three-dimensional (3D). Their interconnections and their role will be described in the following paragraphs.

There are also several 2D families, which are radial bifurcations of x1 and thus part of the x1-tree, but play a less important role in the morphology of the models. They are described in a separate table (Table 3). The ‘t’ families are related with the 3:1 and the ‘q’ with the 4:1 radial resonance region.

Besides the orbits related to the x1 family, we find the x2 and x3 families and their 3D relatives as well. They exist for the same energy intervals as the families of the x1-tree, but their projections on the equatorial plane are elongated along the minor axis of the bar. They are described below.

### 4.2 Families x1, x2 and x3

The characteristics of the x1 and the x2-x3 families in model A1 (Fig. 1) have the typical geometry of the characteristics of 2D Ferrers bars \[Athanassoula 1992a\]. Due to the vertical instabilities, however, x1 becomes unstable over several $`E_j`$ intervals, and not only at the radial 3:1 resonance region, as in the 2D case. In Fig. 1 and in all characteristic diagrams hereafter we draw the unstable regions in light-grey. We observe that the decreasing part of the x1 curve, below the local maximum at the radial 4:1 resonance region ($`E_j0.21`$), is almost everywhere light-grey, indicating that the family is unstable there. The curve at about $`E_j0.17`$ turns back towards lower energies, remaining after that continuously unstable.

The morphological evolution of the x1 orbits is the one expected from the 2D case and is given in Fig. 2. The numbers at the upper right corners of the individual frames correspond to the $`E_j`$ value of the orbit. The orbits are chosen along the characteristic curve starting from the lower values of the Jacobi constant; the orbits in Fig. 2h,i,j belong to the decreasing branch. Except for the instability zones related to the 3:1 resonance all other unstable parts of x1 appear only in the 3D case. As mentioned in section 2.1, the families introduced at the instability strips by bifurcation inherit the kind of stability of the parent family, i.e. of x1. Thus, the instability gaps on the x1 characteristic are covered by the stable orbits of the families born after the corresponding S$``$U transitions. So for almost every energy $`E_j`$ there exists a stable orbit of the x1-tree. As we mentioned in section 2, the 3D bifurcated families are in general characterized by four initial conditions $`(x_0,\dot{x}_0,z_0,\dot{z}_0)`$ so that a $`(E_j,x_0)`$ characteristic diagram cannot provide all the essential information. For this reason we prefer to follow the dynamical evolution of the orbits using stability diagrams.

These diagrams frequently become complicated, but they have the big advantage of giving in a straightforward way the interconnections of the various families, thus becoming a very useful tool in the hunting of periodic orbits.

The evolution of the stability indices $`b_1`$ and $`b_2`$ for x1 are given in Figs. 3, 4, and 5 for successive energy intervals. The arrows denote bifurcated families at the bifurcating points and show the direction of the stability index associated with the S$``$U or U$``$S transition. We observe that the variation of the

index which in Fig. 3 has the larger values for $`E_j<0.38`$ brings in the system the 3D families x1v1, x1v2 etc., while the variation of the other index brings in the families associated with the radial instabilities. The latter remain on the equatorial plane. The variation of their stability indices will in turn bring new families due to vertical and radial instabilities.

The feature depicted in Fig. 4 is typical of the stability diagrams of many of our models. We call this kind of evolution of the stability indices a ‘bow’. The $`b_1`$ and $`b_2`$ curves do not break anywhere, but they evolve in a continuous, rather complicated way, changing direction twice. This ‘bow’ area corresponds to the bend, or elbow, in the characteristic at about $`E_j0.227`$ (Fig. 1), and the complicated evolution of the stability indices happens as we move towards lower $`E_j`$ values along the characteristic curve of x1 at this area. In Fig. 4 one can follow the evolution of $`b_1`$ and $`b_2`$ by following the evolution of both the numbers and the nearby arrows. The lowest value of the stability index at ‘5’, not included in the Figure (indicated only with a dashed arrow outside of figure frame), is $`55`$.

A significant change in the way the 3D bifurcations of x1 are introduced in the system happens at the instability zone found just beyond the local maximum of the $`(E_j,x)`$ characteristic close to the radial 4:1 resonance. As we see in Fig. 3, 4 and 5, the 3D families are bifurcated at S$``$U and U $``$S transitions, where the corresponding stability index intersects the $`b=2`$ axis. Moving on the characteristic towards corotation, before reaching the decreasing branch, a bifurcating family at an S$``$U transition is a stable 3D family with initial conditions $`(x,z,\dot{x},\dot{z})=(a,b,0,0)`$, where $`a,b𝐑`$ and $`a,b0`$. On the other hand, the family bifurcated at the U$``$S transition, is (initially) simple unstable and has initial conditions $`(x,z,\dot{x},\dot{z})=(c,0,0,d)`$, with $`c,d𝐑`$ and $`c,d0`$. This means that the family introduced in the system as stable is a bifurcation at $`z`$, and the simple unstable family a bifurcation at $`\dot{z}`$. For the set of families associated with the vertical 5:1 resonance, on the decreasing branch of the characteristic, this sense of bifurcation is reversed. Namely we have the bifurcation in $`z`$ at the U$``$S transition (x1v8) and the bifurcation in $`\dot{z}`$ at S$``$U (x1v7).

In Fig. 5 we plot the last part of the stability diagram of the x1 family, corresponding to energies higher than $`0.2`$.

As can be seen from the characteristic diagram of Fig. 1, this includes most of the decreasing part of the characteristic, the bend at $`E_j0.173`$ and the part that goes towards lower energies. This part (roughly for $`0.22<E_j<0.173`$), has negative $`x`$ values starting soon after the bend. Heavy arrows and numbers in increasing order on and next to the stability curves in Fig. 5 indicate the evolution of the indices as we move along this part of the characteristic. As we can see most parts are unstable, the short stable parts being drawn with heavy lines. After the turning point, at $`E_j0.173`$, the upper curve, moving now towards lower $`E_j`$ values, is stable until $`E_j0.181`$, then has a part with values smaller than $`2`$ and then reenters the stability region for $`E_j0.197`$. The lower stability curve, however, reaches absolutely large negative values. Thus, the family is always unstable in the parts where $`x<0`$. It is easy to understand how the negative $`x`$ values are introduced by following the evolution of the x1 orbit morphology as we move along the characteristic (Fig. 2). As one moves along the decreasing part of the characteristic (Fig. 2h $``$j), the four apocentra of the orbits develop loops, whose size increases strongly as the energy increases. Already for the orbit in Fig. 2j the loops have become so large, that the sides of the orbit along the bar major axis nearly touch. As we continue along the characteristic they will touch and then cross, so that $`x`$ becomes negative.

Let us now present the evolution of the x2-x3 loop in the 3D case. As seen in Fig. 1, the situation with the x2-x3 characteristic is exactly like in 2D. The stability indices form also loops, as the $`b_1`$ and $`b_2`$ indices of x2 and x3 join each other in pairs (Fig. 6).

In 2D models, the families x1 and x2 are in general the only simple periodic stable families at the x2-x3 area. This is not necessarily the case in 3D models. E.g. in this model, as we can see in Fig. 6b, the 3D family x1v1 has been bifurcated as stable just before the point $`E_j=0.36`$, while close to $`E_j=0.29`$ the family x1v3 is introduced in the system. So the situation at the x2-x3 area is more complicated, since we have there four simple periodic stable families. Since the x2-x3 stability indices form a bubble they have no further intersections with the $`b=2`$ axis and there are no further bifurcations of other simple periodic x2-like families. Both families, however, have tangencies with the $`b=2`$ axis. At these points, as mentioned in the introduction, families of the same kind of stability, but with double multiplicity, will be bifurcated. The one bifurcated from the stable family x2 is interesting. If we put its $`x`$ initial values on the characteristic diagram (Fig. 1), we obtain the extra branch emerging from the x2-x3 loop, pointed with the curved arrow and characterized as ‘x2mul2’. The energy range over which it is stable is indicated with a double arrow above the $`b=2`$ axis in Fig. 6a. Its morphology is given in Fig. 7. The $`(x,y)`$ projection is typical of an x2 orbit, the $`(x,z)`$ one is a fish-like figure

reflecting the double multiplicity of the family, while the $`(y,z)`$ projection offers a shape that could produce a tiny boxy structure in the central region of the bar (note the scale on the axes). The $`(x,z)`$ projection can also offer a boxy structure, if one considers together with every orbit its symmetric with respect to the $`z`$ axis. This, however, is elongated along the minor axis of the bar as will be discussed in paper III. The morphology of this family shows that the model clearly can support in its face-on projection the presence of stellar rings in the x2-x3 area. This, however, is a thick ring structure extending outside the equatorial plane.

### 4.3 The main 3D families

The x1 S$``$U transition at about $`E_j0.36`$ (Fig. 6b), generates the 3D family of periodic orbits x1v1. This family is related with the presence of the vertical 2:1 resonance. It has a stable part close to the bifurcating point, then it has a complex unstable part after an S$``$$`\mathrm{\Delta }`$ transition, and becomes again stable at about $`E_j0.253`$. We have found x1v1 as stable up to $`E_j0.147`$.

The morphological evolution of the family x1v1 is given in Fig. 8. This family corresponds to the z2 family of Hasan et al. (1993) and its orbits have been associated with the appearance of the peanut shaped bulges by Combes, Debbasch, Friedli et al. (1981). Indeed, due to the symmetry of the potential with respect to the equatorial plane, one can find all 3D families in pairs. Thus for x1v1 we will have the smile ($``$) and frown ($``$) types of the $`(y,z)`$ edge-on projections coexisting at a given energy, and the same holds for the $`(x,z)`$ projection. The $`(x,y)`$ projections of the 3D orbits follow in general the morphology of the corresponding x1 orbit of the same energy. As we said in the introduction, the importance of a family of stable periodic orbits is limited as the individual orbits grow in $`|z|`$. The x1v1 orbit for $`E_j=0.2`$ in Fig. 8c exceeds both in its $`(x,z)`$ and its $`(y,z)`$ projections the height of 2 kpc and this means that it cannot contribute significantly to the density of the galactic disc. Its spatial extent on the other hand indicates that this orbit could be used to populate the bulge area.

The U$``$S transition at $`E_j0.343`$ generates the family x1v2 which we followed until $`E_j0.173`$. It remains totally unstable and ends after a U $``$ D $``$ $`\mathrm{\Delta }`$ sequence. It thus doesn’t play any important role in the dynamics of the system.

Family x1v3 (Fig. 9) is stable and its orbits keep low $`|z|`$ values roughly in the interval $`0.293<E_j<0.221`$. It then ends with a S$``$$`\mathrm{\Delta }`$ transition. This family is similar to the z1 family of Hasan et al. (1993). We note that both x1v1 and x1v3 provide useful orbits in the system before their S$``$$`\mathrm{\Delta }`$ transition. This behaviour is also seen in the 3D thick spiral model in Patsis & Grosbøl (1996). Complex instability helps introducing abrupt drops in the density of given features of a model (in our case the peanut), since it stops abruptly the existence of the family responsible for their appearance without bringing new stable families in the system. On the other hand, in cases where a stable family donates its stability to a bifurcation we have a smooth morphological evolution, which can give smooth density profiles in the galaxies. Both x1v1 and x1v3 do not have any intersections or tangencies with the $`2`$ axis and for this reason they do not bifurcate other families with the same multiplicity.

The next bifurcated family is x1v4. This is bifurcated from x1 after a U$``$S transition. We would thus have expected it to be unimportant, since its parent family, x1, is unstable at the bifurcating point. This is the typical behaviour in such cases and we have seen it already happening for x1v2. x1v4 is introduced at about $`E_j0.278`$. One of the two stability indices, let us call it $`b_1`$, remains in the interval $`2<b_1<2`$, while the other, $`b_2`$, goes to negative values smaller than $`2`$. For larger $`E_j`$ values, however, $`b_2`$ increases and for about $`E_j0.224`$, both indices come in the stability zone, i.e. we have $`2<b_{1,2}<2`$. The detailed description of this complicated evolution is beyond the scope of the present paper and does not add anything to the important information that family x1v4 brings stable representatives in the system for $`E_j>0.224`$. The x1v4 family remains stable up to $`E_j0.149`$ where it becomes simple unstable. Its stability indices fold and the family continues existing towards smaller energies. The morphological evolution of x1v4 can be seen in Fig. 10. In Fig. 10a we give the three projections of an unstable orbit close to the bifurcating point from which the family emanates, while in Fig. 10b and

Fig. 10c we give two stable orbits, for energies $`E_j>0.224`$. The last one is for $`E_j=0.206`$ and we see that already the orbit reaches $`|z|`$ values close to 2 kpc away from the equatorial plane. For each orbit of this family there is also a symmetric one with respect to the equatorial plane. If only one of the two is populated, this would give rise to an asymmetric warp-like shape. Populating them both restitutes of course symmetry. The stable orbits of this family enhance the bar, but they deviate substantially from the equatorial plane.

### 4.4 Families at the 3:1 radial resonance

As we already saw, one of the two stability indices bifurcated the 3D families x1v1, x1v2, x1v3 and x1v4, by its intersections with the $`b=2`$ stability axis. The intersections of the second stability index with this stability axis introduce in the system planar 2D orbits. The first family is bifurcated after a S$``$U transition at $`E_j0.244`$, i.e. in the 3:1 resonance region. We call it t1 and it is stable (Fig. 11).

It bridges exactly the instability zone of the x1 in the S$``$U$``$S transition, i.e. its stability indices together with those of the x1, form a bubble \[Contopoulos 1986\]. t1 exists for approximately $`0.244<E_j<0.218`$ and at $`E_j0.218`$ can be considered as an inverse bifurcation222Inverse bifurcation is a non-linear phenomenon encountered in Hamiltonian systems, according to which the bifurcated family, instead of evolving towards the same direction as the parent family, changes direction. It thus extends for the same energies as the parent family before the transition and has the kind of stability of the parent family after the transition. \[Contopoulos 1985\]. of x1. At $`E_j0.214`$, just beyond the ‘bow’ area, the same stability index has another intersection with the $`b=2`$ axis and x1 bifurcates another 2D family, t2. Several 2D and 3D 3:1 type families, related to each other and with x1, are introduced in the interval $`0.214<E_j<0.20`$. Let us briefly mention that, besides t1 (in Fig. 11) and t2, we found a third 2D 3:1 family, t3, which is stable for $`0.2065<E_j<0.2005`$, although it is introduced in the system as simple unstable for $`E_j0.205`$. The morphology of the three 2D families t1, t2 and t3 is given in Fig. 12, and their stable energy intervals in Table 3.

For the lower energies, the t1 orbits have only one loop, which is located on the $`y`$ axis, as the example shown in the left panel of Fig. 12. For higher energies they develop two more loops, symmetric with respect to the $`y`$ axis, and roughly equal in size to the first one. Since the orbits of family t1 are symmetric with respect to the $`y`$ axis, for every orbit we should have also its symmetric with respect to the $`x`$ axis. Combining the two, as in the left panel of Fig. 12, we obtain a shape that is elongated along the bar major axis and resembles the morphology of the x1 orbits with loops, at least for the energies where the orbits have only one loop. The extent of such orbits along the $`y`$ axis reaches up to 4 kpc, i.e. two thirds of the way to corotation.

t2 brings in the system three-dimensional families of periodic orbits with stable representatives. It bifurcates the family t2v1 at $`E_j0.209`$, which in turn bifurcates t2v1.1 at $`E_j0.205`$. The t2v1 family provides stable orbits to the system for $`0.209<E_j<0.207`$ and the t2v1.1 family for $`0.205<E_j<0.203`$. Triangular-like t2-type orbits have characteristic peaks at the sides of the bar, like the peak of the orbits at $`x4`$ in the $`(x,y)`$ projections of Fig. 13. They are near but not always on the minor axis of the bar and their presence can lead to local enhancements of the density at the area between the bar and the L4,5 points.

For any energy in the interval $`0.214<E_j<0.20`$ there are almost always stable 3:1-type orbits of one or the other family. Together with t1, they affect the dynamics of the bar in this region. We note that the 3:1 orbits bifurcated from x1 are very common in all barred potentials and have both in 2D and 3D dynamically only local importance. Orbits of type t2v1 and t2v1.1 have been found even in the early $`N`$-body simulations of 3D bars (Figure 5 in Miller & Smith 1979). The loops of t3 on either side of the major axis of the bar are not of equal size. As can be seen by careful inspection of the t3 orbit in Fig. 12, the loop at the right side of the major axis is slightly bigger than the one to the left. Thus, morphologically, t3 is a kind of asymmetric t1, since for larger energies t1 develops loops which are symmetric with respect to the major axis, besides the one along the major axis.

### 4.5 The last part of the x1-tree

There are two more 3D bifurcations of x1 close to the local maximum of the characteristic at $`E_j0.205`$. It is x1v5 (bifurcated at $`E_j0.213`$ and being stable until $`E_j0.172`$) and x1v7 (bifurcated at $`E_j`$ -0.205, just beyond the peak of the characteristic). The family x1v7 and its bifurcation x1v7.1 provide stable orbits for $`0.205<E_j<0.18`$ and $`0.175<E_j<0.17`$. Nevertheless, the part of this family that contributes to the density of the bar is limited by the fast increase of $`|z|`$ with the energy. Fig. 14 and 15 show the morphology of these families.

In the same region we encounter two more 3D bifurcations of x1, namely the families x1v6, introduced at $`E_j0.211`$ (Fig. 4), and x1v8 introduced at $`E_j0.1925`$ (Fig. 5). Both are born after an U$``$S transition of x1 and remain always unstable. We note that the representatives of x1v5, x1v6, x1v7 and x1v8 families are morphologically similar to those of the B$`z_2`$, B$`\dot{z}_2`$, B$`\dot{z}_3`$ and B$`z_3`$ families of Pfenniger (1984) respectively.

As we have seen, x1 is mostly unstable in the decreasing branch beyond the local maximum at the radial 4:1 gap and the morphology of the orbits at this branch is in general rectangular-like with loops in the corners. There are several families bifurcating from this branch and their orbits have, as already noted for other families, a morphology similar to that of x1 in the region. The 2D families q2 and q3 provide stable asymmetric orbits, two examples of which are given in Fig. 16. We also have one 3D bifurcating family, x1v9, a member of which is shown in Fig. 16. This family also has an asymmetric stable bifurcation for a short energy interval. No stable members of these families can be found outside the interval $`0.186<E_j<0.1808`$. To this we should add the small intervals of stability provided by x1 itself (cf. Fig. 1 and Fig. 5).

Finally, for the sake of completeness we give in Fig. 17 the morphology of the three 3D families, members of the x1-tree, which remain always unstable although they exist for large energy intervals. As we have seen in the corresponding paragraphs they are the families x1v2, x1v6 and x1v8.

## 5 Further families

### 5.1 Orbits around L4 and L5

Another important ‘forest’ of families is the group of the banana like orbits. Here we find the usual planar long and short period banana-like orbits \[Contopoulos & Grosbøl 1989\]. The long period orbits are coming in the system in a large variety of families all of which have stable parts for $`0.1984<E_j<0.1944`$. The stability indices of these orbits exhibit a complicated behaviour having several tangencies and intersections with the $`b=2`$ and $`b=2`$ axes. This brings many families in the system by bifurcation. The family found for lowest $`E_j`$ values is ban1 (Fig. 18a) which is born at $`E_j0.1984`$, followed by ban2 (Fig. 18b,c) that appears at a slightly greater energy value - which in turn bifurcates ban2.1 (Fig. 18d) at $`E_j0.1972`$ \- and ban3 (Fig. 18e,f) introduced in the system at $`E_j0.1983`$. The most important of the planar orbits with stable parts, is ban4 (Fig. 18h), because it is stable over the largest energy interval ($`0.1982<E_j<0.1955`$). It exists for $`E_j>0.1982`$ and it is not bifurcated at this point from any of the families existing for lower energies (ban1, ban2, ban2.1, ban3, ban3.1). From ban4 bifurcates the 2D family ban4.1 (Fig. 18g). The stability indices of ban4 have a complicated behaviour which is typical of a collision of bifurcations \[Contopoulos 1986\]333Collisions of bifurcations happen when both $`b_1`$ and $`b_2`$ are exactly equal to $`2`$ or $`2`$ for a particular set of the control parameters. In order to observe a collision we need to vary continuously a control parameter of our model (i.e. to consider successive individual models), and for all these cases to follow the evolution of the stability indices as a function of $`E_j`$. This practically means that we vary two parameters. If it happens that $`b_1=b_2=2`$ (or 2) for a critical set of the control parameters, then we will observe a change in the interconnections between parent and bifurcating families, before and after the collision. This may also change the general behaviour of the dynamical system.. Approaching $`E_j=0.1955037765`$, the ban4 orbits shrink to L4 (or L5), and beyond this point the short period orbits (spo) grow in size and take their bean-like shape (Fig. 18h and i respectively).

We have also found three 3D families of periodic orbits with stable parts. ban3v1 (Fig. 19a), a bifurcation of ban3 at $`E_j0.1982`$, is initially marginally stable, having one of its two stability indices almost equal to $`2`$, but for $`E_j>0.1962`$ the index become clearly larger than $`2`$. At $`E_j0.1947`$ the two indices join each other and we have a S$``$$`\mathrm{\Delta }`$ transition. ban4v1 (Fig. 19b), a bifurcation of ban4 at $`E_j0.1976`$, is almost everywhere marginally stable in the interval $`0.1976<E_j<0.1944`$. For $`E_j>0.1944`$ it is always complex unstable. ban3v1 and ban4v1 extend to very large $`E_j`$ values, but as complex unstable. Since we have S$``$$`\mathrm{\Delta }`$ transitions there are no bifurcating families and this is the mechanism that terminates the trapping of material around banana-like orbits in our 3D bars. Finally banv1 (Fig. 19c) is introduced in the system at $`E_j0.1957`$ as stable and remains stable up to $`E_j0.1954`$. This family is not obviously related to any other banana-like orbit. Since it is a 3D family we name it banv1.

### 5.2 Orbits around L1 and L2

The L1, L2 Lagrangian points are known to be always unstable \[Binney & Tremaine 1987\]. Around them we find a family of planar periodic orbits we call $`\mathrm{}_1`$. It appears for $`E_j`$ values larger than the one corresponding to L1, the morphology of its orbits resembles that of the spo orbits rotated by $`\pi /2`$, and their periods are of the order of the epicyclic period. Close to the L1 energy and for $`E_j<0.168`$ these orbits are unstable. For $`E_j>0.168`$, however, $`\mathrm{}_1`$ has both stability indices between $`2`$ and 2 and the family becomes stable. Orbits of this family can be found only by starting with initial conditions on the major axis of the bar. For this reason they had not been previously found, since in previous studies searches for periodic orbits started only with initial conditions on the $`y=0`$ axis. In Fig. 20, we plot a few stable orbits of $`\mathrm{}_1`$ and their symmetrics with respect to the $`x`$ axis for $`E_j>0.168`$.

These stable orbits do not support the bar since they are elongated parallel to the minor axis. Nevertheless, they are of physical interest since they support motion parallel to the minor axis, contribute to the exchange of material between regions inside and outside corotation and are able to influence the dynamics in the region between bar and spirals in barred spiral galaxies. The streaming at the apocentra of the $`\mathrm{}_1`$ orbits could support arc-like features beyond the end of the bar.

For larger energies the $`\mathrm{}_1`$ orbits can be observed shifted towards the $`x`$ axis (minor axis of the bar), at about $`E_j0.12`$ they cross the $`x`$ axis and after a short unstable zone they fall on the retrograde family x4 as stable.

### 5.3 Orbits outside corotation

Beyond corotation we find the usual planar families \[Contopoulos & Grosbøl 1989\]. Most of their members display loops. We also find several 3D families with stable parts. As an example we give the family depicted in Fig 21, which is a bifurcation of the planar family called x1(1) by Contopoulos & Grosbøl (1989). The vertical extent of the 3D orbits we found beyond corotation is in general small.

Let us also mention some 2D families, orbits of which are given in Fig. 22. They have been calculated starting with initial conditions on the major axis of the bar as in the case of the $`\mathrm{}_1`$ family and have thus not been described in previous papers. All of them have large stable parts. These orbits have two interesting properties. First, some of them could support motion close to the end of the bar parallel to its minor axis, at radii shorter than the corotation radius. Second, they could efficiently transport material from the outer parts of the disc, e.g. from a distance close to 20 kpc from the center, to the central regions of the bar (e.g. Fig. 22a). This is particularly true for orbits as those shown in Fig. 22a and c.

## 6 Conclusions

In this paper we have made an extensive study of both the 2D and 3D periodic orbits in a fiducial model representative of a barred galaxy. We report on the stability and morphology of the main families. Our main conclusions are:

1. So far the x1 orbits were considered the backbone of bars. This, however, can only be the case for 2D bars, since the x1 can only populate the $`z`$=0 plane. For 3D bars the backbone is the x1, together with the tree of its 3D bifurcating families. Trapping around these families will determine the thickness and the vertical shape of galaxies in and around the bar region. Major building blocks for the 3D bars can be supplied also by families initially introduced as unstable. Thus the family x1v4, introduced in the system after a U$``$S transition, is a basic family, giving stable representatives for large energy intervals in the system.
2. The $`(x,y)`$ projections of the 3D families of the x1-tree retain in general a morphological similarity with their parent family at the same energy. This has important implications for the morphology of a galaxy since it introduces building blocks which have similar morphology as the x1 orbits, but have a considerable vertical extensions. Especially at the regions close to the bifurcating points the $`(x,y)`$ morphology of a x1v$`n`$ family is not only geometrically similar, but actually very close to the morphology of the corresponding x1 orbit.
3. The way the 3D families of the x1-tree are introduced in the system at an instability strip determines the importance of the bifurcations in $`z`$ or $`\dot{z}`$. Particularly in the present model, all 3D families of the x1-tree at the increasing part of the characteristic which are bifurcated in $`z`$ are introduced in the system as stable. On the other hand the stable family associated with the 5:1 vertical resonance (x1v7), bifurcated at the decreasing part of the x1 characteristic, beyond its local maximum, is bifurcated in $`\dot{z}`$. Whether the stable family of the last S$``$U$``$S transition is the bifurcation in $`z`$ or $`\dot{z}`$ determines in a large degree the model’s morphology at its outer parts.
4. The radial 3:1 resonance region provides in the system several 2D and 3D stable families. Their role, however, is locally confined, as in 2D models.
5. 3D orbits elongated along the minor axis of the bar can be given by bifurcations of the planar x2 family.
6. We have found several families of 3D banana-like orbits around L4,5. Their extent is always restricted by a S$``$$`\mathrm{\Delta }`$ transition.
7. Stable families found beyond corotation circulate material between the outer parts of the system and regions as far inwards as 1 kpc. This contributes to the mixing of the elements in a disc galaxy.

The families of periodic orbits we described up to now are indeed the basic families of a 3D Ferrers bar. As we explore the parameter space, however, their properties change, while new important families may appear and play a crucial role. A notable example is z3.1s, a family related to the z-axis orbits along the rotational axis, which will be described in paper II. However, these are rather particular cases and are not encountered in every model.

## Acknowledgments

We acknowledge fruitful discussions and very useful comments by Prof. G. Contopoulos. We thank the referee for useful suggestions that allowed to improve the presentation of our work. This work has been supported by E$`\mathrm{\Pi }`$ET II and K$`\mathrm{\Pi }\mathrm{\Sigma }`$ 1994-1999; and by the Research Committee of the Academy of Athens. Ch. Skokos and P.A. Patsis thank the Laboratoire d’Astrophysique de Marseille for an invitation during which essential parts of this work have been completed.


format: text

text:
# Birkhoff’s Theorem and the Void Matter Density

## 1. Introduction

Voids are structures noted for their pronounced lack of galaxies; concentrations of matter are thought to take the form of a luminous “filamentary” structure, so that together they resemble a “foam” of roughly spherical voids, edged by filaments and walls, and sprinkled with richer concentrations of galaxies where the filaments intersect. While galaxies provide an important tracer of mass in the universe, it is possible that there may be significant amounts of undetected matter in the form of highly ionized clouds, or perhaps a significant background of dark and gaseous matter, either within the filamentary structure, or distributed within the voids. Observations of the low-redshift Ly$`\alpha `$ cloud population have shown that clouds cluster about galaxies (e.g., Tripp et al., 1998; Chen et al., 2001), but there are also many clouds which appear to be isolated.

In Manning (2002) (hereafter Paper 1), a primary cloud catalog (Penton et al., 2000) was separated into complementary void and non-void catalogs through the use of a parameter which assesses cloud environments by weighting the effects of galaxy mass and distance. This parameter, the scalar tidal field, which was summed over galaxies within $`7.5\mathrm{h}_{75}`$ Mpc of a cloud, is optimal for this purpose, as it may easily be calculated from galaxy catalogs. In addition, there is a physical basis for supposing tides affect cloud distributions, since it is the tide which determines the lower limit of cloud density which is stable to dynamical disruption. In Paper 1, a dimensionless form of the tidal field $`𝒯`$ was evaluated at the locations of each of a catalog of low-redshift H I absorbers (Penton et al., 2000), and used to characterize the clouds as members of a void, or non-void environment, according to whether that field was greater or less than some limiting tide $`𝒯_{lim}`$. Equivalent width distribution functions (EWDF) are calculated for each sub-catalog (see Eq. 44, Paper 1). These EWDFs are well-approximated as linear in $`\mathrm{log}d𝒩/dz`$, *vs.* $`\mathrm{log}𝒲`$, where $`𝒩=n(N_{HI})`$, and $`𝒲`$ is the rest equivalent width (EW) in mÅ. One may vary $`𝒯_{lim}`$, calculate the EWDFs of each sub-catalog, and derive the weighted linear fits thereof, based on the number of clouds contributing to the calculation of the EWDF at a given EW. Ditsributions are fitted to the equation, $`\mathrm{log}d𝒩/dz=𝒞+𝒮\mathrm{log}(𝒲/63\mathrm{m}\mathrm{\AA })`$. It was found that the trends of fitting parameters for void and non-void catalogs as a function of $`𝒯_{lim}`$ appear linear at large, or small $`𝒯_{lim}`$ (see Fig. 1), but undergo a non-linear change within a relatively small range of $`𝒯`$. As it happens, there is a corresponding change in the distributions of Doppler parameters at a similar range of $`𝒯`$ for clouds sorted by $`𝒯`$ (see the differential histogram, Fig. 8 in Paper 1). It appears certain that two distinct populations have been isolated, and that a transition zone lies within the range, $`1.3<\mathrm{log}𝒯<0.7`$.

Apart from the interest in the detailed nature of these clouds, and what they say about the nature of voids (paper in preparation), the very fact that they can be sepated says something about the density of voids. For, by Birkhoff’s theorem (Birkhoff, 1923), if voids have a low density, they should behave like a low-density universe, and suffer less deceleration, and have larger local expansion parameters, than the mean. With simulations I calculate the void expansion parameter as a function of the void matter density relative to the average. If the smallness of the range of $`𝒯`$ over which the transition of void to non-void clouds occurs can be used to constrain the outfall velocity, then this can be used to constrain the the void matter density. This is, in short, my plan.

The cosmology assumed for this paper is a standard flat lambda model with $`h=0.75`$. The total matter density of voids is either referred to as $`\mathrm{\Omega }_V=\overline{\rho }_V/\rho _{crit}`$, or as $`f_0=\mathrm{\Omega }_V/\mathrm{\Omega }_m`$, where $`\mathrm{\Omega }_m=0.3`$ is initially assumed. As noted in Paper 1, the large discovered line density of void clouds appears to require that clouds are to a significant degree self-gravitating and discrete. I treat them as such herein.

## 2. Constraints on Cloud Outfall Velocity

Peculiar velocities in galaxy clusters cause a phenomenon known as the “finger of God”. A similar error in attributed position may result when void clouds have large outfall velocities. Large outfall velocities will smear the distinction between void and non-void clouds, as attributed positions are based on the redshift and the mean Hubble constant (see §1). Figure 1 shows the trend of linear fitting parameters for EWDFs defined by various $`𝒯_{lim}`$. It appears to show a fairly small amount of smearing, suggesting that the transition zone (dotted vertical lines) is $`0.13<\mathrm{log}𝒯<0.7`$. How can this observation be used to constrain the outfall velocity? The void filling factor as a function of $`𝒯_{lim}`$ provides the means. Tidal fields were calculated along the pathlength containing the absorption systems (see §5 of Paper 1). Figure 2 shows the fractional pathlength with $`𝒯𝒯_{lim}`$. As a measure of the probability of a random spot in the universe having a tide less than $`𝒯_{lim}`$, this is interpreted as the fractional *volume* with $`𝒯𝒯_{lim}`$. Figures 4 in Paper 1 show that tidal fields increase rapidly at void edges. By raising the tidal field upper limits, we are increasing the volume of what we think of as voids. The bracketing filling factors which correspond to the above range of $`𝒯_{lim}`$ are, $`0.75f_V0.91`$ (see Fig. 2), with a mean $`f_V=0.86`$, corresponding to $`𝒯=0.1`$.

Now, consider an average, fiducial void, which, at $`𝒯_{lim}=0.1`$, has a radius of 15 Mpc ($`11.25\mathrm{h}^1`$ Mpc) (e.g., Lindner et al., 1995). The range of uncertainty in the position of the transition zone corresponds to radii of $`r_V=15(f_V/0.86)^{1/3}=14.33`$, and 15.396 Mpc, for void filling factors $`f_V=0.75`$ and 0.91, respectively; an average variation of $`\pm 530\mathrm{kpc}`$. This uncertainty in the true fiducial radius corresponds to a line of sight (LOS) velocity error of $`H_0\delta r\pm 40`$ $`\mathrm{km}\mathrm{s}^1`$, and would roughly reproduce the same range of uncertainty apparent in the transition zone. Admittedly, many clouds will be moving along vectors at large angles to the LOS, and so even with high outfall velocities, they could have low LOS peculiar velocities. But if clouds are randomly distributed in the void relative to our LOS, and radially outfalling, an average of half of the clouds would have to be moving at an angle less than 60 degrees from the LOS, so that the observed LOS outfall velocity at the edge would be greater than half of the *radial* outfall velocity for half of these clouds. However, if clouds and mass are uniformly distributed in voids, then half of the clouds will be at a radius $`r0.794r_V`$, so that most clouds would have an outfall velocity $`v_{out}>80\%`$ of the maximum. Therefore, we expect that half of the clouds will have a peculiar outfall velocity along the LOS $`>0.4`$ times the radial outfall velocity at the void edge. Since the observed range of the LOS positional error corresponds to a velocity error of $`40`$ $`\mathrm{km}\mathrm{s}^1`$, this would seem to limit the true radial outfall velocity at the void edge to be $`v_{out}<40/0.4=100\mathrm{km}\mathrm{s}^1`$. However, since we know that some of the error in the determination of the transition is plausibly due to peculiar velocities of non-void clouds in proximity to galaxies (see §5.2, Paper 1), this is probably an over-estimation of the upper limit on $`v_{out}`$.

## 3. The Outfall Velocity from Voids

Voids are thought to grow more rapidly than the scale factor (e.g., Bertschinger, 1985; Piran, 1997). This can be conceptually understood by comparing the evolution of the expansion parameter for a flat (the overall universe), and an open universe (the void). For the former,

$$H(z)=H_0\sqrt{\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }},$$
(1)
and the latter, (Scott et al., 2000),

$$H_V(z)=H_V(0)\sqrt{(f_0\mathrm{\Omega }_mz+1\mathrm{\Omega }_\mathrm{\Lambda })(1+z)^2+\mathrm{\Omega }_\mathrm{\Lambda }},$$
(2)
where $`H_V`$ reflects a Hubble constant in the void, and $`f_0`$ is defined in §1. One might think that as $`z`$ approaches infinity, their respective expansion parameters would approach the same value so that the ratio of their current values would simply be,

$$_H\frac{H_V(0)}{H_0}=\left(\frac{\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }}{(f_0\mathrm{\Omega }_mz+1\mathrm{\Omega }_\mathrm{\Lambda })(1+z)^2+\mathrm{\Omega }_\mathrm{\Lambda }}\right)^{1/2},$$
(3)
so that a“peculiar” outfall velocity could be determined by subtracting off the mean Hubble flow,

$$v_{out}=(_H1)H_0r_V,$$
(4)
where $`r_V`$ is the radius of the void. The dashed line in Fig. 3 shows the functional form of Eq. 4 as a function of $`f_0`$.

However, while the theorem says that the lower-density void acts like a separate “universe”, our reconing of it, in terms of the correspondence between lookback time and redshift, is not straightforward, for the lookback time to a given redshift is a function of cosmology, and a flat universe is younger than an open universe with the same $`H_0`$: we cannot so easily derive our answer.

However, we may confidently approach this subject by simulating the void as a “negative” top-hat with perturbation $`\delta `$. A 1-dimensional Lagrangian hydro/gravity code (Thoul & Weinberg, 1995) was acquired courtesy of one of its authors, (A. A. T.), and adjusted for a flat, $`\mathrm{\Lambda }`$ cosmology (see §1). A more detailed description of the code and some of its broader uses will be presented in a subsequent paper (Manning, in preparation). Linearized analysis of voids dictate that the appropriate initial velocity perturbation should be $`v=(1\delta /3)H_V(z)r`$ (Bertschinger, 1985; Peebles, 1993), where $`H_V`$ is the local, void expansion parameter. Simulations of the dark matter component were performed using this velocity perturbation, with a range of uniform underdensities $`\delta `$ embedded in the otherwise flat universe, beginning at $`z=50`$. A $`\delta =0`$ case was run to normalize $`H_0`$. The expansion parameters at $`z=0`$ were derived. The outfall velocity was calculated by subtracting $`H_0`$ from $`H_V`$, and multiplying the result by the 15 Mpc of our fiducial void. Figure 3 shows the results (solid line) , where $`v_{out}`$ is the outfall velocity. For void densities $`f_0<0.2`$, outfall velocities are $`v_{out}>370\mathrm{km}\mathrm{s}^1`$. However, this analysis (§2) finds that $`v_{out}<100\mathrm{km}\mathrm{s}^1`$, and thus, $`f_0>0.75`$. We may apparently conclude that the Hubble flow in voids is very close to the average Hubble flow, and that therefore the difference in density between voids and the mean is minimal.

## 4. Implications for $`\mathrm{\Omega }_m`$

Using the mean filling factors $`f_V=0.86`$ for voids and $`f_F=0.14`$ for filaments, we may analyze the distribution of matter implied by the above. The basic equation for this analysis is,

$$\mathrm{\Omega }_m=f_Vf_0\mathrm{\Omega }_m+f_F\mathrm{\Omega }_F,$$
(5)
where $`f_0\mathrm{\Omega }_m`$ is substituted for $`\mathrm{\Omega }_V`$. Solving for $`\mathrm{\Omega }_F`$, when $`\mathrm{\Omega }_m=0.3`$ and $`f_00.75`$, we find,

$$\mathrm{\Omega }_F<\frac{\mathrm{\Omega }_m(1f_Vf_0)}{f_F}0.761,$$
(6)
which is absurd, since gravitational clustering could not occur unless $`\mathrm{\Omega }_F>1.0`$. It is apparent that $`\mathrm{\Omega }_F`$ must be much greater. Introducing commonly accepted values into Eq. 6, for instance, $`\mathrm{\Omega }_m=0.35`$ and $`f_0=0.18`$, the mean of the range $`0.1<f_0<0.26`$ (e.g., Cen & Ostriker, 1999; El-Ad & Piran, 1997), we find $`\mathrm{\Omega }_F2.2`$. If we assume this more reasonable value of $`\mathrm{\Omega }_F`$, then it turns out that using the above constraint on $`f_0`$, then

$$\mathrm{\Omega }_m=\frac{f_F\mathrm{\Omega }_F}{1f_Vf_0}>0.86.$$
(7)
That is, a low outfall velocity from voids mandates a “high” $`\mathrm{\Omega }_m`$. But is the outfall velocity for large $`\mathrm{\Omega }_m`$ the same function of $`f_0`$? When the 1-D code was run for $`\mathrm{\Omega }_m=0.85`$, a graph consistent with that for $`\mathrm{\Omega }_m=0.3`$ resulted (dashed line, open points, Fig. 3), so the same constraints maintain in a high $`\mathrm{\Omega }_m`$ universe as with $`\mathrm{\Omega }_m=0.3`$.

## 5. Discussion and Summary

### 5.1. Impact on Filaments

First, let us consider the effects of the deposition of energy in volved in the outfall of void clouds onto filaments. If clouds are diffuse, and the kinetic energy of outfall is transformed entirely into random motions, then $`T=\mu m_Hv_{out}^2/3k`$, where $`k`$ is the Boltzmann constant, and $`\mu 0.59`$ is the mass of an average particle. A large velocity $`v_{out}=400\mathrm{km}\mathrm{s}^1`$ would produce a plasma of temperature nearly $`4\times 10^6`$ K, similar to that predicted by Cen & Ostriker (1999) for filaments, suggesting that absorbers should have Doppler parameters $`b>250\mathrm{km}\mathrm{s}^1`$, much larger than observed $`b`$-values of Ly$`\alpha `$ lines of non-void clouds ($`b60\mathrm{km}\mathrm{s}^1`$, Paper 1). For $`v_{out}=100\mathrm{km}\mathrm{s}^1`$, diffuse clouds are converted to a temperature $`T=2.4\times 10^5`$ K; $`b62\mathrm{km}\mathrm{s}^1`$. Interestingly, *Far Ultraviolet Spectroscopic Explorer* observations (Shull et al., 2000) find actual $`b`$values may be roughly half that of the observed Ly$`\alpha `$ line, the excess being attributed to bulk motions. Doppler parameters of order $`30\mathrm{km}\mathrm{s}^1`$ may be explained if clouds are not diffuse, but held by dark matter halos, and hence centrally condensed. In this case, only gas with a density of order that of the filament will be stripped, leaving the central cloud to proceed relatively unmolested (Murakami & Ikeuchi, 1994), minimizing the thermal broadening of the cloud. Shear-induced vorticity (Manning, 1999), will contribute to the Doppler broadening of the cloud, and ram pressure will heat it, so that these absorbers should have greater Doppler parameters than void clouds, as observed (Paper 1, §5.6). According to the scenario in Manning (1999), the induced vorticity of clouds, which are falling into dissipated gaseous envelopes of galaxies and groups of galaxies, introduces an orderly cycle of compression and cooling which results in elevated central H I equivalent widths as the cloud approaches mass concentrations. This inverse correlation between cloud rest equivalent widths and projected galactocentric radii has been observed (Lanzetta et al., 1995; Tripp et al., 1998; Chen et al., 1999, 2001).

### 5.2. Implications for Cosmology

The conclusion that voids have a high matter density depends heavily on the successful separation of void, from non-void clouds. That the Doppler parameter differential histograms, and trends in EWDF fitting parameters both implicate the same transition zone suggests that this is a robust result. Outside the transition zone $`1.3<\mathrm{log}𝒯<0.7`$, the trends seen in Fig. 1 suggest distinct cloud characteristics with few interlopers. Within this range, however, the trends change in a manner consistent with the effects of peculiar velocities of non-void clouds about galaxies, compounded with that of the outfall of clouds from voids with $`v_{out}<100\mathrm{km}\mathrm{s}^1`$ ($`H_V<1.1H_0`$). The simulations have shown that for $`v_{out}<100\mathrm{km}\mathrm{s}^1`$, the ratio of the void density to the mean is $`f_0>0.75`$, requiring $`\mathrm{\Omega }_m`$ to be of order 1. Under the strictures of the flat Lambda model, this requires $`\mathrm{\Omega }_\mathrm{\Lambda }0`$, and $`\mathrm{\Lambda }0`$.

The current standard model postulates two forms of positive energy which, under the analysis of the multipole fluctuations in the CMB, appear to imply that $`\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }1`$ (e.g., de Bernardis et al., 2000). Independent assessments of $`\mathrm{\Omega }_m`$ are all based on concentrations of matter, usually involving dynamical tests, or biasing estimates (e.g. Fukugita et al., 1998; Verde, 2001), and result in estimates $`\mathrm{\Omega }_m0.3`$. However, since these methods are insensitive to a possible background, or additional quantities of matter in void regions, these estimates are necessarily *under*-estimates. As we have seen, constraints on the velocity of outfall from voids imply a rather large “background”.

On the other hand, supernova Ia studies strongly imply an exponential expansion of the universe (e.g., Riess et al., 1998; Filippenko, 2001). How is this to be reconciled? What has not been shown by observation or theoretical derivations is that what is created in conjunction with the exponential expansion, and what enables the universe to be flat, is a “dark” energy. The conclusion that $`\mathrm{\Omega }_m>0.85`$ and $`\mathrm{\Lambda }0`$ can only mean that a dark energy is not formed, but in its stead, ordinary matter, dark and baryonic.

It is unclear how one might adjust the parameters of the void evolution simulations to be consistent with these conditions, for if matter is forming in voids, the tendency of the mass density to rapidly decline is reduced. A careful numerical modeling of void cloud absorbers as remnants of sub-galactic perturbations (paper in preparation) may shed light on this problem (Manning, in preparation).

I thank Ann Thoul for the 1-D code, and Martin White, Mike Cai, and Carl Heiles for useful discussions. I am very grateful to Nahum Arav for helpful advice on the manuscript. I received financial support of NSF grant #AST-0097163 and the UC Department of Astronomy.


text:
# Effects of initial axion production and photon–axion oscillation on type Ia supernova dimming

## I Introduction

Type Ia Supernovae (SNe Ia) are very interesting objects in observational cosmology bruno . Recent observations of distant ($`z1`$) SNe Ia have raised this interest enormously because the data indicate that they are fainter than would be expected for a decelerating universe SCP ; HZT . The dimming seems to be independent of the wave length of the emitted light. The standard interpretation of this result is that our universe is accelerating at present as a result of a non-vanishing cosmological constant. Several alternatives to this idea have been suggested. One idea suggests that the light from distant SNe Ia can be partially absorbed by the dust which may be present on the light path dust . Ideas involving extra dimensions have also been suggested extradim . The effect of cosmic evolution on distant SNe Ia was studied in evolve .

Recently, Csáki, Kaloper and Terning (CKT) ckt1  proposed another explanation of the dimming due to oscillation of photons into very light axions in the background magnetic field of the intergalactic plasma. The statistical significance of this result and possible generalizations of the mechanism are discussed in Ref. erlich .

The intergalactic plasma produces an energy dependent photon–axion oscillation and thus a frequency dependent dimming deffayet ; ckt2 ; MBG . This is in contrast to the achromatic dimming caused by a non-vanishing cosmological constant. While in the future the energy dependence may be the best tool to discriminate between these two mechanicsms, the present data are not sensitive enough ckt2 ; MBG . Based on this fact, it is important to further investigate the CKT mechanism.

In all these works ckt1 ; erlich ; deffayet ; ckt2 ; MBG  it was implicitly assumed that only photons are emitted from the supernova (SN) while axions are produced only via oscillation in the intergalactic medium. If this assumption failed the dimming would be less significant. As an example consider the case of polarized photons, where the initial axion flux is equal to the initial photon flux (with polarization state parallel to the external magnetic field). Then, since the photon to axion and axion to photon oscillation probabilities are equal ($`P_{\gamma a}`$ = $`P_{a\gamma }`$) the oscillation will cause no effect at all. More generally, in the case of unpolarized photons and/or randomly distributed background magnetic field, the photon intensity on Earth is given by

$$I_\gamma (y)=I_\gamma (0)P_{\gamma a}(y)\left[I_\gamma (0)2I_a(0)\right]$$
(1)
where $`y`$ is the distance traveled by the photons, $`I_\gamma (0)`$ $`[I_a(0)]`$ is the initial photon \[axion\] flux, $`I_\gamma (y)`$ is the final photon flux and $`P_{\gamma a}(y)`$ is the oscillation probability. Clearly, if the initial axion flux is significant it must be taken into account.

In this paper we estimate the initial axion flux for the relevant energy range. We consider axion production in the SN, axion conversion in the SN magnetic field and axion conversion in the magnetic field of the host galaxy. We find that it is unlikely that the initial axion flux is significant. We thus conclude that neglecting the initial axion flux does not pose a problem to the axion–photon mixing explanation of the dimming.

## II Photon-axion mixing formalism

In a background magnetic field axions and photons oscillate into each other sikivie ; raffelt  due to the axion–photon coupling

$$_{int}=\frac{a}{M}\stackrel{}{E}\stackrel{}{B},$$
(2)
where $`a`$ is the axion field, $`\stackrel{}{E}`$ is the electric field of the photon, $`\stackrel{}{B}`$ is the background magnetic field and the mass scale $`M`$ characterizes the strength of the axion-photon interactions. In Ref. ckt1  it was found that $`M=4\times 10^{11}`$GeV gives the best fit to the data and we will use this value in the analysis. In order to have mixing a transverse external magnetic field is needed. Moreover, only the plane wave photon polarization state parallel to the external field can mix with the axion.

Photon–axion oscillation is induced by the following mass-squared matrix

$$^2=\left(\begin{array}{cc}\omega _p^2& i\mu \\ i\mu & m^2\end{array}\right).$$
(3)
where $``$ is the photon energy,

$$\mu =\frac{B}{M},$$
(4)
$`B|𝐁|`$ is the projection of the magnetic field on the direction of the photon polarization and $`m`$ is the axion mass. In Ref. ckt1 , the authors found that $`m10^{16}\mathrm{eV}`$ is needed. For the purposes of this paper, as long as $`m^2\mathrm{max}(\mu ,\omega _p^2)`$, we do not need the exact value of the axion mass ($`m10^{16}\mathrm{eV}`$ is small enough). The plasma frequency is given by

$$\omega _p^2=\frac{4\pi \alpha n_e}{m_e},$$
(5)
where $`m_e`$ is the electron mass and $`n_e`$ is the number density of electrons in the plasma. We omitted the Euler-Heisenberg term heisenberg  $`\omega _{EH}^2=(7\alpha /90\pi )(B/B_{cr})^2^2`$ (where $`B_{cr}=4.41\times 10^{13}\mathrm{G}`$ is the critical field strength raffelt ), which contribute to $`_{11}^2`$, since the magnetic fields we are considering are much smaller than $`B_{cr}`$ and the effect of this term is negligible.

For a constant magnetic field the oscillation probability of photons into axions is given by

$$P_1=\frac{\mu ^2}{k^2}\mathrm{sin}^2\left(\frac{kx}{2}\right)=\frac{(\mu x)^2}{4}\left(\frac{\mathrm{sin}(kx/2)}{kx/2}\right)^2,$$
(6)
where $`x`$ is the distance traveled by the photon, $`k=2\pi /L_{Osc}`$ and the oscillation length is given by

$$L_{Osc}=\frac{4\pi }{\sqrt{(\omega _p^2m^2)^2+4\mu ^2^2}}.$$
(7)
Note that the oscillation probability is bounded by

$$P_1\frac{(\mu x)^2}{4}.$$
(8)
Of particular interest to us is the case where the constant magnetic field domain size ($`L_{dom}`$) is small, namely when $`kL_{dom}1`$. Expanding in this small parameter one finds that the bound of Eq. (8) (with $`x=L_{dom}`$) is saturated, namely $`P_1=(\mu L_{dom})^2/4`$. In particular, in this case the conversion is independent of the axion mass, the photon energy and the plasma frequency. This is the reason that the CKT mechanism does not produce energy dependent dimming as long as the plasma frequency is small enough such that the above approximation is valid ckt2 .

Next we consider the case where the magnetic field is not constant. Within a set of simplifying assumptions (see the Appendix for the derivation) the total photon-axion conversion over many domains is given by

$$P_{\gamma a}(y)=\frac{1}{3}\left(1e^{y/L_{decay}}\right),$$
(9)
such that $`y=NL_{dom}`$ is the total length traveled by the photons in $`N`$ domains and

$$L_{decay}=\frac{2L_{dom}}{3P_1},$$
(10)
where $`P_1`$ is given in (6). Using (8) we find the following bound

$$P_{\gamma a}(y)\frac{(y\mu )^2}{8N}.$$
(11)
As discussed in the appendix there are several assumptions that enter the derivation of (9). Still, we do not expect significant deviation from the bound (11) in more realistic cases.

We also consider cases where the mean free path of photons, $`L_{mfp}`$, is not very large. In such cases it takes an average of $`N=3R^2/L_{mfp}^2`$ random steps to escape from a region of radius $`R`$. Assuming that the magnetic field is constant over a region of order $`L_{mfp}`$, the conversion probability is given by Eq. (9) with

$$y=NL_{mfp}=\frac{3R^2}{L_{mfp}},L_{decay}=\frac{2L_{mfp}}{3P_1},$$
(12)
Using (11) we get the bound

$$P_{\gamma a}(y=R)\frac{3(R\mu )^2}{8}.$$
(13)
Note that the bound does not depend on the photon mean free path.

In the following we will make use of the bounds (8), (11) and (13). Of course, they are useful only when the characteristic scale set by the magnetic field ($`1/\mu `$) is much larger than the relevant scales in the problem ($`x`$, $`y`$ and $`R`$). As it turns out this is indeed the situation in the cases we study.

## III Photon–axion mixing inside Supernova Ia

We now turn to study axion–photon mixing inside the supernova and determine under which physical conditions there is an appreciable photon–axion oscillation probability. SNe Ia are identified by the variations of their apparent luminosity with time. The intrinsic peak luminosity of these SNe varies within $`20\%`$, and based on the time variation even this variation can be corrected for. This is the reason that they are used as “standard candles” for various cosmological tests. The spectral analysis of SNe Ia shows the absence of hydrogen. Supernovae of this type are believed to result from the thermonuclear explosion of White Dwarfs (WDs) when their masses are very close to the Chandrasekhar limit sarkar ; pbpal .

In order to estimate the axion–photon oscillation probability inside SNe Ia we consider a simple model pinto . The supernova is assumed to be a sphere of uniform density with an initial radius $`R_010^9\mathrm{cm}`$ that expands with an outer velocity of $`v=c/30=10^9\mathrm{cm}\mathrm{s}^1`$. Photons are emitted uniformly throughout the volume of the supernova and the energy distribution follows Planck’s law pinto ; eastman . The peak luminosity occurs within 10 to 20 days after the explosion. The progenitor stars are magnetized objects. Some WDs are known to host magnetic fields of intensity ranging between $`10^5`$ and $`10^9`$G. The magnetic fields during the supernova event are not very well understood. It is possible that due to the turbulent motions prior to the explosion the WD goes through a stage which might dramatically change the intensity of its initial magnetic field. In particular, if the effect is to increase the magnetic field its value can be as large as $`B_010^{11}`$G bolometric . After the explosion of the WD, the ejecta undergo a large expansion. In a homologous expansion all components of the field decrease like

$$\frac{B}{B_0}=\frac{(R_0)^2}{(R)^2}\frac{(R_0)^2}{(vt)^2},$$
(14)
where the approximation is valid for $`vtR_0`$. This background magnetic field can be chaotic which is a likely result of the turbulent motion prior to the explosion.

We first consider the production of axions in the SN envelope due to photon conversion in its magnetic field. The photon mean free path at peak luminosity is not well known. Ref. eastman , for example, estimates that $`L_{mfp}`$ is in the range $`10^610^{14}\mathrm{cm}`$. The SN radius at that time, which is of order $`10^{15}\mathrm{cm}`$, is larger than the mean free path. The conversion probability at peak luminosity is therefore bounded by Eq. (13). Using (14), we get

$$P_{\gamma a}\frac{3}{8}\left(\frac{B_0R_0^2}{Mvt}\right)^23\times 10^8\left(\frac{B_0}{10^{11}\mathrm{G}}\right)^2\left(\frac{10\mathrm{days}}{t}\right)^2$$
(15)
where for the last step we used $`M=4\times 10^{11}\mathrm{GeV}`$, $`R_0=10^9\mathrm{cm}`$ and $`v=c/30`$. We learn that for values of $`t10`$days and $`B_010^{11}`$G the conversion probability is negligible. Our conclusions would not hold in the unlikely case that initial magnetic fields in SNe Ia reach values significantly higher than $`10^{14}`$G.

Let us now discuss other processes by which axions could in principle be produced inside the SN. It is well known that axion production is potentially large in the first stage of type II supernovae explosions Raffelt:wa . One would then like to check whether a similar situation occurs in the type Ia SNe explosions. However, this is not the case since the temperature and the density of matter in the plasma of SNe Ia at its peak luminosity are relatively very small. Let us elaborate on this.

The potentially relevant processes are the ones that occur at the typical energy scale of a few eV. The axion luminosity due to these processes scales as $`n_B^pn_e^rT^s`$, where $`n_B`$ is the baryon density, $`n_e`$ the electron density, $`T`$ the temperature, and $`p`$, $`r`$, and $`s`$ are process dependent nonnegative numbers. We consider the following processes: axiorecombination — i.e. the reaction, where axions are emitted when electrons and ions form a bound state; axion bremsstrahlung from the electrons and nuclei in the plasma; Compton scattering with axion instead of a photon in the final state; and the axion emission from electrons via the conversion of longitudinal plasmons. The luminosity due to the axiorecombination scales as $`n_Bn_eT^{3/2}`$Dimopoulos:1986kc , the luminosity due to the bremsstrahlung from electrons scales as $`n_Bn_eT^{5/2}`$Raffelt:wa . The axion bremsstrahlung emission from nuclei in the eV energy range is suppressed by additional powers of $`m_e/m_N`$. The luminosity due to the Compton scattering scales as $`n_eT^6`$Raffelt:wa , and the luminosity due to the plasmon conversion scales as $`n_e^2T`$Mikheev:1998bg . The density of plasma in SNe Ia at the time of peak luminosity is of order $`10^{13}\mathrm{g}\mathrm{cm}^3`$ and the temperature is about $`2\times 10^4K`$ eastman . To get a rough estimate of the axion production in the SNe Ia we compare these values with the typical density $`\rho 10^{13}\mathrm{g}\mathrm{cm}^3`$ and the temperature of several MeV inside the core of SNe II Raffelt:wa . Multiplying by the volume of the SNe explosion remnant and taking $`n_{e,B}R^3`$ roughly constant we arrive at a suppression factor of about $`10^{40}10^{50}`$. This expectation is borne out by a more detailed calculation, where it is found that the axion flux due to these processes in the SNe Ia remnant is suppressed by about 50 orders of magnitude compared to the photon flux. We conclude that axion production by these mechanisms can be safely neglected.

## IV Photon–axion mixing in the host Galaxy

Next we consider photon axion conversion in the host galaxy. The details of the background magnetic field of the galaxy are not well understood. We first assume that far away galaxies are similar to the better known nearby galaxies. The size of the galactic disk, which is also the coherence length of the smooth component of the magnetic field is of order $`1\mathrm{kpc}10^{26}\mathrm{eV}^1`$. For a typical galaxy, Ref. carlson  estimates a smooth toroidal component of the magnetic field of magnitude about $`2\mu `$G and a random component with a magnitude of about $`5\mu `$G and a coherence length of order $`50\mathrm{pc}`$. The mean electron density in the host galaxy is taken to be $`n_e=0.03\mathrm{cm}^3`$ carlson ; harari . With these parameters we find that the mean free path of photons is much larger than $`1\mathrm{kpc}`$ and thus they basically do not interact after they leave the SN. For the smooth toroidal component of magnetic field, using $`B=2\mu `$G, we have $`\mu 10^{28}\mathrm{eV}`$. Using (8) with $`x=L_{dom}=1\mathrm{kpc}`$ we conclude

$$P_{\gamma a}\stackrel{<}{}\mathrm{\hspace{0.25em}10}^4.$$
(16)
The random component of the magnetic field is somewhat larger in strength. However, due to smaller coherence length its contribution to the conversion probability is likely to be even smaller. For example, using (11) with $`L_{dom}=50\mathrm{pc}`$, $`y=1\mathrm{kpc}`$ and $`B=5\mu `$G we find

$$P_{\gamma a}\stackrel{<}{}\mathrm{\hspace{0.25em}10}^5.$$
(17)
Far away host galaxies may have considerably different magnetic field than what we find in nearby galaxies. One or two order of magnitude enhancement of the magnetic field could result in a significant conversion. (In order for such conversion to be realized also the plasma frequency has to be small.) Yet, there is no well motivated reason to believe that this is the case. Thus it is likely that the conversion probability in the far away host galaxy is small.

## V Conclusions

We have studied axion production in type Ia supernovae and photon–axion oscillations in their magnetic fields and in the magnetic fields of their host galaxies. If such axion production is significant it can disfavor the photon–axion oscillation explanation of the observed dimming of light from far away type Ia supernovae since then the dimming is less significant, and it may also become energy dependent. We found that all these effects are very small. Our result is robust in the sense that it does not rely on many unknown parameters of the SNe and their host galaxies. For example, we did not have to use exact values of the plasma frequency or the photon mean free path. We had to use only the rough size of the star, its temperature and its magnetic field. We therefore conclude that axion production in the supernovae and their host galaxies has a negligible effect on the CKT explanation of the supernovae dimming.

###### Acknowledgements.

We are grateful to Borut Bajc, Adam Burrows, Arnon Dar, Yossi Nir, Yael Shadmi and Eli Waxman for helpful discussions. J.Z. thanks the Technion theory group for hospitality while this work was in progress. Y.G. was supported in part by the Israel Science Foundation under Grant No. 237/01-1 and by the United States–Israel Binational Science Foundation (BSF) through Grant No. 2000133. J.Z. was supported in part by the Ministry of Education, Science and Sport of the Republic of Slovenia.

## Appendix A Conversion in varying background

In this appendix we discuss photon to axion conversions in the varying magnetic field background. To simplify the discussion it is assumed that photons and axions traverse $`N`$ domains of equal length. In each domain the magnetic field $`B`$ is assumed to be homogeneous with a discrete change in $`B`$ from one domain to the other, i.e. the situation is far from adiabatic conditions. The component of the magnetic field perpendicular to the direction of flight is assumed to have random orientation but is equal in size in each domain. The domains can be (not necessarily) separated by an interval where the conversion is negligible.

Consider an initial state $`c_1(0)|\gamma _1+c_2(0)|\gamma _2+c_a(0)|a`$. The two photon polarization states, $`|\gamma _{1,2}`$, correspond respectively to photons parallel and perpendicular to the magnetic field in the first domain. The initial photon and axion fluxes are then $`I_\gamma (0)|c_1(0)|^2+|c_2(0)|^2`$ and $`I_a(0)|c_a(0)|^2`$ respectively. In the $`n`$-th domain the magnetic field is tilted by an angle $`\theta _n`$ compared to the magnetic field in the first domain

$$|\gamma _{}^n=c_n|\gamma _1+s_n|\gamma _2,|\gamma _{}^n=s_n|\gamma _1+c_n|\gamma _2,$$
(18)
or

$$c_1(y)=c_nc_{}^n(y)s_nc_{}^n(y),c_2(y)=s_nc_{}^n(y)+c_nc_{}^n(y),$$
(19)
where $`c_n\mathrm{cos}\theta _n`$ and $`s_n\mathrm{sin}\theta _n`$. The magnetic field mixes photons polarized parallel to the magnetic field with the axions according to the mass matrix given in (3). We denote the transition elements for one domain $`t_{\gamma \gamma }=\gamma _{}^n(L_{dom})|\gamma _{}^n(0)`$ and similarly for $`t_{\gamma a}`$. The values of the transition elements are equal in each domain since the magnitude of the magnetic field has been assumed to be the same everywhere. The transition probability for photon to axion oscillation in one domain is $`P_1=|t_{\gamma a}|^2`$ with $`P_1`$ given in (6), while the survival probability of photon is $`1P_1=|t_{\gamma \gamma }|^2`$. At the end of the $`n`$-th domain the photon and axion fluxes are

$`I_\gamma (n+1)`$ $``$ $`(1P_1c_n^2)|c_1(y_n)|^2+(1P_1s_n^2)|c_2(y_n)|^2+P_1|c_a(y_n)|^2+\mathrm{}`$ (20)
$`I_a(n+1)`$ $``$ $`P_1c_n^2|c_1(y_n)|^2+P_1s_n^2|c_2(y_n)|^2+(1P_1)|c_a(y_n)|^2+\mathrm{}`$ (21)
where dots represent terms that are proportional to $`\mathrm{cos}\theta _n`$, $`\mathrm{sin}\theta _n`$ or $`\mathrm{cos}\theta _n\mathrm{sin}\theta _n`$. We have defined $`y_n=(n1)L_{dom}`$. Namely, the coefficients $`c_{1,2,a}`$ are taken at the beginning of the $`n`$-th domain.

Next we assume that the transition probability in one domain is small i.e., $`P_11`$, and that the direction of the magnetic field is random, i.e., $`\theta _n`$ is a random variable with $`\theta _{n+1}\theta _n𝒪(1)`$. In this limit $`c_n^2`$ and $`s_n^2`$ can be replaced by their average value $`1/2`$, while $`c_n`$, $`s_n`$ and $`c_ns_n`$ are averaged to zero. Using the fact that $`I_\gamma (n)|c_1(y_n)|^2+|c_2(y_n)|^2`$ and $`I_a|c_a(y_n)|^2`$ we arrive at

$`\left(\begin{array}{c}I_\gamma (n+1)\\ I_a(n+1)\end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}1\frac{1}{2}P_1& P_1\\ \frac{1}{2}P_1& 1P_1\end{array}\right)\left(\begin{array}{c}I_\gamma (n)\\ I_a(n)\end{array}\right)=`$
$`\left(\begin{array}{cc}\frac{2}{3}+\frac{1}{3}\left[1\frac{3}{2}P_1\right]^{n+1}& \frac{2}{3}\frac{2}{3}\left[1\frac{3}{2}P_1\right]^{n+1}\\ \frac{1}{3}\frac{1}{3}\left[1\frac{3}{2}P_1\right]^{n+1}& \frac{1}{3}+\frac{2}{3}\left[1\frac{3}{2}P_1\right]^{n+1}\end{array}\right)\left(\begin{array}{c}I_\gamma (0)\\ I_a(0)\end{array}\right).`$
We can then use the fact that the number of domains is large and replace $`\left(13P_1/2\right)^{n+1}`$ with the limiting expression $`\mathrm{exp}\left[(3P_1y)/(2L_{dom})\right]`$ to arrive at the final expressions

$`I_\gamma (y)`$ $`=`$ $`I_\gamma (0)P_{\gamma a}(I_\gamma (0)2I_a(0)),`$ (24)
$`I_a(y)`$ $`=`$ $`I_a(0)+P_{\gamma a}(I_\gamma (0)2I_a(0)),`$ (25)
with

$$P_{\gamma a}=\frac{1}{3}\left(1e^{y/L_{decay}}\right),$$
(26)
and $`L_{decay}=(2L_{dom})/(3P_1)`$. In cases where the oscillation length is much larger than the domain size the conversion probability $`P_1`$ saturates the upper bound $`P_1=(\mu L_{dom})^2/4`$ and we arrive at the expression $`L_{decay}=8/(3\mu ^2L_{dom})`$ given also in ckt1 .