text:
# Orbital dynamics of three-dimensional bars: IV. Boxy isophotes in face-on views

## 1 Introduction

Boxy isophotes are a typical feature at the end of the bars of early type (SB0, SBa) barred galaxies seen not far from face-on. Typical examples can be found in Athanassoula, Morin, Wozniak et al. (1990) (NGC 936, NGC 4314, NGC 4596), in Buta (1986) (NGC 1433), in Ohta, Hamabe and Wakamatsu (1990) (NGC 2217, NGC 4440, NGC 4643, NGC 4665) and in many other papers. Loosely speaking, the shape of these isophotes is rectangular-like. Their main characteristic is that their small sides, at the largest distance of the isophotes from the galactic center, are roughly parallel to the bar minor axis (Athanassoula 1984, Athanassoula et al. 1990, Elmegreen 1996). In particular, Athanassoula et al. (1990) use generalized ellipses to fit the isophotes of the bars in a sample of early-type strongly barred galaxies, thus describing quantitatively the result that the shapes of the isophotes are rectangular-like rather than elliptical-like. Fig. 1 is a DSS image of NGC 4314 in $`B`$ and demonstrates a typical case of a galaxy with boxy isophotes at the end of its bar.

We observe that the last isophotes of the bar are indeed rectangular-like. Beyond the area of the boxy isophotes starts the spiral structure of this galaxy.

There is a correspondence between the morphology of the boxy isophotes of the early type barred galaxies and the isodensities encountered in snapshots of several $`N`$-body models of bars. This was shown for the first time in a simulation that was run specifically for this purpose (Athanassoula et al. 1990, Fig. 7). In that paper this correspondence was underlined by measuring the rectangularity in the same way as in the observations. Since then there have been several snapshots in $`N`$-body simulations reproducing this morphological feature (see e.g. Shaw, Combes, Axon et al. 1993; Friedli & Benz 1993; Debattista & Sellwood 2000). Recently Athanassoula & Misiriotis (2002) in their $`N`$-body models describe this feature also quantitatively. We can thus conclude that both observations and numerical models clearly show that boxiness close to the end of the bars is a very frequently encountered phenomenon, and thus it should be related to the standard dynamical behaviour in such systems.

Early calculations of orbits in 2D static potentials have underlined the presence of rectangular-like periodic orbits at the 4:1 resonance region (Athanassoula, Bienayme, Martinet et al. 1983; Contopoulos 1988; Athanassoula 1992a). They are either orbits of the x1 family on the decreasing part of the characteristic – towards lower $`x`$ values – in type-2 4:1 resonance gaps111for the nomenclature of the gaps see Contopoulos (1988), or Contopoulos & Grosbøl (1989)., or orbits at the ‘4:1 branch’ beyond the type-1 gap \[Contopoulos 1988\]. These orbits, whenever they exist, support outer boxiness on the face-on views of the models. Their mere presence, however, is not sufficient to explain the observed morphology. Explanations based on the presence of rectangular-like planar 2D x1 orbits suffer from the following problems:

* The range of the Jacobi integral222We will hereafter refer to the Jacobi integral as the ‘energy’. over which one finds stable rectangular-like x1 orbits (in type-2 gaps) or 4:1 orbits (in type-1 gaps) is in general narrow (Contopoulos and Grosbøl 1989). Furthermore, these orbits develop loops at the four corners, whose size increases considerably with energy, while the near-horizontal sections of the orbit approach the minor axis (as in Fig. 3d in Athanassoula 1992a). This happens for energies only a little larger than the energy at which the orbits become rectangular-like. Orbits with loops cannot easily reproduce the observed boxiness. Thus, we have only a very small energy interval with useful orbits (see Fig. 8 in Athanassoula et al. 1990).
* Poincaré sections for an energy value within the small energy interval where rectangular-like orbits are stable, show that the size of the stability area is very small (see e.g. Fig. 21 in Patsis et al. 1997a). This renders the trapping of many non-periodic orbits around stable boxy periodic orbits rather difficult. It is characteristic that Patsis, Efthymiopoulos, Contopoulos et al. (1997b) used dynamical spectra in order to trace tiny stability islands of rectangular-like orbits in a 2D Ferrers bar potential.
* In many models the orbital rectangles are not sufficiently elongated. They frequently are rather square-like, even in cases of strong bars (as in Fig. 3c in Athanassoula 1992a).

Besides the papers of Contopoulos on non-linear phenomena at the 4:1 resonance region in 2D bars in the ’80s, and the work of Athanassoula and collaborators on the morphology of bar orbits in the early ’90s, where the problem is explicitely stated, not much work has been done on this issue. Nevertheless, from figures of orbits in models of 3D bars (Pfenniger 1984; 1985) one can infer that problems like the squareness of the rectangular-like orbits persist even if one considers also the third dimension of the bars. The morphology of single orbits, however, does not determine the boxiness of the isodensities of a model.

The goal of the present paper is to reconsider the problem of rectangular-like isodensities in the framework of 3D orbital structure models, and, more specifically, to examine the contribution of the third dimension to the face-on orbital profiles. Based on the models presented in Skokos, Patsis & Athanassoula (2002a,b; Papers I and II respectively), we investigate the parameters which favour the presence of boxy outer isophotes in the face-on profiles of 3D bars. As in most orbital structure studies, the models we present are not self-consistent. As underlined in Papers I and II, our goal is to study the orbital behaviour, and therefore the morphological changes, as a function of the model parameters. For this we consider even extreme cases, in order to make the effects clearest. We thus determine the main parameters that influence the boxiness of a model. We base our results not on the morphology of single orbits, but on collective appearance when orbits of more than one family are taken into account simultaneously. In parallel to this we study, as counterexamples, cases where the appearance of rectangular-like orbits close to the end of the bar is excluded. Finally, we discuss the correspondence between boxiness observed in edge-on profiles and boxiness observed in the middle of the bars when viewed face-on.

## 2 Face-on profiles

In this series of papers (Papers I, II, and Patsis, Skokos & Athanassoula 2002, hereafter Paper III) we study the basic families in a general model composed of a Miyamoto disc of length scales A=3 and B=1, a Plummer sphere bulge of scale length 0.4 and a Ferrers bar of index 2 and axial ratio $`a:b:c=6:1.5:0.6`$. The masses of the three components satisfy $`G(M_D+M_S+M_B)=1`$ and are given in Table 1. The length unit is 1 kpc, the time unit is 1 Myr and the mass unit is $`2\times 10^{11}M_{}`$. In the present paper we examine two additional models. One of them, model B2, is characterized by a very strong bar, whose mass is 40% of the total, and whose remaining parameters are as in model B (Paper II). The other additional model has different axial ratios than the rest of our model bars. It has $`a:b:c=6:1:0.6`$, instead of $`a:b:c=6:1.5:0.6`$, as all others. It is used to study the contribution of families related to the radial 3:1 resonance to the appearance of boxy isophotes at the end of bars. The rest of its parameters are as in model A1, and we thus call it A1b. The basic properties of the models we present in this paper are summarized in Table 1.

A fundamental conclusion of Papers I and II is that essentially the backbone for building 3D bars is the x1 tree of families of periodic orbits. The tree consists of the x1 planar orbits and of its 3D bifurcations at the vertical resonances. Only in one case (family z3.1s in model B) did we find a family which supports the bar without being introduced in the system after a x1 bifurcation. In the present paper, we use these families in order to present the face-on views of the skeletons of the models. We use for this purpose sets of weighted orbits as in Paper III. As we explained in Paper III, in order to build a profile of weighted orbits for a model, we first calculate a set of periodic orbits and pick points along each orbit at equal time steps. We keep only stable representatives of a family. The ‘mean density’ of each orbit (see §2.2 of paper III) is considered as a first approximation of the orbit’s importance and is used to weight the orbit. We construct an image (normalized by its total intensity) for each calculated and weighted orbit, and then, by combining sets of such orbits, we construct a weighted profile. The selected stable orbits are equally spaced in their mean radius. The step in mean radius is the same for all families in a model. We underline the fact that the orbital profiles we present throughout the paper comprise only stable periodic orbits.

As we have seen, the edge-on orbital profiles (paper III) are of stair-type, which means that the families building the outer parts of the bars have the lowest $`\overline{|z|}`$. We are thus here mainly interested in orbits remaining close to the equatorial plane, since these orbits will contribute more to the surface density at the end of the bars.

### 2.1 Model A1

The fiducial case model A1 offers, also for the face-on structure of the models, a typical example of the contribution of the individual families to the observed face-on orbital morphology. Fig. 2 shows the weighted profiles of all contributing families. It is evident by simple inspection that

the orbits of the 2D family x1 (Fig. 2a) are the most important, mainly because they have stable representatives over a large energy range. Nevertheless, the projections of the 3D families depicted from Fig. 2c to Fig. 2h play an important role. Contrarily, the x2 orbits (Fig. 2b) affect only the central parts of the system. As we can see, boxy features are related to the families x1v1 (Fig. 2c), x1v4 (Fig. 2e), x1v5 (Fig. 2f), x1v7 (Fig. 2g) and x1v9 (Fig. 2h). The family x1v3 (Fig. 2d), on the other hand, has always orbits with elliptical-like projections, thus in a way plays a complementary role, together with x1, providing building blocks for the elliptical-like part of the bar.

The orbits with the rectangular-like projections in model A1 face one of the main problems for explaining the boxy isophotes, i.e. they are less elongated than necessary. In general the boxy isophotes of the early type bars are not squares. The rectangular-like orbits in model A1 have also shorter projections than the x1 ellipses on the semi-major axis. An exception is x1v9 (Fig. 2h). However, the four loops of these orbits are larger than can be admitted by the corresponding shapes of the isophotes of real galaxies Furthermore, the x1v9 family contributes, because of its stability, only over a narrow energy interval ($`0.185<E_j<0.182`$).

It is interesting to note that the boxy x1v1 orbits, responsible in many models for peanut-shaped edge-on profiles, provide to the face-on view of the system a ‘bow tie’ structure. However, in this particular case at least, these orbits remain confined well inside corotation and can be responsible only for an inner boxiness in a galaxy and not for boxy isophotes at the end of the bars.

Fig. 3 combines all x1 related families, i.e. the 2D x1 family and its 3D bifurcations, as well as the x2 and the long-period banana-like orbits.

We can see that the orbits support a bar with a semi-major axis of length about 0.75 of the corotation radius, so that the ratio of the corotation radius $`R_c`$ to the orbital length of the bar $`a_o`$ will be $`R_c/a_o=1.33`$. The longest orbits along the bar are elliptical-like. This however, is not an obstacle to forming more rectangular-like bars, since the elliptical-like orbits extend only little beyond the rectangular-like ones. This extra extent could be suppressed if the outermost periodic elliptical-like orbits were not populated, or it could be masked by orbits trapped around the periodic rectangular-like ones. In the latter case, the isophotes (or rather isodensities) would be more elongated than the corresponding rectangular-like orbits.

We note that in a galaxy or in an $`N`$-body simulation, not all families included in the figure should be necessarily populated. In our non-self-consistent models we try to identify structures in order to seek their corresponding features in galactic images and snapshots of self-gravitating models.

### 2.2 Model A2

The slow rotating bar in model A2 brings new morphological features. In this case square- or rectangular-like orbits play a minor role, as we realize from Fig. 4. Only x1 orbits (Fig. 4b) at the decreasing branch of the characteristic (paper II) and part of the x1v1 family, after the

S$``$$`\mathrm{\Delta }`$ transition (Fig. 4c), provide such stable orbits to the system. For both these families, however, the size of the orbits and the energy range over which they exist do not allow them to play a major role in the orbital structure of the model. In Fig. 4b we choose the contrast of the image such as to allow us to see the loops of the x1 square-like orbits. For even larger energies the loops of the x1 square orbits become huge, and finally the orbits become retrograde. Such orbits are not depicted in Fig. 4b.

Model A2 has two main features: First, the loops along the bar major axis, which are brought in the system by the x1 family and its x1v3 and x1v4 bifurcations. Second, the almost circular (and/or square-like) projections on the equatorial plane of the families x1, x1v4, and x1v5 (Fig. 4b, f, h respectively). The face-on morphology supported by this slow rotating bar is depicted in Fig. 5.

In this figure we combine all orbits of the families presented in Fig. 4. The result is a bar with loops along the major axis surrounded by almost circular orbits. We note that the bar-supporting orbits extend to a distance about 9 from the center (corotation in this case is at 13.24), i.e. in this case $`R_c/a_o=1.5`$.

### 2.3 Model A3

If the bar rotates fast (model A3), the face-on orbital structure changes significantly. The x1 family dominates once again, as we can see in Fig. 6a. Its orbits remain always elliptical-like, but now they support a bar with length 0.95 of the corotation radius, i.e. $`R_c/a_o=1.05`$. The dynamics at the radial 4:1 resonance region are crucial for the morphology of the model (paper II). In this case the families q0 and x1v8 provide the system with rectangular-like orbits quite elongated along the major axis of the bar. Family q0 can be found in two branches, symmetric with respect to the bar’s major axis. We can, however, consider orbits of only one of its two branches if we want in our weighted profiles a non-rectangular parallelogram-like morphology. In order to obtain a desired morphology in our response models, we chose any combination of stable orbits from all available families. In the present case both families q0 and x1v8 have orbits at least as long as the x1 orbits. In addition, families x1v1 and x1v5 also have orbits with boxy projections on the equatorial plane, but reaching distances from the center 1.2 and 3.2 respectively (Fig. 6c,e). That means that, if all boxy families are populated, we can have boxy isophotes in several scales in this model. Even the hexagonal-like projections of the orbits of family x1v3 (Fig. 6d) contribute to the boxiness of the model with their sides which are parallel to the bar major axis.

The face-on view of this family has also the ‘bow tie’ appearance.

The total effect when considering all orbits is given in Fig. 7,

where we take into account all families of Fig. 6. In Fig. 7a we give all weighted orbits together, while in Fig. 7b we apply a gaussian filter in order to show clearly, in a first approximation, the shapes of features, which could be supported by the orbits in the density maps of the models. Fig. 7 tells us that the fast rotating bar model, has a boxy bar with $`R_c/a_o1.05`$.

The fast rotating bar case offers the opportunity to study models with a substantial non-axisymmetric force near corotation. In this way we compensate for the standard shortcoming of Ferrer’s 333 We use Ferrers bars because we are not aware of any other models which are realistic and analytic, and do not have this shortcoming. ellipsoids, namely that their force drops too fast as the radius increases and approaches corotation. Thus A3 allows us to cover also cases with substantial non-axisymmetric forcing near corotation.

### 2.4 Model B

Model B is a model without radial or vertical 2:1 resonances. As we have seen in paper II, the first vertical bifurcation of x1 is x1v5. Thus, the families that build the bar are x1, x1v5/x1v5, x1v7 and the family z3.1s. We remind that this latter family plays a significant role in the orbital behaviour of this particular model (see Paper II) and it was found as a bifurcation of the z-axis orbits when the latter are considered as being of multiplicity 3, i.e. the orbits are repeated three times. The z3.1s family is not related to the x1-tree. The weighted face-on profiles of the above mentioned families are given in Fig. 8. We see that all of them have boxy representatives in their projections on the equatorial plane.

In Fig. 9 we have the orbital face-on view for model B obtained by overplotting the weighted orbits of all families together. In this model $`R_c/a_o1.4`$.

### 2.5 Model D

The strong bar model D is, after the fast rotating bar model A3, the second case where we have periodic orbits that can contribute to sufficiently elongated rectangular-like boxy isophotes. They are, however, considerably less elongated than the rectangular orbits in A3. As we see in Fig. 10, the main contributors are now families x1, x1v5/x1v5 and x1v7.

Again x1v1 orbits have ‘bow tie’ face-on projections, and support a feature of corresponding morphology (Fig. 10b). This feature occupies the main part of the bar, but it does not reach its end, which is at about a distance 0.9 of the corotation radius.

The composite face-on orbital profile of model D is given in Fig. 11.

We see that –although the boxy orbits are present and, if populated, could characterize the morphology of the model– the dominant feature is the loops of the x1 orbits along the major axis. This morphology is enhanced by the loops of the projections of the 3D bifurcations of x1, x1v5/x1v5 and x1v6 (Fig. 10). It should be noted, however, that only the loops of x1v6 and some of the x1 loops extend beyond the rectangular outline. If we omit orbits with such loops we get, for this model, $`R_c/a_o1.23`$, while if we include them we get $`R_c/a_o1.08`$.

### 2.6 Very massive bars

In all the models we studied so far the rectangular-like orbits, if they exist, are rather square-like, except for the fast rotating bar case (model A3), where they are more elongated. The second parameter, after the pattern speed, which made the rectangular-like orbits more elongated is the increase of the relative mass of the bar. For this reason, we studied several more models with increased bar mass fraction and calculated the orbital stability of each individual case. We stopped when we reached models with very large intervals of instability of the x1 family and its 3D bifurcations. From a sequence of models starting with the values of the parameters of model B and with an increasingly large fraction of mass in the bar, we found that this fraction is about 50% of the total mass. It was evident that models with bar masses about 40% of the total mass offered families of periodic orbits with enough stable representatives to support sufficiently elongated rectangular-like isophotes at the end of the bars. E.g. model B2, which differs from model B only in that the mass of the bar GMB=0.4 instead of 0.1 as in model B, has rectangular-like orbits with the ratio of their projections on the major axis of the bar to their projections on the minor axis ($`p_{max}/p_{min}`$) close to 2. This can be seen in Fig. 12.

In Fig. 12a we have orbits from the families x1, x1v5/x1v5 and x1v7, omitting rectangular-like orbits with loops. The blurred image in Fig. 12b includes also orbits of the t1 family (see Paper I) and stable orbits of family x1v6. The latter is initially bifurcated as unstable. For larger energies, however, it has stable orbits with face-on projections with loops along the major axis of the bar. The t1 orbits, bifurcated at the radial 3:1 resonance, introduce in the system some interesting features if one considers both their branches which are symmetric with respect to the minor axis of the bar at a given energy. The t1 orbits are responsible for the guitar-like feature we observe in Fig. 12b closer to the major axis than the outer rectangular-like orbits. For model B2 we have $`R_c/a_o1.04`$.

### 2.7 Support by 3:1 families

The contribution of t1 to the features we found in model B2 gave us the incentive to investigate the contribution of families bifurcated at the radial 3:1 resonance in all our models. We found that in general orbits from the t1 family may support motion roughly parallel to the minor axis of the bar, but not at the end of the bar. Fig. 12b describes the kind of contribution t1 orbits can offer to the overall boxiness of face-on profiles of barred galaxies. We give also an example of the combination of t2 orbits (see Paper I) in a model in which this contribution is pronounced. It is model A1b, a model that differs from the fiducial case (model A1, Paper I) only in the axial ratios, which in model A1b are $`a:b:c=6:1:0.6`$ instead of $`a:b:c=6:1.5:0.6`$. In Fig. 13 we observe

stable x1 orbits with loops along the major axis of the bar, and x1 orbits with a more or less rectangular like shape. The two groups are separated by an empty region caused by an instability region of the x1 family. On top of the x1 orbits we overplot stable orbits of both branches, symmetric with respect to the major axis of the bar, of the t2 family. The t2 orbits contribute to the boxiness of the model by enhancing the borders of the empty region and by forming a rectangular-like feature indicated in Fig. 13 by four black arrows. Model A1b has $`R_c/a_o1.35`$.

## 3 Discussion

We examined the face-on views of all models of Papers I and II and we found families that could support outer, as well as inner, boxy isophotes. In the present paper we do not refer explicitly to model C, since its orbital behaviour does not differ essentially from that of model A1. Furthermore, we use composite profiles of two additional models to demonstrate the effect of very massive barred components and the possible role of families related to the radial 3:1 resonance.

The first quantity that can be used to compare the bars of our models with the bars of real galaxies concerns the length of the bar built by the orbits with respect to corotation radius. Indeed observations (for a compilation see e.g. Athanassoula 1992b, Elmegreen 1996 and Gerssen 2002) as well as hydrodynamical simulations (Athanassoula 1992b), have shown that the ratio $`R_c/a_o`$ takes values only in a restricted range, namely $`R_c/a_o=1.2\pm 0.2`$. This is in agreement with 2D orbital models by Contopoulos (Contopoulos 1980). All our response models have $`R_c/a_o`$ ratios within this range, except for model A2 which has a slow rotating bar and $`R_c/a_o=1.5`$. Table 2 summarizes the $`R_c/a_o`$

ratios for the models we studied in the present paper.

A second quantity we can compare with galaxies and snapshots of $`N`$-body models is the ratio $`p_{max}/p_{min}`$ for the rectangular-like orbits. Eye estimates show that the boxy outer isophotes at the end of the bar of strongly barred galaxies have a ratio $`p_{max}/p_{min}`$ typically larger than 2 (e.g. NGC 936, NGC 4314, NGC 4596). Exceptional cases can be found in the literature but are not many (see e.g. the not-rectified images of NGC 1415 in Garcia-Barreto & Moreno (2000)). A more quantitative study has been made by Athanassoula et al (1990). They find that the axial ratios of the isophotes near the end of the bar are, for all the galaxies in their early-type strongly-barred galaxy sample, considerably larger than 2. On the other hand, in most of our models the boxy orbits are quite square-like. Ratios of boxy isophotes larger than 2 have been mainly found in the fast rotating bar case, where $`p_{max}/p_{min}2.5`$. The strong bar model D has $`p_{max}/p_{min}1.4`$, while in model B2, where the mass of the bar is 40% of the total mass, we attained a $`p_{max}/p_{min}`$ ratio close to 2. In the rest of the models the ratio was less than 1.4.

Let us underline that, for studying the elongation of the rectangular-like isophotes observed in barred galaxies, one needs to combine orbits of several families instead of evaluating the properties of a single orbit or a single family. This is a result of the three dimensional character of our models and of the fact that we have bars built by orbits belonging mainly to the x1-tree. This effect is indicated by isodensities we plotted on some blurred images of our models. Such a typical case is given in Fig. 14, for model B.

The small sides of the rectangular-like structure outlined by the isophotes reflect mainly the contribution of the rectangular-like orbits. However, the two large sides, parallel to the major axis of the bar, are due to the overlapping of orbits of the x1, and z3.1s families.

Our models show that it is the pattern speed that mainly determines the elongation of the outer boxy orbits,. The orbits building the bar of model A3 are rectangular-like and could make a bar with the geometry encountered in the early-type bars with outer boxy isophotes. It is just the increase of $`\mathrm{\Omega }_b`$ that brought in the system the 2D family q0 and the 3D, stable x1v8 orbits. We note the morphological similarity of the q0 orbits with the outer boxy isophotes of NGC 4314 (Quillen, Frogel & Gonzalez 1994), which have the form of a non-rectangular parallelogram (cf. Fig. 7b with Fig. 1). We note that in model A3 the rectangular-like bar ends closer to corotation than the bars of the models for which the longest orbits are the x1 with loops along the major axis. This is consistent with the result found in Patsis et al. (1997a) for the pattern speed of NGC 4314, in which the boxy isophotes are very close to the corotation radius indeed. It is also obvious that, apart from the increase of the pattern speed, the increase of the strength of the bar favoured the elongation of the rectangular-like orbits.

Model A2, which is the slowest rotating case, is the other extreme of the face-on morphologies we encounter in our models. In this model, planar boxy orbits and projections of orbits on the equatorial plane are square-like, and very little elongated along the major axis of the bar (Fig. 4c of the present paper, and Fig.5 in paper II). The development of loops at the apocentra of the x1 elliptical-like orbits, as well as in the projections on the equatorial plane of the x1v3 and x1v4 families, together with the almost circular x1 orbits and the projections of the x1v4 and x1v5 families, give to model A2 a characteristic face-on morphology with the circular-like orbits surrounding all other orbits in the face-on projection. Such a morphology could be linked to the existence of inner rings in barred galaxies.

Inner boxiness of the face-on profiles, much closer to the center than the corotation region, is associated mainly with the x1v1 orbits, i.e. with the 3D family born at the vertical 2:1 resonance. As we have seen in paper III, in several models this family is responsible for peanut-shaped orbital structures in the edge-on views of the models. The face-on orbital skeletons of the models we present here show in another clear way what we noticed in paper III, namely that the x1v1 family builds boxy- or peanut-shaped features which do not approach corotation. Inner boxiness is not rare in galactic bars. Typical examples of boxy isophotes in the middle of the stellar bars are NGC 3992 and NGC 7479 (Wilke, Möllenhoff & Matthias 2000). The x1v1 orbits, or at least their representatives with the largest energy values (see Fig. 2c, Fig. 6c and Fig. 10b), are boxy in their face-on views, but they also have a characteristic ‘bow-tie’ morphology. Inner boxiness is also supported by the z3.1s orbits, at large values of the energy, in model B.

A ‘bow-tie’ morphology in models of the kind we study here can also be introduced by rectangular-like orbits at high energies. These orbits (or more precisely their face-on projections) develop loops, which, in some cases, stay close to the orbits that build the bulk of the bar. This is e.g. the case of the x1v5 family in model B (Fig. 8b). When the rectangular-like orbits develop loops they also become of ‘bow-tie’ shape. Since their deviations from the equatorial plane increase with energy, their projections occupy almost the same area as the rectangular-like orbits of the same family, for lower energies, do. In model B, for $`E_J0.185`$, integrating even chaotic orbits for a small number of dynamical times (of the order of 30) we get a bow tie morphology. This is shown in Fig. 15.

Non-periodic elliptical-like orbits trapped around the x1 family not far from the center of the bar, as well as the projection of orbits trapped around the three-dimensional x1 bifurcations would fill more densely the area between the regions A and B (and symmetrically between C and D), than between the regions A and C (and B and D). This results from the elliptical-like shape of the projections of these orbits on the equatorial plane and the fact that they have their apocentra between the regions A and B (and C and D). We note also the contribution of the x1v1 orbits to a ‘bow-tie’ morphology of face-on profiles, as mentioned previously.

Due to the presence of these families a bow-tie morphology, like the one found by Barnes & Tohline (2001), is not excluded from the usual barred morphology. Nevertheless, we did not find orbits of a single family that can both enhance a rectangular-like structure at the end of the bar without loops at its four corners, and simultaneously enhance a bow-tie morphology. Single families with bow-tie profiles in their face-on projections do not extend to the end of the bar (x1v1 in Figs. 2c, 6c, 10b and x1v3 in Fig. 6d).

## 4 Conclusions

In this paper we discussed boxiness in the face-on views of 3D models, and, in general, their face-on orbital structure. We examined a large variety of possible 3D orbital behaviour, that could contribute to the boxiness in face-on views of barred galaxies. Our main conclusions are:

1. Boxiness in the face-on views of 3D barred models is an effect caused by the co-existence of several families, each contributing appropriate stable orbits. The morphology of boxy isodensities/isophotes is not necessarily similar to the morphology of individual stable, rectangular-like orbits. In some cases (Fig. 14) the iso-contours can be narrower than the orbits.
2. In models with the right morphological parameters, we find appropriate building blocks to account for the rectangular-like isophotes or isodensities seen in early type barred galaxies and in some $`N`$-body simulations.
3. In 3D models the family of the planar x1 orbits is subject to vertical instabilities, and thus in several cases it has considerable instability strips at the 4:1 resonance region. This should be an obstacle for x1 planar orbits to account for the rectangular-like shape of bars, since unstable periodic orbits can not trap regular orbits around them. However, at the instability regions of the x1 we find other stable families, whose ($`x`$,$`y`$)-projected orbital shapes are, at least near their bifurcations, very similar to those of the x1. They have orbits which are stable over large energy intervals and also have ($`x`$,$`y`$)-projected shapes that can enhance a rectangular-like bar outline. Thus the inclusion of the third dimension in the models enhances the possibility of rectangular-like isodensities.
4. There are families of periodic orbits that support boxiness in the outer bar regions, as well as families that support boxiness in somewhat more inner parts. The standard families belonging to the former category are the stable representatives of the x1 orbits (close to the radial 4:1 resonance), and of the families x1v5, x1v7 and x1v9. The families q0 and x1v8 play a major role in model A3 and thus could be essential for fast rotating bars. Nevertheless, the geometry of a boxy feature becomes evident mainly in weighted profiles, where orbits from one or more families are considered. Inner boxiness is associated mainly with the x1v1 family and, in model B, with family z3.1s.
5. Orbits of the families related to the 3:1 resonance (t1 and t2) may contribute in some cases to the boxiness of the profiles. The t1 family can do this by supporting motion parallel to the minor axis of the bar in intermediate distances between the galaxy center and corotation, and the t2 family by enhancing the sides of the boxes parallel to the major axis.
6. The consideration of several families of orbits for building a profile may lead to boxy features close to the end of bars, with $`p_{max}/p_{min}`$ ratios different from the corresponding ratios of individual orbits or families of orbits.
7. An essential conclusion of our investigation is that outer boxiness is favoured by fast bar pattern speeds, while in the slow-rotating model the bar is surrounded by almost circular orbits. These are the two extremes of an orbital behaviour that changes as the pattern speed varies from one model to the other. The near-circular orbits could be building blocks for inner rings.
8. The fast rotating bar has a length 0.95 of its corotation radius, while the slow one only 0.68. This indicates that boxy bars end close to their corotation, while the end in slow rotating bars might be associated with $`n:1`$ resonances of lower $`n`$ values.
9. The rectangular-like orbits in models with faster bars are more elongated than the corresponding orbits in models with slower bars. When modeling individual barred galaxies, for which we know a priori that in general $`1<R_c/a_o<1.4`$, the elongation of the orbits can be used to give estimates of, or at least set limits to, the bar pattern speed, without the use of kinematics, which is not always available.
10. The second most efficient way of stretching rectangular-like orbits is to increase the mass of the bar. This mechanism, however, is restricted by the fact that in very strong bars the families of periodic orbits which support the boxy face-on profiles are unstable over large energy intervals.
11. Outer boxy, ‘bow tie’ morphology is possible in some models by combining orbits of several families. A necessary condition for this is to have 4:1 type orbits with loops close to the main bar. Inner ‘bow tie’ morphology can be due to the presence of x1v1 orbits (Fig. 2c, Fig. 6c, Fig. 10b).

## Acknowledgments

We acknowledge fruitful discussions and very useful comments by G. Contopoulos and A. Bosma. We thank the anonymous referee for valuable remarks, which improved the paper. 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. ChS and PAP thank the Laboratoire d’Astrophysique de Marseille, for an invitation during which, essential parts of this work have been completed. ChS was partially supported by the “Karatheodory” post-doctoral fellowship No 2794 of the University of Patras. All image processing work has been done with ESO-MIDAS. DSS was produced at STSI under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope.


format: text

text:
# Dynamical Criterion for Quantum Chaos : Entropy Production in Subsystems

## Abstract

We examine the conjecture that entropy production in subsystems of a given system can be used as a dynamical criterion for quantum chaos in the latter. Numerical results are presented for finite dimensional spin systems as also for the quantum baker’s map. Of especial importance is the power spectrum of the entropy production which gets progressively more and more broad-banded as the degree of correlation in the Hamiltonian matrix is made to decrease.

More than twenty years after the discovery of dynamical localisation , possible signatures of quantum chaos are still under active exploration. While investigations have lent overwhelming support to a set of sharp conjectures and predictions relating to stationary aspects of the problem like level statistics and eigenvector distributions (see, e.g., ), corresponding dynamical criteria of quantum chaos are only recently beginning to take shape. One recalls that a simple-minded search for chaotic features in quantum time evolution of a system along lines similar those in classical evolution should yield negative results in consequence of the linearity of the Schrödinger equation. A more fruitful approach, on the other hand, would be to focus on reduced states of subsystems. Recent investigations have focussed on the reduced entropy of subsystems  as an important indicator in the context of quantum chaos. For a set of relevant results on reduced entropy as a measure of correlation in subsystems, see

The stationary criteria of level statistics and eigenvector distribution share the satisfactory feature that they are specifically quantal in nature, having no direct reference to the semiclassical limit. Dynamical criteria like the spreading of wavepackets examined in earlier investigations (see, e.g., ) on quantum chaos are not entirely satisfactory on this count since they refer to time scales pertaining to the semiclassical limit. By contrast, the dynamic criterion of entropy production in subsystems retains its relevance in the deep quantum regime as well as in the semiclassical one.

Quantum chaos has acquired a special relevance in the context of quantum information processing and quantum computation  in that it reflects in the manner in which decoherence and entanglement arise in a system of qbits. Entropy production can be used in this context as a natural measure of decoherence, and it seems appropriate to look at various features of entropy production in order to formulate a sufficiently satisfactory and universal dynamical criterion for quantum chaos. The degree of entanglement, which is represented by entropy production under certain circumstances, has also been shown to be sensitive to the presence of classical chaos. The present paper folows this line of approach and aims to point out that features of entropy production are relevant in characterising quantum chaos, one such useful feature being the time dependence of entropy production as revealed by its power spectrum.

More precisely, starting from an initial state $`\rho (0)`$ of a system $`S`$ made up of subsystems $`S_1`$ and $`S_2`$, we let the system evolve for a time $`t`$ and obtain the reduced state $`\rho _R(t)`$ of $`S_1`$ on taking partial trace over $`S_2`$. We then look at the time dependence of $`\rho _R`$. Even when one starts from a pure state $`\rho (0)`$, $`\rho _R`$ turns out more often than not to correspond to a mixed state and there takes place, in general, a decoherence with the passage of time. As mentioned above, one can measure the extent of decoherence by the reduced von Neumann entropy

$$s_R(t)=trace(\rho _R(t)ln\rho _R(t)).$$
(1)
We compare features of time dependence of $`s_R(t)`$ for systems conforming to the predictions of random matrix theory (RMT) with the corresponding features for a ‘regular’ (see below) syestm. We present evidence that while the former is characterised by a broad-band power spectrum, the latter involves sparsely distributed peaks, indicative of a quasiperiodic time dependence.The system we consider is the simple one of an interacting assembly ($`S`$) of $`n`$ spins or, in other words, an $`n`$-qbit system. The state space has dimension $`N=2^n`$, and the Hamiltonian determining the time evolution of $`S`$ is an $`N\times N`$ matrix. While quantum chaos in the true sense is expected to emerge only in the limit $`N\mathrm{}`$, we consider finite (and rather low) values of $`N`$ here so as to present in as simple a context as possible the features we aim to underline in this note. Finite dimensional considerations are, moreover relevant in numerous contexts relating to quantum computation.

Thus, we consider, first, an $`N\times N`$ real symmetric Hamiltonian matrix $`H_c`$ with randomly distributed matrix elements for which the stationary features predicted by RMT are conformed to and, next, a second Hamiltonian matrix $`H_r`$ that is equivalent to the Harper Hamiltonian on a torus. In the qbit representation the Hamiltonians correspond to strong and long-range interactions among the qbits. Alternatively, either of these can also be interpreted to represent the Hamiltonian of a chain of $`N`$ fermions of the form

$$a_{ij}(c_{i}^{}{}_{}{}^{+}c_j+c_ic_{j}^{}{}_{}{}^{+}),$$
(2)
where $`a_{ij}=a_{ji}`$, and $`c_i`$, $`c_{i}^{}{}_{}{}^{+}`$ denote fermion annihilation and creation operators. For $`H_c`$ , the coefficients $`a_{ij}`$ are chosen as randomly distributed matrix elements (subject to reality ad hermiticity) while, for $`H_r`$, we take

$$H_r=\gamma _1T+\gamma _2V.$$
(3)
Here $`\gamma _1`$ , $`\gamma _2`$ are real constants,

$$T_{i,i+1}=T_{i+1,i}=T_{1N}=T_{N1}=1/2(i=1,\mathrm{}N1),$$
(4)
and

$$V_{ij}=\mathrm{cos}(2\pi j/N)\delta _{ij}.$$
(5)
This represents a Harper’s Hamiltonian

$$H_r=\gamma _1\mathrm{cos}(2\pi p/P)+\gamma _2\mathrm{cos}(2\pi q/Q),$$
(6)
on a torus of area $`PQ`$ with $`N`$ basic states ($`\mathrm{}=PQ/2\pi N`$) and with torus boundary conditions $`\omega _1=\omega _2=0`$. The parameters $`\gamma _1,\gamma _2`$ have been set at (resp.) $`0.5,2.5`$ in our computations. In the fermion-chain interpretation, $`H_r`$ represents a Hamiltonian in the $`1`$-fermion subspace for the chain. However, in the following, we refer to $`H_c`$ and $`H_r`$ simply as Hamiltonians of an interacting spin system.

With the Hamiltonian of the system $`S`$ defined as above, we focus on a subsystem $`S_1`$ made up of, say, $`p(<n)`$ specified spins and look at the reduced density matrx (RDM) given by

$$\rho _R(t)=Tr^{(np)}[e^{\frac{iHt}{\mathrm{}}}\rho (0)e^{\frac{iHt}{\mathrm{}}}],$$
(7)
where $`H=H_c`$ or $`H_r`$ as the case may be, and where $`Tr^{(np)}`$ indicates partial trace with respect to states of the remaining $`(np)`$ number of spins. We then calculate $`s_R(t)`$ of Eq.(1). We present below numerical results for $`n=8`$ ($`N=256`$), and $`p=5`$, i.e., a subsystem made up of five spins. For convenience of presentation, we represent all the quantities as dimensionless; in other words, we use arbitrary but consistent units. Fig. 1 gives a comparison of the typical time-dependence of $`s_R(t)`$ for (a) $`H_c`$ and for (b) $`H_r`$ with the same value of $`\mathrm{}`$. For the sake of comparison we subtract a constant term from each matrix element of $`H_c`$, $`H_r`$ so as to reduce to zero the mean of all the matrix elements.

Referring to Fig. 1 (see also Fig. 6) we note the following distinctive features of $`s_R(t)`$ : (i) there is an initial ‘transient’ regime in which $`s_R(t)`$ increases more rapidly for $`H_c`$ (almost instantaneous in the scale of Fig. 1, but discernible in Fig. 6(a) for the quantum baker’s map) than for $`H_r`$; (ii) after the initial transient regime there arises a ‘steady’ state in which $`s_R(t)`$ fluctuates according to some invariant pattern; however, the mean value around which the fluctuation takes place is, in general, larger for $`H_c`$ than for $`H_r`$; (iii) the fluctuation is more erratic for $`H_c`$ than for $`H_r`$ \- while the latter appears to be quasiperiodic in nature, the former is random by comparison.

The contrast between the modes of fluctuation in the two situations under study is seen more clearly in Fig. 2 where we present the power spectrum of the fluctuations in $`s_R`$. One finds that the power spectrum for $`H_c`$ is a broad-band one while that for $`H_r`$ has a relatively smaller number of sharp peaks (the power spectra have been obtained with the help of a standard FFT package, and the symmetry about the mid-point of the frequency range has no physical relevance, being an artifact resulting from the program). It may be mentioned that the qualitative features indicated above persist against changes in the value of $`\mathrm{}`$ and the dimension of the subsystem considered. Details will be presented elsewher.

In order to set in perspective the distinction between $`H_r`$ and $`H_c`$ one can compute the nearest neighbour level spacing distribution for the two Hamiltonians. One observes in Fig. 3 that the level spacing distribution for $`H_c`$ is close to the Wigner distribution implied by RMT, while that for $`H_r`$ deviates to a large extent from the Wigner distribution and is, in fact, closer to the Poisson distribution followed by systems close to integrable classical systems. In other words, there exists a discernible correlation between stationary features implied by the RMT and dynamic features of entropy production in subsystems. Put differently, entropy production can indeed be taken as a valid dynamical criterion for quantum chaos - a criterion, moreover, that makes no direct reference to classical time evolution. While the latter feature is apparent from the manner in which the criterion has been formulated, one also finds support for this assertion from the fact that features of entropy production mentioned above appear even in situations with higher values of $`\mathrm{}`$ (not shown in our figures) compared to the value used here.

A quantum system represented by a Hamiltonian with randomly distributed matrix elements is, in a sense, analogous to classical systems exhibiting ‘hard’ chaos. By contrast, analogs to ‘regular’ classical systems would be ones with Hamiltonian matrices having correlated matrix elements. In between, there exist Hamiltonian matrices with matrix elements having various intermediate degrees of correlation. It would be interesting to see if features of entropy production in subsystems of such systems is, in some sense, intermediate between those of the systems considered above, analogous to the fact that systems characterised by ‘soft’ chaos in the semiclassical limit exhibit level spacing distributions intermediate between the Wigner and Poisson distributions(see e.g.,).

With this in view we modify our Hamiltonian $`H_r`$ by replacing a given fraction of matrix elements with randomly distributed ones. To be more precise, we construct a family of Hamiltonian matrices $`H(f)`$ depending on the parameter $`f(0<f<1)`$ where increasing values of $`f`$ correspond to greater degrees of correlation between the matrix elenments. $`H(f)`$ is obtained from $`H_r`$ by replacing $`(H_r)_{mn}`$ with randomly distributed numbers for $`m,n`$ satisfying $`|mn|>fN`$, where $`f`$ is some chosen fraction. Plots of $`s_R(t)`$ for decreasing values of $`f`$(not shown in our our figures) indicate that decreasing values of the degree of correlation indeed correspond to a gradual transition in features of entropy production from those of regular systems (Fig. 1(b)) to chaotic ones (Fig. 1(a)), being intermediate between the two. This is further corroborated from Figure 4 where we display the corresponding power spectra for two chosen values of $`f`$. One observes that with decreasing $`f`$ the power spectrum of entropy production becomes more and more broad-banded.

The fact that the dynamic criterion of chaos indicated in this paper including, in particular, the power spectrum of entropy production, is well correlated with the static criteria implied by RMT is borne out by looking at eigenvector distributions of $`H(f)`$ for various values of $`f`$.

The results are presented in Fig. 5 in which the residual parameters  for the eigenvectors of $`H(f)`$ are displayed, once again for a chosen pair of values of $`f`$. The residual parameter for any given eigenvector is defined as

$$r=\frac{1}{N}((p_kp_{k}^{}{}_{}{}^{(0)})^2)^{\frac{1}{2}},$$
(8)
where the set of numbers $`p_k`$ represent the probability distribution for the elements of the eigenvector and $`p_{k}^{}{}_{}{}^{(0)}`$ stands for a Gaussian distribution with the same mean and variance. Low values of $`r`$ for all the eigenvectors implies a correspondingly high degree of allegiance to the predictions of the RMT. One observes from the figures that decreasing values of $`f`$ indeed corresponds to an increasing degree of conformity to the stationary features implied by RMT.

The above features of reduced entropy for $`H_c`$ are also in evidence for the quantum baker’s map . Fig. 6 shows (a) the time-dependence of $`s_R`$, and (b) the power spectrum, which are of the general nature of Fig.s 1(a), 2(a) (as opposed to 1(b), 2(b)) respectively. Interesting details of entropy production for the quantum baker’s map, as also the quantum cat map , both known to be chaotic in the classical limit, will be presented elsewhere.

In summary, we have presented evidence that features of entropy production, including the power sprctrum in the steady state, in subsystems of a given system can be used as dynamical criteria of quantum chaos in the latter. Such dynamical criteria complement the stationary criteria discussed in the literature. The features presented show a progressive transition from regular to chaotic systems when the degree of correlation in the underlying Hamiltonian matrix is made to decrease.

While entropy production in subsystems can be loooked upon as an ‘internal’ dynamical criterion for chaos, dynamical features of an ‘external’ nature have also been discussed in the lterature in various related contexts, namely the distinctive features of chaotic systems in interaction with or perturbed by external ones.

Other features of entropy production in subsystems, including fractal characteristics, will be presented in a future publication.


text:
# The mechanical relaxation study of polycrystalline MgC1-xNi3

## I Acknowledgments

This work is supported by the National Science Foundation of China (Grant No. 0204G141) and by the Ministry of Science and Technology of China (NKBRSF-G19990646).

Figure Captions

Figure 1. Reduced resistance $`R/R_{300K}`$ of a superconducting ($`x0`$) and a non-superconducting ($`x0.2`$) MgC1-xNi3 samples vs temperature $`T(K)`$.

Figure 2. Internal friction ($`Q^1`$) and reduced modulus ($`\mathrm{\Delta }Y/Y_{81K}`$) of the mechanical relaxation spectra of a superconducting and a non-superconducting MgC1-xNi3 samples from liquid nitrogen temperature to room temperature.

Figure 3. The internal friction spectra ($`Q^1`$) of the P2 peak for the superconducting MgC1-xNi3 sample with the resonant frequency 2605 and 4098 Hz, respectively.

Figure 4. A fit of Cole-Cole law (solid line) to the P2 peak. The dashed line is the prediction of the Debye theory.


text:
# A hybrid model for chaotic front dynamics: From semiconductors to water tanks

## Introduction.

Moving fronts are the source of complex patterns in very different areas of physics CRO93 , chemistry KAP95a ; MIK94  and even biology SHE97 ; MUE02 . While the propagation of single fronts is well understood for many systems, only little is known on how front generation and annihilation processes or collisions between fronts influence the possible bifurcation scenarios.

For the sake of concreteness we illustrate our method with a semiconductor superlattice system, which is well studied experimentally and theoretically ESA74 ; KAS95 ; HOF96 ; PAT98 ; WAC02 ; BON02  and shows a wide range of complex front patterns (cf. Fig. 1) including chaos BUL95 ; LUO98b ; AMA02a . Similar front patterns occur in many other systems SCH01 , e.g., spatially continuous semiconductor models describing impurity impact ionization breakdown CAN01  or globally coupled heterogeneous catalytic reactions GRA93a . Our aim is to point out a striking connection of those front patterns with discrete hybrid models of $`n`$ tanks such as they occur, for instance, in industrial production processes CHA93a , and derive a simple description of the front generation and annihilation dynamics in terms of an $`n`$-tank model which allows for an analytical treatment of the bifurcation scenarios. In particular, an unusual chaotic scenario will be characterized by this generic model. We stress that the presented methods do not depend on the pecularities of the superlattice model, and are expected to be equally applicable to other front systems with global coupling.

In order to formulate a generic model for front propagation, we will first develop a front model based on the spatio-temporal dynamics in a semiconductor superlattice. The dynamical variables are the positions of the accumulation and depletion fronts which evolve with a velocity determined by the number of fronts present in the system. The number of fronts changes according to simple rules which mimic generation of fronts at the emitter and the annihilation of fronts. Such models are quite general and similar ideas have been sketched for front propagation in bulk GaAs CAN01 . It will be shown that this front model is identical to an $`n`$-tank switched arrival system CHA93a ; PET03 , if one identifies the spacing $`x_i`$ (see Fig. 1) between the leading depletion and the trailing accumulation front with the filling heights of the tanks. The tanks drain at a constant rate $`\mu `$, while one of the tanks is connected to a server and filled at a rate $`\lambda =\stackrel{~}{n}\mu `$ (Fig 2), where $`\stackrel{~}{n}`$ is the number of tanks which are not empty. In the present model the server will switch to a different tank, if this tank is empty and the filling level of the currently filled tank is larger than $`x_{min}`$. In the simplest case $`n=3`$, the tank–model transforms into a one–dimensional iterated map which is very similar to the well known tent–map. We will show that the front dynamics in semiconductor superlattices as simulated in AMA02a  can largely be explained by such a one–dimensional map. In particular, atypical border collision bifurcations occur BAN99 . The general type of models we consider here is called hybrid since they contain continuous as well as discrete state variables.

## Modeling the front dynamics.

The superlattice dynamics is described by the electron densities $`n_m(t)`$ in the quantum wells, labeled by $`m=1,2\mathrm{}N`$. The current density $`J_{mm+1}(n_m,n_{m+1},F_m)`$ between adjacent wells is evaluated by the sequential tunneling model, where $`F_m<0`$ is the electric field between these two wells WAC02 . A microscopic calculation of $`J_{mm+1}`$ yields an N–shaped current density vs. electric field characteristic as depicted in Fig. 3. In the electrically neutral regions of the superlattice, there are two stable states for fixed current density $`J`$, one at a low field $`F^l(J)`$ and one at a high field $`F^h(J)`$ (Fig 3). Since in the following $`J`$ is varied only in a small range, we make the approximation that both fields do not depend on the current density, and set $`F^l(J)=0`$, $`F^h(J)=F^h>0`$.

The emitter boundary condition (at well $`m=1`$) is characterized by its intersection point $`J_c`$ with the characteristics between neutral wells. For small $`J_c`$ propagating inhomogeneous field distributions alternating between $`F^l`$ and $`F^h`$ occur, see Fig. (1). For a transition from a low field domain $`F_m=0`$ to a high field domain $`F_n=F^h`$ $`(n>m)`$ a negatively charged electron accumulation front located at the domain boundary, is required by Poisson’s equation. Let us denote by $`a_i`$ the distance of its center of charge from the emitter. The index $`i=1\mathrm{}N_a`$ labels different accumulation fronts such that $`a_i<a_{i+1}`$, where $`N_a`$ is the total number of accumulation fronts in the system. Similarly, a transition from a high field domain to a low field domain can be attributed to an electron depletion front whose center of charge is located at $`d_i`$ with $`i=1\mathrm{}N_d`$, and $`N_d`$ the number of depletion fronts. Since we only allow for alternating accumulation and depletion fronts, we have $`N_dN_a\{1,0,1\}`$.

The applied voltage $`U`$ between emitter and collector defines the partial length $`L_h=U/F^h`$ of the superlattice which is at high field. This imposes a global constraint on the front positions by

$$L_h=\underset{i=1}{\overset{N_d}{}}d_i\underset{i=1}{\overset{N_a}{}}a_i\text{mod }L.$$
(1)
Here $`L`$ denotes the total length of the superlattice. The expression $`\text{mod }L`$ in (1) means that $`L`$ has to be added if $`a_{N_a}>d_{N_d}`$ such that $`L_h[0,L]`$.

The front velocities are given functions of the total current $`\dot{a}_i=v_a(J)`$, $`\dot{d}_i=v_d(J)`$, as shown in Refs. CAR00 ; AMA01 . Differentiating eq. (1) yields $`v_d(J)/v_a(J)=N_a/N_d`$. This determines the current $`J(N_a/N_d)`$, which is a monotonically increasing function, as $`v_d(J)`$ and $`v_a(J)`$ are increasing and decreasing, respectively AMA02a . Rescaling the time such that $`v_a+v_d=2`$ gives

$$v_a=\frac{2N_d}{N_a+N_d}v_d=\frac{2N_a}{N_a+N_d}$$
(2)
Front injection at the emitter is governed by $`J_c`$. For $`J<J_c`$ the region close to the emitter is pinned at a low field WAC02 . Hence an accumulation front is injected, if the preceding front is a depletion front. Due to the finite width and buildup time of the fronts new fronts can not be injected arbitrarily fast. Phenomenologically we therefore introduce a distance parameter $`p_h`$, and assume that front generation is suppressed while $`d_1<p_h`$, where $`d_1`$ is the position of the first depletion front. In the same fashion, a depletion front is injected if $`J`$ rises above $`J_c`$ and $`p_l<a_1`$, where $`p_l`$ is a parameter describing the minimum front distance for depletion front injection. As the current $`J`$ is a monotonic function of $`N_a/N_d`$, the condition $`JJ_c`$ can be transformed into $`N_a/N_dr_c`$ with a new critical parameter $`r_c`$.

The processes which reduce the number of fronts are collisions of two fronts of opposite polarity and the running out of fronts at the collector. Both processes affect $`N_a/N_d`$ and potentially trigger generation of a new front at the emitter.

Let us summarize the rules for front dynamics:

1. The positions of the accumulation fronts $`a_i`$ for $`i=1\mathrm{}N_a`$ and depletion fronts $`d_i`$ for $`i=1\mathrm{}N_d`$ evolve according to $`\dot{a}_i=v_a`$ and $`\dot{d}_i=v_d`$ with the velocities (2) until one of the following rules applies.
2. If $`N_a/N_d<r_c`$ and $`p_h<d_1<a_1`$ then increase $`N_a`$ by one, re-index $`a_ia_{i+1}`$ for all $`i`$ and set $`a_1=0`$ (injection of accumulation front).
3. If $`N_a/N_d>r_c`$ and $`p_l<a_1<d_1`$ then increase $`N_d`$ by one, re-index $`d_id_{i+1}`$ for all $`i`$ and set $`d_1=0`$ (injection of depletion front).
4. If $`a_i^{}=d_j^{}`$ for any $`i^{},j^{}`$ then decrease $`N_a`$ and $`N_d`$ by one, re-index $`a_{i+1}a_i`$ for $`ii^{}`$ and $`d_{j+1}d_j`$ for $`jj^{}`$ (annihilation of front pair).
5. If $`a_{N_a}>L`$ decrease $`N_a`$ by one (accumulation front hits collector).
6. If $`d_{N_d}>L`$ decrease $`N_d`$ by one (depletion front hits collector).

Note that the voltage parameter $`L_h`$ only enters in the initial condition for the front positions (see Eq. (1)).

In the following we will restrict ourselves to the case $`r_c<1`$ and $`p_l=0`$. Rule FII does not apply as long as $`N_a>r_c/(1r_c)`$ and $`N_a`$ can then only decrease. Consequently $`r_c`$ imposes the following limits on the number of fronts:

$`N_dn;N_a`$ $``$ $`n1;`$ (3)
where $`n`$ is the largest integer less than $`1/(1r_c)+1`$. The injection of an accumulation front is immediately followed by the injection of a depletion front (FIII). This detaches a high field domain from the emitter and opens a new one. Furthermore we restrict ourselves to the case where no fronts reach the collector (i.e. $`L>nL_h`$). This corresponds to the situation in Fig. (1) and we have $`N_a=N_d1`$.

## The n-tank model.

Under the restrictions given above we may choose the lengths of the high field domains (cf. inset of Fig.1) $`x_1=d_1`$ and $`x_i=d_ia_{i1}`$ for $`i=2\mathrm{}N_d`$ as new dynamical variables governed by the rules:

1. The lengths evolve according to $`\dot{x}_i=\mu +\lambda \delta _{i1}`$ with $`\mu =2/(2N_d1)`$ and $`\lambda =N_d\mu `$.
2. If $`N_d<n`$ and $`x_1>p_h`$ then increase $`N_d`$ by one, re-index $`x_ix_{i+1}`$ for all $`i`$ and set $`x_1=0`$.
3. If $`x_i^{}=0`$ then decrease $`N_d`$ by one, re-index $`x_{i+1}x_i`$ for all $`ii^{}`$.

The voltage constraint (1) is expressed in the new variables by $`_kx_k=L_h`$.

It is very illustrative to interpret this model in the language of an $`n`$-tank switched arrival system (cf. Fig.2). A domain may be viewed as a tank with fluid content $`x_i`$. One tank is filled at rate $`(N_d1)\mu `$, while $`N_d1`$ tanks are drained at rate $`\mu `$, and $`nN_d`$ tanks are empty and waiting for filling, because the currently filled tank has not yet reached the minimum filling level $`x_{min}=p_h`$.

Let us first discuss the limiting case $`p_l=p_h=0`$. If a tank is empty at time $`t_m`$ (i.e. $`x_k(t_m)=0`$), the rules TII and TIII apply instantaneously and the system starts to fill this tank ($`x_1(t_m^+)=0`$). In this case no tank has to wait for filling. Apart from the re-indexing procedures we obtain a switched arrival system CHA93a . As shown in SCH95h  the system is chaotic for all $`n>2`$ and has a constant invariant probability measure. By sampling the dynamics at the switching times $`t_m`$ we obtain the Poincaré map $`P:\stackrel{}{x}(t_m)\stackrel{}{x}(t_{m+1})`$:

$$\stackrel{}{x}(t_{m+1})=\stackrel{}{x}(t_m^+)+\dot{\stackrel{}{x}}\mathrm{\Delta }t_m$$
(4)
where $`\dot{\stackrel{}{x}}=(\lambda \mu ,\mu ,\mathrm{},\mu )`$ is a constant vector and $`\mathrm{\Delta }t_m=t_{m+1}t_m=\mathrm{min}_{j1}(x_j(t_m^+)/\mu )`$. Generally this is an $`n`$ dimensional mapping which acts due to the normalization $`_ix_i=L_h`$ on the $`n2`$ dimensional boundary of a regular $`n`$-simplex embedded in $`^n`$. For $`n=3`$ we can reduce this map to a one-dimensional map of the interval $`[0,L_h]`$ onto itself: From the rules we obtain $`x_1(t_m^+)=0`$, $`x_2(t_m^+)=x_1(t_m)`$, and $`x_3(t_m^+)=L_hx_1(t_m)`$. If $`x_1(t_m)<L_h/2`$, the tank $`2`$ empties after $`\mathrm{\Delta }t_m=x_1(t_m)/\mu `$ and we find $`x_1(t_{m+1})=2x_1(t_m)`$. Otherwise, if $`x_1(t_m)>L_h/2`$, the tank $`3`$ empties after $`\mathrm{\Delta }t_m=(L_hx_1(t_m))/\mu `$ and we find $`x_1(t_{m+1})=2(L_hx_1(t_m))`$. Thus the dynamics is that of the well known tent map: $`x_1(t_{m+1})=L_h2|x_1(t_m)L_h/2|`$.

Now we consider the case $`n=3`$ with a finite minimum filling level $`p_h>0`$. Since by TII $`x_1`$ must be larger than $`p_h`$ for switching, the Poincaré map is now only defined on the interval $`[p_h,L_h]`$. For $`x_1[p_h,L_hp_h/2]`$ the map remains unchanged from the tent map discussed before. If the switch to an empty tank occurs at $`x_1(t_m)>L_hp_h/2`$, the third tank will be empty before the newly filled tank has reached the minimum filling level $`p_h`$. Switching is therefore suppressed until $`x_1(t_{m+1})=p_h`$. Thus for $`n=3`$ we obtain (cf. inset of Fig.4):

$$x_1(t_{m+1})=\mathrm{Max}\{(L_h2|x_1(t_m)L_h/2|),p_h\}$$
(5)
The bifurcation structure of this map with respect to $`L_h`$ is unusual (Fig.4) but has similarities to the flat-topped tent map GLA94 ; WAG02 . Due to the horizontal segments the map cannot show chaotic behaviour although it undergoes several bifurcation cascades. A trajectory on the tent map segments would explore the attractor of the tent map and will therefore either fall on one of the flat segments or stay at an unstable periodic orbit contained in the remaining part of the tent map. Once arrived at one of the flat segments the trajectory continues on a stable periodic orbit. By analyzing the $`k`$-th return map of the modified tent map we find that windows of period $`k`$ are found inside the intervals $`L_h/p_h(\frac{2^k1}{4j2},\frac{2^k+1}{4j2})`$ for $`jI_k\{j|j2^{k2}1\}`$. Thus higher and higher stable periodic orbits are arising in smaller and smaller intervals.

## Discussion.

In order to compare the bifurcation diagram of the modified tent map (Fig.4) with the bifurcation diagram of the superlattice model (inset of Fig. 3), we note that the voltage $`U`$ and the position of front annihilation in the superlattice model correspond to the parameters $`L_h`$ and $`x_1(t_m)`$ in the map, respectively. We can therefore identify the periodic windows found in the Poincaré map with the periodic pattern observed in the front model. We note that the bifurcation scenarios agree very well up to voltages of about $`1.1`$V, where a period-3 window occurs (inset of Fig. 3). In particular, the remarkable spider-like structure at $`0.9`$V is found for $`L_h/p_h=3`$ in the bifurcation diagram of the map. The reason for this spider is that the contracting flat region is mapped onto the unstable fixed point of the map. Further prominent features in both bifurcation diagrams are the horizontal lines given by the first, and higher iterated images of the minimal filling height. Thus the character of the bifurcation scenario of complex front dynamics is explained by considering a one-dimensional iterated map with only one free parameter $`L_h/p_h`$.

The minimum filling height corresponds to the minimum domain width. For a real superlattice system this quantity will slightly depend on the history of the dynamical evolution. Therefore the slope of the flat segments is not exactly zero but attains a small finite value $`\alpha `$. In this case a period-$`k`$ orbit is unstable if $`|2^k\alpha |>1`$. Close to the spider-like points where periodic orbits with large period $`k`$ occur, we can then find a finite chaotic regime, as observed in our simulations with the microscopic model.

## Conclusions.

We have introduced a hybrid tank model as a simplified model for front dynamics in nonlinear systems with a global constraint. This simplified model provides an explanation for the complex accumulation and depletion front dynamics in semiconductor superlattices. In this unexpected analogy, the tank filling levels correspond to the widths of the high field domains in the superlattice, while the position of the server determines which high field domain is connected to the emitter. Although this analysis has started from a superlattice model, the presented methods are very general and are expected to work successfully for other systems which provide a global coupling similar to Eq. (1).

## Acknowledgments.

This work was partially supported by DFG in the framework of Sfb 555 and the Volkswagen Foundation (Grant No. I/76 279 - 280). We thank N. Janson, A. Pikovsky and J. Schlesner for stimulating discussions.


text:
# Limit current density in 2D metallic granular packings

## 1 Introduction

Among the wide panel of physical properties exhibited by granular materials, the electrical properties are rather spectacular since they are the origin of radiotelecommunication. Let us remind that the first reported antenna has been found by Branly in 1891 branly ; this antenna was essentially composed of a glass tube containing metallic filings. It is not yet understood why a packing of metallic filings becomes a better conductor when an electromagnetic wave hits the system. In the 19th century, Calzecchi-Onesti has also found that a similar insulator-conductor transition occurs when a high current flows through such a system onesti . Thermal analysis vdb  has shown that the Calzecchi-Onesti transition is due to the soldering of grains. The authors have observed a strong increase of the temperature along certain electrical pathes. This indicates the welding of the beads at the contacts. The pathes being better conducting shunt the remain of the packing and an electric transition occurs when a continuous weakly resistive path is established. On the other hand, Bonamy thesisbon  has found that the $`iV`$ curve of one contact between two beads is continous : no sharp transition occurs. Finally, in apl , it was shown that the $`iV`$ diagram of a lead bead heap is a combinaison of an exponential behavior of the voltage due to slow microsoldering of the beads followed by a sharp transition as soon as a critical current is reached. This behavior is only obtained in specific conditions : (i) the filings have to be loose packed and (ii) to be a little oxydized at the grain surface. Those conditions are rather fuzzy and need quantification. Experiments should be proposed to obtain a complete picture of the electrical behaviors of the granular matter.

It is well known that a 3D granular packing is characterized by a complex system of arches along which force lines take place. A good way for investigating this network is to study how an electrical current $`i`$ can flow through the system. Indeed, Hertz hertz  observed that electrical contacts are enhanced by pressure; electrical currents tends to follow the paths which are created by the lines of mechanical forces. From this point of view, the electrical paths form an image of the force chains.

In a previous paper apl , a 3D granular packing has been studied. The voltage $`V`$ has been recorded as a function of the injected current intensity $`i`$. High memory effect has been pointed out and two regimes, i.e. a conducting and a non-conducting one, have been found to occur depending on the electrical history of the system. The transition between both regimes is sharp and the so-called Calzecchi-Onesti transition has been found. This transition was interpreted in term of the dielectric breakdown of the oxyde layer as in semiconductor devices mich . The voltage $`V`$ behavior versus the current $`i`$ was found to be well fitted by an exponential law like

$$V=V_0(1\mathrm{exp}(i/i_0))$$
(1)
where $`V_0`$ is the saturation voltage obtained at high injected current and $`(V_0/i_0)`$ is the initial resistance of the packing.

Several questions rise from those experiments especially about the origin of the reversibility of the transition. The microsoldering of beads is surely expected to be an irreversible process, thus a strong memory effect behavior should be highlighted.

This paper concerns further investigation of the $`iV`$ diagram of a 2D metallic grains packing. For that purpose, a simple system has been considered. It consists in a packing of 1 mm thick planar pentagons for which edges are 10 mm wide (Fig.1). Such a 2D system has been studied in ref.epje ; physA  as far as compaction properties are concerned. It has been shown that the formation of arches inside the packing is responsible for the anisotropy of the force lines in granular packings. The impossibility for pentagons to be packed into in a compact configuration, whence form intrinsic arches, is thus very interesting to study.

In the present work, we report the analysis of the voltage behavior during current cycles. This allows us to describe the $`iV`$ plane as supporting two distinct regimes. Contrary to bulk materials, granular packing is not only represented by one single curve in the $`iV`$ plane but by several trajectories depending on the history of the injected currents in the system. Memory effects have been pointed out. The points of the $`iV`$ plane where such a granular system may exist define an area limited by the so-called irreversibility line and the intrinsic resistance line. The irreversibility line is a line along which the electrical resistance is weakened when the current increases.

In Sect. 2, the experimental setup and conditions are described. The next sections are devoted to results and interpretations. Conclusions are drawn in Sect. 5.

## 2 Experimental setup

About 200 aluminium pentagons are disposed in a plexiglas vertical vessel of $`210\times 200`$ mm2. In order to avoid overlaps between pentagons and, in so doing, to keep 2D conditions, the vessel width is 1 mm thick (Figure 1). The packing is built by droping the pentagons one by one between vertical planes. No further compaction process has been performed. The electrical contacts are set on two particular pentagons. Those contacts are soldered and the measure of the resistance is performed in a two-wires configuration : the lowest measured global resistance of the packing is about 20 $`\mathrm{\Omega }`$ which is much higher than that of the current leads.

In order to avoid electromagnetic perturbations, the vessel and the measurements instruments are placed in a Faraday cage. Electrical contacts are connected by coaxial cables to a Keithley K2400 current source which allows to measure the voltage. During the experiment, we chose to modify the current intensity $`i`$ and to measure the voltage $`V`$.

## 3 Results

The injected current is set to $`i=sign10^q`$ A, $`q`$ being the running parameter in the following description. A maximum current value is set by fixing $`q_{max}`$. Starting from the minimum current $`q=6`$, the current is increased up to $`q=q_{max}`$. It is then decreased to $`q=6`$ and to “negative” values. To sum up, a cycle is defined as $`q`$ varies as follow

$`sign=1:`$ $`q=6q_{max}`$
$`q=q_{max}6`$
$`sign=1:`$ $`q=6q_{max}`$ (2)
$`q=q_{max}6`$
such that, a cycle is achieved. A new value for $`q_{max}`$ is fixed and another cycle is made. The maximum values have been successively set to 4, 3, 2.5, 2, 1.5 and finally 1 which corresponds to the maximum current intensity that the source can deliver (100 mA). In Figure 2a, 2b and 3, the same recorded $`iV`$ curves are shown. Since the range of current intensities is broad, Figure 2b has been represented in a log-log plot.

In order to describe a complete cycle, let us follow the arrows indicated on the Figure 2a, 2b and 3. Starting from point A, while the current increases, the voltage increases following the law (materialized by a solid line in Figure 2a and 2b)

$$V=\left[V_0(1\mathrm{exp}(i/i_0))\right]+R_ri$$
(3)
where $`R_r`$ is the intrinsic electrical resistance of the packing. The system reaches the point B (10-1.5,1.6). The current is then decreased (arrow 2); the voltage decreases following a linear law $`V=R_ai`$ where $`R_a`$ is a fitting coefficient, till the point B’ (-10-1.5, -1.6) (Figure 3). Next, the current is increased till 10-1 A (arrows 3). The voltage first follows the same linear law till point B. After that, it follows the law of Eq.(2) till point C (10-1,2.4). The injected current $`i`$ is then decreased (arrows 4) : the voltage linearly decreases as $`V=R_a^{}i`$, where $`R_a^{}`$ is lower than $`R_a`$ and follow the same linear law as the current is increased again to $`10^6`$ A (arrow 5).

A close look at Figure 3 shows that the feature still exists for low injected currents (see Figure 4). The electrical resistance is very high : about 25 k$`\mathrm{\Omega }`$; while the structure of the $`iV`$ curve is similar to Fig.2(a). This feature would not be observed whether the starting injected current is higher than in this example. That shows the high hysteretic behavior of such system and that different scales are in presence in a granular packing. Moreover a very small jump occurs at $`i=10^4`$ A. This should be attributed to a breakdown of the reverse diode as observed in a 3D system apl . This breakdown will be explained herebelow.

## 4 Interpretations

In order to summarize our results, Fig. 5 is a sketch of Figure 2a. The irreversibility line is described by Eq.(3). It is thought that microsoldering process occurs along this line since electrical arcs can be seen during the experiment. The intrinsic resistance line is obtained when the contacts can be neglected. The resistance of a aluminium wire, representative of the soldered grains, would be then obtained. Those two lines define in the $`iV`$ plane an area included between both lines. This zone localizes the accessible points by the system. Ohm’s law is of application as indicated by the slope of straight lines in Figure 2b. Let us consider the point ($`i_1`$,$`v_1`$) from this zone. The straight line which links this point to the origin is given by $`V=(v_1/i_1)i`$ (Ohm’s law). This line intersects the curve of Eq.(2) at ($`i_c`$,$`v_2`$). The point ($`i_1`$,$`v_1`$) can be thus reached by increasing the current following the irreversibility line to the point ($`i_c`$,$`v_2`$) and then by decreasing the current following $`V=(v_1/i_1)i`$. That also means that the Ohm law is applicable in the range of the injected currents \[$`i_c`$,$`i_c`$\]. In this range, the electrical resistance behaves as a conducting bulk material.

The different behaviors met in the $`iV`$ plane find their origin in the intergrain nature. A large number of publications concern the modelization of contact network. They consider either a resistors network last ; kirk ; takayasu ; dearcangelis ; roux ; pennetta  or diodes network hinrichsen ; inui  or mixing resistors and diodes network redner ; bigalke . Some of those models may reproduce the dielectric breakdown and consider non-linear resistance. Considering the chemical composition and the geometry of a contact, i.e. an oxyde layer between two metallic pentagons, the non linearity of the electrical resistance can be introduced by considering a granular packing as a network of diodes. A contact is then represented by two opposite diodes in parallel. In order to be complete, irreversible processes have to be introduced in the model since memory effect behavior occurs in the system.

We propose to model contacts as if they are made by two opposite diodes in parallel plus a resistor $`R`$ in parallel. Diodes have a behavior described by Eq.(3). Resistors are characterized by a section $`S`$ and a resistivity $`\rho `$ and follow Ohm’s law. The current will preferencially flow through either the diode or the resistance according to the relative resistance of both components.

Consider the system as being a mean-field media characterised by an average resistivity $`\rho `$. The resistance $`R`$ of the contact can be written using Pouillet law as $`R=\rho \mathrm{}/S`$ where $`\mathrm{}`$ is the effective length of an electrical path between the electrodes. When $`V=V_0`$ at the saturation, it can be rewritten according to Pouillet law : $`V_0=\rho \mathrm{}i/S`$. The ratio $`i/S`$ equals to the conductivity $`\sigma `$ multiplied by the limit electrical field $`E_0=V_0/\mathrm{}`$ (about 5 V/m). In other words, $`i/S`$ equals to $`J_{\mathrm{}}`$ which is hereby called a limit current density. This quantity remains constant in the range of high currents. This means that an increase of current intensity corresponds to an increase of the section of the contact resistor so as to keep $`J_{\mathrm{}}`$ constant. The grains knit together following the rule $`S=J_{\mathrm{}}/i`$. Within this approximation, the section increases of a factor 3 between the point B and C in Fig.2a. This microsoldering process is of course irreversible. When the current intensity decreases, the current follows the path defined by those enhanced contacts which shunts the diodes. The Ohm’s law is from then on application as shown in Fig. 2a and 5.

The very low current feature shown in Fig.4 exhibits rather the same behavior as described here above. When the same arguments are applied to this, a larger $`J_{\mathrm{}}`$ is found since $`E_0=2.5`$ V/m in this case. Nevertheless, a breakdown occurs at 10-4 A (Fig.4). This is similar to the one in a 3D system as in ref.apl . It can be explained as a Zener or avalanche breakdown of the reverse diode when the voltage across a contact reaches the critical value, i.e. the Calzecchi-Onesti transition.

## 5 Conclusions

The above measurements show that a weak electrical contact in granular metallic packings has to be modelized as three components in parallel : two opposite diodes and a resistor. This resistor modelizes the microsoldering and is characterized by a limit current density. When this limit density is reached, the equivalent section of the resistor increases, the grains weld together. This is an irreversible process. The model allows to explain the hysteretic loops in $`iV`$ curves. Moreover electrical properties of metallic grains packing are strongly determined by the electrical history of the system. Fine structures can be pointed out or not according to injected current cycles.

## Acknowledgements

S.D. is a F.N.R.S Scientific Research Worker. We would like to thank Prof. H.W. Vanderschueren for the use of MIEL facilities. We want also to thank Dr. M. Houssa (U. Marseille 3) and Dr. R. Goffaux (SEO, Luxemburg) for comments.


text:
# Pulse Shape Discrimination in the IGEX Experiment

## 1 Introduction

The nuclear Double Beta Decay (DBD) is an unique laboratory to investigate the nature and properties of the neutrino . The neutrinoless decay mode, if it exists, would provide an unambiguous evidence of the Majorana nature of the neutrino, its non-zero mass, and the non-conservation of lepton number. After the definitive confirmation that neutrinos have indeed non-zero mass, as the solar and atmospheric neutrino oscillation results imply , the neutrinoless Double Beta Decay (DBD) has become a most relevant subject of research because it is a process able to provide, in a relatively short time, the neutrino mass scale and its hierarchy pattern. The current best sensitivity limits of the effective Majorana electron neutrino mass -derived from the neutrinoless half-life lower bound- stands around $`m_\nu 3001300`$ meV  (the dispersion range is due to the uncertainties in the evaluation of the nuclear matrix elements), whereas most of the neutrino mass models indicate that the mass region where the Majorana nature of the neutrino would be resolved lies two orders of magnitude below ($``$ a few meV). To achieve the sensitivity required for such new objectives, it will require a large number of double beta emitter nuclei, a very low background and a sharp energy resolution in the Q-value region, and/or methods to disentangle signal from noise.

A typical example of this type of search is IGEX. The International Germanium EXperiment (IGEX) has been running in the Canfranc Underground Laboratory at a depth of 2450 m.w.e. in a search for the neutrinoless Double Beta Decay. Details and results of the experiment can be found in ref. . Three Germanium detectors (RG1, RG2 and RG3), of $``$2 kg each, enriched to 86% in 76Ge were used.

This paper presents the development and application of one technique used to reduce part of the radioactive background by discriminating it from the expected signal by comparison of the shape of the pulses (PSD) of both types of events. As an example, the method has been applied to the data recorded by two Ge detectors of the IGEX Collaboration, which has produced one of the two best current sensitivity limits for the Majorana neutrino mass parameter .

The rationale for PSD is quite simple: in large intrinsic Ge detectors, the charge carriers take 300 - 500 ns to reach their respective electrodes. These drift times are long enough for the current pulses to be recorded at a sufficient sampling rate. The current pulse contributions from electrons and holes are displacement currents, and therefore dependent on their instantaneous velocities and locations. Accordingly, events occurring at a single site ($`\beta \beta `$-decay events for example) have associated current pulse characteristics which reflect the position in the crystal where the event occurred. More importantly, these single-site events (SSE) frequently have pulse shapes that differ significantly from those due to the background events that produce electron-hole pairs at several sites by multi-Compton-scattering process, for example (the so-called Multi-Site Events (MSE)). Consequently, pulse-shape analysis can be used to distinguish between these two types of energy depositions: DBD events belong to the SSE class of events and will deposit energy at a single site in the detector while most of the background events belong to the MSE class of events and will deposit energy at several sites.

The IGEX detectors have modified preamplifier electronics to route and record the current pulses at a very early stage of preamplification, thus producing unique high-bandwidth pulse shape signals. Furthermore, to develop PSD techniques it would be highly desirable to obtain an earlier signal, even before it passes through the few unavoidable electronic components at the first stage of the detector preamplifier, resembling as much as possible the displacement current of the detector. This allows the development of algorithms that do not depend strongly on the preamplifier electronics in use. To this end, the transfer function of the preamplifier and associated front-end stage has been measured for each detector. This allows the reconstruction of the displacement current and an easy comparison to computed pulse shapes.

The paper is organized as follows: in Section 2, the properties of the digitized experimental pulses are shown and the work performed to understand and reproduce their shapes is presented. In Section 3, the PSD method is described. Finally, Section 4 displays the results of this analysis for the IGEX data.

## 2 Experimental and calculated pulses

The data acquisition system of the IGEX experiment is based on standard NIM electronics, each Ge detector having an independent electronic chain. Preamplifiers were modified for pulse shape analysis and each preamplifier fast-pulse output is routed to a LeCroy 9362 digital oscilloscope (800 MHz analog bandwith). The digitized pulse signal covers a total time of 1 $`\mu `$s using 500 points; it is worth noting that the time resolution, of about 100 ns (as inferred from the width of the peaked features or the fall time), limits the ability to resolve nearby features in the pulse such as lobes or discontinuities characteristic of a multiple-site interaction signal. Figure 1 shows the main features of the digitized pulses. These output pulses are taken at the very first stage of the amplification chain, but even so, there is an unavoidable instrumental distortion due to the preamplifier. This has been studied to determine its transfer function $`h(t)`$. To take into account this distortion in the pulse shape analysis, either the calculated pulse $`i(t)`$ is folded with the transfer function,

$$o(t)=_{\mathrm{}}^+\mathrm{}i(\tau )h(\tau t)𝑑\tau $$
(1)
or the experimental signal unfolded. The latter allows the recovery of some information lost because of the instrumental distortion . The transfer function of the preamplifiers, $`h(t)`$, depicted in figure 2, has been directly measured as the response of the preamplifier for a narrow $`\delta `$-like signal. Studies were also made following several methods including analog simulation of the preamplifier circuits and the analysis of the shapes of selected populations of experimental pulses. It was observed that the shape of the rise and the fall of the folded pulses is the same that the shape of the transfer function for those events in which the energy is released in some particular regions of the crystal; in particular, the left side of the transfer function can be deduced by studying pulses of events in the inner and lower part of the crystal, while the right side is derived from those produced in the outer and upper region. Measurements and estimates for the transfer function are found to be in quite good agreement.

The pulse shapes of the output signals can be reproduced numerically. An energy deposition in a Ge crystal produces a proportional number of electron-hole pairs, which move towards the electrodes. The induced current $`i`$, and consequently the electric pulse taken from the detector, is the sum of the contributions due to each type of charge carriers:

$$i(t)=i_\mathrm{e}(t)+i_\mathrm{h}(t)$$
(2)
The current is calculated as :

$$i_{\mathrm{e}(\mathrm{h})}=q_{\mathrm{e}(\mathrm{h})}\stackrel{}{E}_\mathrm{w}\stackrel{}{v_{\mathrm{e}(\mathrm{h})}}(\stackrel{}{E})$$
(3)
where $`\stackrel{}{E}_\mathrm{w}`$ is the so-called weighting field444The weighting field is the field that would appear in the crystal if a unity voltage was applied to the sensor electrode (assuming no impurity in the crystal). In principle, the drift velocity $`v`$ of charge carriers is proportional to the electric field $`E`$: $`v(r)=\mu E(r)`$, where $`\mu `$ is the mobility in the material medium. However, at high values of the electric field, a saturation velocity is reached. The dependence of the velocity on the electric field has been obtained empirically :

$$v(E)=\frac{\mu E}{(1+(\frac{E}{E_{\mathrm{sat}}})^\beta )^{1/\beta }}$$
(4)
The values commonly used for the parameters involved in expression 4 (the mobility $`\mu `$, the numerical parameter $`\beta `$ and the electric flield for saturation $`E_{\mathrm{sat}}`$) are summarized in table 1 .

The solution for the electric field ($`\stackrel{}{E}(\stackrel{}{r})=\stackrel{}{}\varphi (\stackrel{}{r})`$) in the crystal is derived from the Poisson equation, applied for the depleted and non-depleted regions:

$$^2\varphi (\stackrel{}{r})=\frac{\rho (\stackrel{}{r})}{ϵ},\mathrm{depleted}$$
(5)
$$^2\varphi (\stackrel{}{r})=0,\mathrm{non}\mathrm{depleted}$$
(6)
$$\varphi (\stackrel{}{r}_{\mathrm{i}nt})\varphi (\stackrel{}{r}_{\mathrm{e}xt})=V_0$$
(7)
The electric field depends on the geometry of the crystal (true coaxial or closed-end), the supplied voltage $`V_0`$ and on the residual space charge density $`\rho `$. The parameters $`r_{\mathrm{int}}`$ and $`r_{\mathrm{ext}}`$ correspond to the locations of the internal and external electrodes of the detector. A solution (depicted in figure 3) has been obtained for the cases of depleted and non-fully depleted, closed-end crystals of the IGEX experiment by using an iterative method of calculation. For a non-fully depleted detector, the boundary between the depleted and non-depleted regions in the crystal is unknown a priori. An invalid solution is obtained when using the boundary condition of eq. (7), because a region with a field having inverse direction appears (which means that this region is not depleted). Therefore, the boundary condition is changed by using the deduced limit for the depleted region instead of $`r_{\mathrm{int}}`$; then, a new solution is found for the electric field and, consequently, a new limit for the non-depleted region is obtained. This procedure is repeated until the solution does not change significantly between two successive iterations.

Once the electric field is known, the pulse shapes can be calculated. Some examples for SSE are shown in figure 4. Different radial positions have been considered in each plot for several vertical coordinates; the effect of the instrumental distortion is shown on the right plots. To reproduce MSE pulses, the shapes due to each individual interaction should be properly added weighted by the fractional amount of energy released in the crystal, obtained by Monte Carlo simulation.

## 3 PSD Method

The PSD method we have used consists in counting the number of lobes of the pulses and rejecting those events having more than two significant lobes or peaks. A SSE pulse is expected to have at most two lobes, one due to electrons and the other due to holes. Experimental pulses are first unfolded using the transfer function of the preamplifier. Then, to detect lobes a ”mexican-hat” filter $`F`$ of the proper width $`\sigma `$ is applied to the pulse. In fact, this filter is the second derivative of a gaussian:

$$F(t)=\frac{\sigma ^2t^2}{\sigma ^4}\times \mathrm{exp}(\frac{t^2}{2\sigma ^2})$$
(8)
and the filtered signal has a null mean value where there is no lobe in the original signal and a peak where a lobe is present. Therefore, it is straightforward to reject all the events having more than two lobes. Figure 5 shows the results of applying this method to four different pulses.

This robust method is nearly model-independent. Its effectiveness has been evaluated on calculated SSE. A test population of 2000 SSE pulses was generated for each detector. The locations were randomly chosen, uniformly distributed in the volume of the crystals. Each calculated pulse was folded with the proper preamplifier transfer function, then scaled to unit height, and a variable amount of gaussian noise was added, to reproduce experimental pulses of different energies and noise levels. Finally, the number of lobes of each pulse was obtained by applying ”mexican-hat” filters derived from gaussians of different widths. After the analysis, the fraction of pulses with three o more lobes was retained. This is an estimation of the expected number of misidentified SSE pulses with this technique. The results are shown in figure 6. When applied to a mixed set of SSE and MSE pulses, to obtain the best discrimination results, a narrow filter should be preferred over a broad filter, but it would produce a large fraction of misidentifications. Conversely, a safer, broader filter would not find as many lobes as a narrower filter, and its discrimination power would be smaller. As a compromise, we choose the filter with a characteristic width of 60 ns, thus keeping the misidentification error for the calculated SSE pulses under 5% for all pulses in the whole noise range considered. Notice that those MSE in which the separation of the interactions in the (r,z) plane is too small or the amounts of energy deposited at each point are extremely unbalanced will not be rejected by this technique.

This method of counting the number of lobes of the pulses has been applied also to a 22Na calibration spectrum and to a set of data taken following a large intrusion of radon in the shielding. Fig. 7 illustrates the reduction in the case of 22Na spectrum for detector RG2 and Fig. 8 the case of radon for detector RG3. A comparison of the results of the method for the cases of background, 22Na calibration and radon is shown in Fig. 9.

## 4 Results

The background reduction technique described above has been applied to the IGEX data whose pulse shape was recorded (those events whose pulse shape was not available are conservatively considered SSE). Table 2 summarizes the results (exposure, background levels in the region of interest from 2 to 2.5 MeV before and after the PSD, rejection factor) for each detector. Figure 10 shows the spectra before and after application of the PSD technique to the data of detectors RG2 and RG3. This method results in an efficient rejection leading to a background level (in the best case) of 0.10 c/(keV kg y) in detector RG2, as can be seen in Table 2. This value can be considered as the background limit achieved with the present PSD technique. The overall final background level of the set of detectors RG2 and RG3 together turns out to be of 0.10 c/(keV kg y), in the region of interest.

The PSD analysis has been applied to only 52.51 mole y out of the total 116.75 mole y (8.87 kg y) accumulated in the IGEX experiment. In Table 3 the IGEX data correponding to 8.87 kg.y in 76Ge in the region between 2020 and 2060 keV, in 2-keV bins, are given, with and without application of PSD (see also ref. ). The obtained half-life lower bounds are $`T_{1/2}^{0\nu }1.13\times 10^{25}y`$ for the complete data set and of $`T_{1/2}^{0\nu }1.57\times 10^{25}y`$ for the complete data set with application of PSD to 52.51 mole y. Accordingly, the upper limits on the neutrino mass parameter are 0.38–1.55 eV for the first data set and 0.33–1.31 eV for the second data set . The uncertainties originate from the spread in the values of the calculated nuclear matrix elements.

## 5 Conclusions

A Pulse Shape Discrimination technique to reject the radioactive background in the region in which the double beta decay signal is expected has been developed and applied to the data collected in the IGEX experiment, searching for the neutrinoless double beta decay of 76Ge. A satisfactory understanding of the pulse shapes has been achieved. The method described in this paper is based on the counting of the number of lobes of the pulses, using a proper filter. It has provided a rejection of $``$ 60 % of the events in the region of interest, accepting the criterion that those events having more than two lobes cannot be due to a double beta decay. Accordingly, the improved background levels provided by the PSD technique have allowed the improvement of the limits for the half-life of 76Ge and consequently, the effective electron neutrino mass bound.

## 6 Acknowledgments

The Canfranc Underground Laboratory is operated by the University of Zaragoza under contract No. FPA2001-2437. This research was partially founded by the Spanish Ministry of Science and Technology (MCYT), the US National Science Foundation and the US Department of Energy. The isotopically enriched 76Ge was supplied by the Institute for Nuclear Research (INR), Moscow, and the Institute for Theoretical and Experimental Physics (ITEP), Moscow.