text:
# Multi-soliton energy transport in anharmonic lattices

## Abstract

We demonstrate the existence of dynamically stable multihump solitary waves in polaron-type models describing interaction of envelope and lattice excitations. In comparison with the earlier theory of multihump optical solitons \[see Phys. Rev. Lett. 83, 296 (1999)\], our analysis reveals a novel physical mechanism for the formation of stable multihump solitary waves in nonintegrable multi-component nonlinear models.

Spatially localised solutions of multi-component nonlinear models, multi-component solitary waves, have received a great deal of attention in the last decade. In particular, recent studies in the nonlinear optics  and Bose-Einstein condensation  have shown that, under certain conditions and only in multi-component systems, the formation of dynamically stable localized states and soliton complexes is possible. Unlike their single-component (or scalar) counterparts, multi-component (or vector) solitons possess complex internal structure forming a kind of “soliton molecules”, which makes them attractive, both from the fundamental and applied point of view, as composite and reconfigurable carriers for a transport of spatially localized energy.

Recent discovery of stable multi-component spatial solitons in optics  shed a light on the general physical mechanisms of the formation and stability of multi-component localised states. Such states are often called multihump solitons due to multiple maxima displayed in their intensity profile. Usually, multihump solitary waves appear via bifurcations of scalar solitons when a primary soliton plays a role of an effective waveguide (“potential well”) that traps higher-order guided modes excited in a complimentary field . On the other hand, the multihump solitons can be formed as multi-soliton bound states, when two or more different vector solitons are “glued” together due to balanced interaction between the soliton constituents .

Soliton bifurcations and binding enable the existence of multi-component localized states in many nonlinear models; these include bound states of dark solitons  and incoherent solitons  in optics, and multihump plasma waves . Importantly, multihump solitary waves are also found in higher dimensions .

The experimental and theoretical results on optical solitons  challenge the conventional view on multi-component solitary waves in other fields of nonlinear physics. The main question we wish to address here is: Can stable multihump solitons exist in other important models of nonlinear physics? This is a crucial issue because, so far, the stable multihump solitons have been positively identified only in the nonlinear optical model of Refs.  that is known to possess additional symmetries, which might be the reason for their unique stability.

In this Letter, we demonstrate the existence of dynamically stable multihump solitary states in a completely different (in both the physics and properties) but even more general model that describes the interaction of envelope and lattice excitations, a generalisation of the well-known polaron model. We reveal a novel physical mechanism for the formation of stable multihump solitary waves in nonintegrable multi-component nonlinear models.

Model. Let us consider the continuous model of the energy (or excess electron) transport in an anharmonic molecular chain , described by the system of coupled nonlinear Schrödinger (NLS) and Boussinesq equations:

$`i{\displaystyle \frac{\psi }{t}}`$ $`+`$ $`{\displaystyle \frac{^2\psi }{x^2}}+2w\psi =0,`$ (1)
$`m{\displaystyle \frac{^2w}{t^2}}`$ $`=`$ $`{\displaystyle \frac{^2w}{x^2}}+\mu {\displaystyle \frac{^4w}{x^4}}+\alpha {\displaystyle \frac{(w^2)}{x^2}}{\displaystyle \frac{(|\psi |^2)}{x^2}},`$ (2)
where $`t`$ and $`x`$ are the normalised time and spatial coordinate, correspondingly, $`\psi (x,t)`$ is the excitation wave function, and $`w(x,t)`$ is the chain strain. The system is characterized by three dimensionless parameters: the particle mass $`m`$, the anharmonicity of the chain $`\alpha `$, and the dispersion coefficient $`\mu `$.

Equation (1) appears in a number of other physical contexts including, for example, the interaction of nonlinear electron-plasma and ion-acoustic waves , coupled Langmuir and ion-acoustic plasma waves , interaction of optical and acoustic modes in diatomic lattices , particle theory models , etc.

System (1) is known to be integrable for $`\alpha \mu =6`$ . In this case, it possesses two types of single-soliton solutions: scalar ($`\psi =0`$) supersonic Boussinesq ($``$q) solitons and vector Davydov-Scott (DS) solitons  which can be both subsonic and supersonic. Because of the complete integrability for $`\alpha \mu =6`$, these solitons do not interact with each other. For $`\alpha \mu 6`$, the situation changes dramatically, and it has been recently shown  for the nearly-integrable case ($`\alpha \mu 6`$) that $``$q and DS solitons can form a bound state for $`\alpha \mu >6`$. From the other hand, it is also known that in a weakly anharmonic lattice two subsonic DS solitons can form a bisoliton . However, it remains a mystery what happens when the system (1) is far from its integrable limit and, especially, when the solitons are supersonic. In this Letter we examine the model (1) numerically for arbitrary values of $`\alpha \mu `$ and, employing the concept of soliton bifurcations, demonstrate the origin and exceptional robustness of multihump supersonic stationary solitary waves.

Bifurcation analysis. The stationary solutions of Eq. (1) can be found in the form of the traveling waves

$`\psi (x,t)={\displaystyle \frac{c}{4\sqrt{\mu }}}e^{ivx/2i\mathrm{\Omega }t}\varphi (z),w(x,t)={\displaystyle \frac{c}{4\mu }}u(z),`$ (3)
where $`z=\sqrt{c/4\mu }(xvt)`$, and the constant $`c=(mv^21)`$ is positive for a supersonic velocity, $`v>1/\sqrt{m}`$. Substituting Eq. (3) into Eq. (1), we derive a system of coupled ordinary differential equations

$`\varphi _{zz}\lambda \varphi +2u\varphi =0,`$ (4)
$`u_{zz}4u+gu^2\varphi ^2=0,`$ (5)
where $`g=\alpha \mu `$ is an effective anharmonicity parameter, and $`\lambda =\mu (v^24\mathrm{\Omega })/c`$ is a characteristic eigenvalue of the stationary localized solutions.

Equation (4) has two types of one-soliton solutions: a one-component $``$q soliton

$`\varphi _0(z)=0,u_0(z)=(6/g)\text{sech}^2(z),`$ (6)
which exists for arbitrary values of $`g`$, and a two-component DS soliton

$`\varphi _1=2\sqrt{\lambda (\lambda 1)}\text{sech}(\sqrt{\lambda }z),u_1=\lambda \text{sech}^2(\sqrt{\lambda }z),`$ (7)
which exists only in the integrable case $`g=6`$.

To understand what happens for $`g6`$, we consider the limit $`\varphi /u\epsilon 1`$ and apply a multi-scale asymptotic analysis. In the zeroth order in $`\epsilon `$, $`\varphi =0`$ and Eq. (4) reduces to a nonlinear equation for the component $`u(z)`$ only, with the supersonic $``$q soliton solution (6). In the first order in $`\epsilon `$, we obtain a linear eigenvalue problem for $`\varphi (z)`$ characterized by the effective potential $`u=u_0(z)`$. For a given value of $`g`$, the spectrum of the eigenvalue problem consists of $`N+1`$ discrete eigenvalues $`\lambda _n=\left(Nn\right)^2`$, where $`n=0,\mathrm{\hspace{0.33em}1},\mathrm{}N,`$ and $`N`$ is the integer part of $`(1/2)[\sqrt{1+(48/g)}1]`$. Each cut-off value $`\lambda _n`$ corresponds to a bifurcation point of the node-less scalar soliton $`u_0`$ where a two-component solution with a nonzero component $`\varphi `$ emerges. The latter has $`n`$ nodes and, near the bifurcation point, can be treated as a fundamental (or higher-order) bound mode of an effective potential created by the soliton $`u_0(z)`$. The emerging vector soliton can therefore be characterised by a “state vector” $`|0,n`$, according to the number of nodes in the corresponding components.

It is easy to see that for $`g>6`$ only bifurcations of the $`|0,0`$ state, which corresponds to the DS soliton (7), are possible. First bifurcation of the $`|0,1`$ state occurs for the completely integrable case $`g=6`$ at $`\lambda _1=0`$. In this case, the bifurcation pattern is identical to that of the completely integrable Manakov limit of two coupled NLS equations , namely, the $`|0,0`$ state appears at $`\lambda _0=1`$, and the $`|0,1`$ state appears at $`\lambda _1=0`$.

Weaker anharmonicity (smaller $`g`$) means larger number of possible bound states supported by the effective potential $`u(z)`$, and thus the increasing number of bifurcations. Indeed, the depth of the effective trapping potential is inversely proportional to $`g`$. The $`|0,0`$ state always exists, even for a shallow potential $`u(z)`$.

We now consider in detail the formation of multihump solitons in the cases of weak (subcritical, $`g<6`$) and strong (supercritical, $`g>6`$) anharmonicity, respectively.

Supercritical regime. In the absence of bifurcating higher-order solutions, the multihump solitons are formed only via binding of the $`|0,0`$ vector solitons. The physics of this mechanism is simple. The interaction forces between closely separated fundamental solitons are different for both the $`\varphi `$ and the $`u`$ components. Namely, while two in-phase $`u`$-solitons attract, the two in-phase $`\varphi `$-solitons repel. This allows for the existence of multihump nodeless modes of the field $`\varphi `$ trapped in the multi-well potential $`u`$. Each of such multihump solitons can be considered as a bound state of several $`|0,0`$ DS solitons, with in-phase humps in both components.

It is convenient to represent the solution families as branches on the bifurcation diagram $`Q`$ vs. $`\lambda `$, where $`QQ_u+Q_\varphi u^2𝑑z+\varphi ^2𝑑z`$ is the total soliton power. Typical bifurcation diagram for a supercritical case ($`g=7`$) is shown in Fig. 1. Solid line represents the bifurcating solution $`|0,0`$, and it can be seen on the close-up of the bifurcation region, that the branches representing two- and three-hump solutions start off at the same point $`\lambda =\lambda _0`$ but with the energies approximately equal to that of two or three $`u`$-solitons. Examples of such multihump solitons are shown in Fig. 2, and it is clear that this novel type of soliton bifurcations occurs from a countable set of infinitely separated single solitons. With increasing $`\lambda `$, separation between the humps decreases until all solitons of this type become single-humped (Fig. 2, right column).

Subcritical regime. In this case, bifurcations of the $`u`$-state do not lead to multihump solitons. That is, in sharp contrast to the coupled NLS equations describing vector solitons in nonlinear optics, none of the higher-order states $`|0,n`$ become multihumped in the system under consideration. Although the function $`\varphi `$ does have multiple maxima in its intensity profile, because of the non-self-consistent source for the $`u`$-component, it does not cause significant distortions in the shape of the effective potential, $`u(z)`$, and the total intensity $`I(z)=u^2+\varphi ^2`$ remains single-humped. Typical bifurcation diagram for the case $`g=1`$ is presented in Fig. 3. In this case $`N=2`$, and only bifurcations of the $`|0,1`$ (dot-dashed line) and $`|0,2`$ (solid line) solitons are shown. Corresponding modal profiles of the bifurcating solitons are presented in Fig. 3 (top row).

Similar to the supercritical regime, the multihump solitons can exist only as bound states of the bifurcating $`|0,1`$ or $`|0,2`$ solitons. They appear at the bifurcation points $`\lambda _n`$, and they have energies equal to a number of lower-order solitons “glued” together by the low-amplitude components. The number of nodes that the $`\varphi `$-component has in the composite soliton depends on the number of $`|0,n`$ solitons forming that bound state. Typical examples of such solutions are presented in Fig. 3 (bottom row).

From this analysis, we can conclude that, in this model, the multihump solitons appear as bound states of $`|0,0`$ solitons for any value of anharmonicity parameter $`g`$. In addition, multihump solitons of more sophisticated modal structure are also possible.

Dynamical stability. The second equation of the system (1) is the so-called “ill-posed” (or “bad”) $``$q equation . It possesses an intrinsic linear instability and therefore its reliable numerical solution for non-zero $`\mu `$ is unfeasible. This linear instability is not inherent in the original physical model, and it can be traced to neglecting higher-order spatial derivatives in Eq. (1). In the case of the energy transport in anharmonic molecular chains, Eq. (1) with $`\mu =1/12`$ originates from the following system of discrete equations

$`i{\displaystyle \frac{d}{dt}}\mathrm{\Psi }_n`$ $`=`$ $`\mathrm{\Delta }_2(\mathrm{\Psi }_n)(W_n+W_{n1})\mathrm{\Psi }_n,`$ (8)
$`m{\displaystyle \frac{d^2}{dt^2}}W_n`$ $`=`$ $`\mathrm{\Delta }_2(W_n)+\alpha \mathrm{\Delta }_2(W_n^2)`$ (9)
$`+`$ $`{\displaystyle \frac{1}{2}}(|\mathrm{\Psi }_{n+1}|^2+|\mathrm{\Psi }_n|^2|\mathrm{\Psi }_{n1}|^2|\mathrm{\Psi }_{n+2}|^2),`$ (10)
where $`\mathrm{\Delta }_2(X_n)X_{n+1}+X_{n1}2X_n`$. Discrete functions $`\mathrm{\Psi }_n`$ and $`W_n`$ define, in the continuous limit, the excitation wave function, $`\psi (x,t)`$, and the strain function of the lattice, $`w(x,t)`$.

Therefore, it would be justified to study the dynamics of the stationary solutions of Eqs. (4) numerically by employing the original discrete dynamical system (8). Besides, the argument can be reversed, and such a discrete system can be treated as a regularised numerical discretisation scheme for a general system of coupled NLS and ill-posed $``$q equations.

We investigate the dynamical stability of the multihump solitons for two distinct cases of a subcritical and supercritical anharmonicity discussed above. The condition of a unit norm for the envelope function $`\mathrm{\Psi }_n`$ ($`_n|\mathrm{\Psi }_n|^2=1`$) is satisfied in all cases, and $`v`$ is chosen close to the sound velocity ($`c=0.0244`$) to allow for a smooth discretisation.

In the supercritical regime ($`g>6`$), where multihump solitary waves can be formed through binding of several $`|0,0`$ states together, our numerical simulations indicate that such solitons are stable as long as the separation between the humps is sufficiently large. This property is just opposite to that observed for multihump optical solitary waves . An example of the stable dynamics of a three-hump soliton for $`g=7`$ is shown in Figs. 4(a,b). All solitons of the DS type, i.e. $`|0,0`$ soliton states, presented in Fig. 1 (by solid line) and Fig. 2 (bottom row), exhibit similar stable dynamics.

It is important that the same mechanism of the creation of multihump solitons applies to the subcritical regime. This means that the dynamically stable multi-soliton bound states described above exist also for $`g<6`$. As an example, propagation of a four-hump soliton at $`g=3.005`$ is demonstrated in Figs. 4(c,d). After initial adjusting of the soliton amplitudes (due to the discretisation), only small amplitude breathing occurs \[see Fig. 4(d)\], otherwise the soliton dynamics is stable. In contrast, all bifurcating higher-order solitons are dynamically unstable.

Our results on the robustness and stability of multi-soliton states call for a systematic revision of our understanding of the role of nonlinear localized modes in a number of physical phenomena related to the nonlinear transport in macromolecules  and even artificial nanoscale structures , where the coupling of two (or more) degrees of freedom occurs. How the soliton binding and existence of multi-soliton states modifies the nonlinear kinetics, nonequilibrium thermodynamics , and other properties of the system? These questions remain to be answered.

In conclusion, we have found robust two-component solitary waves in a polaron-type model of the energy transport in anharmonic lattices. We have revealed a novel physical mechanism for the formation of multihump solitons in a discrete anharmonic lattice and demonstrated their dynamical stability. Along with the recent studies on multihump optical solitons , these results call for re-examination of the role of multi-component solitary waves in other fields of nonlinear physics.

We thank A. V. Zolotaryuk and L. Cruzeiro-Hansson for helpful discussions. S. M. acknowledges support from the Australian Research Council. The work of E. O. and Yu. K. is supported by the Performance and Planning Fund of The Australian National University.


format: text

text:
# Wonderful blowups associated to group actions

## 1. Introduction

### 1.1.

Let $`X/`$ be a smooth variety, and for $`n>1`$ consider the product $`X^n`$. For any collection $`I=\{I_k\}`$ of disjoint subsets $`I_k\{1,\mathrm{},n\}`$ with each $`I_k`$ of order at least two, let $`\mathrm{\Delta }_I`$ be the polydiagonal consisting of all $`xX^n`$ with $`x_i=x_j`$ if and only if $`i,jI_k`$ for some $`k`$. The locally closed subvarieties $`\mathrm{\Delta }_I`$ provide a stratification of $`X^n`$ that is stable under the permutation action of the symmetric group $`S_n`$. In , Fulton and MacPherson constructed a smooth variety $`X[n]`$ and a birational morphism $`\pi :X[n]X^n`$ such that $`\pi `$ is an isomorphism over the open set $`\mathrm{\Delta }_{\mathrm{}}`$, and such that the stabilizers of the induced $`S_n`$ action on $`X[n]`$ are solvable. More recently, Ulyanov  constructed a space $`Xn`$ with a stratified birational morphism to $`X[n]`$ such that the stablilizers of the induced $`S_n`$ action are abelian. The variety $`Xn`$ is constructed by systematically blowing up the closures of the preimages of polydiagonals $`\mathrm{\Delta }_I`$ of increasing dimension.

Both $`X[n]`$ and $`Xn`$ are special cases of a general construction of MacPherson and Procesi . Suppose that $`Z`$ is an irreducible smooth variety and that $`S=\{S_i\}`$ is a set of locally closed subvarieties endowing $`Z`$ with a stratification. Suppose further that this stratification is *conical*; this means that every point $`zZ`$ has an analytic neighborhood isomorphic to a product of a trivial stratification of a disc and of a restriction to a disc of a $`^{}`$-invariant stratification of a vector space. A basic example is $`Z=X^n`$ and $`S=\{\mathrm{\Delta }_I\}`$. Then MacPherson and Procesi constructed a family of *wonderful* compactifications of the open stratum. These smooth compactifications map birationally onto $`Z`$ isomorphically away from strata of codimension bigger than one, and have the property that the complement of the open stratum is a normal crossing divisor. Moreover, for any $`S`$ there are canonical minimal and maximal wonderful compactifications. In the example above, these spaces are respectively $`X[n]`$ and $`Xn`$.

### 1.2.

In this paper we explain the connection

$`X[n]`$ $`S_n\text{ acts with solvable stabilizers}`$
$`Xn`$ $`S_n\text{ acts with abelian stabilizers}`$
by investigating wonderful blowups associated to certain stabilizer stratifications of a variety with group action. More precisely, let $`X`$ be a smooth variety and let $`G`$ be a finite group acting on $`X`$. We define two stratifications of $`X`$, the *stabilizer* and $`Y`$ stratifications, and show that they are conical. Our main results (Corollary 3.3, Theorem 4.3, and Corollary 5.3) are that $`G`$ acts on the maximal wonderful blowups associated to each of these stratifications with abelian stabilizers, and on the minimal wonderful blowup associated to the $`Y`$ stratification with solvable stabilizers. The main tool is the abelianization algorithm of Batyrev (, corrected in ).

### 1.3. Notations

We fix a finite group $`G`$ for the rest of this paper, and work in the category of smooth $`G`$-varieties over $``$, or smooth complex $`G`$-manifolds. We always assume that the action of $`G`$ is effective, i.e. that only the identity stabilizes all of any $`G`$-variety $`X`$. If $`X`$ is a $`G`$-variety, then for every point $`xX`$ its stabilizer is denoted by $`\mathrm{Stab}(x)`$. For every subgroup $`HG`$ the set of fixed points of $`H`$ is denoted by $`\mathrm{Fixed}(H)`$. The latter is always a disjoint union of smooth subvarieties of $`X`$.

## 2. Batyrev’s abelianization algorithm

### 2.1.

The goal of this section is to describe a version of Batyrev’s abelianization algorithm . Our version is a bit easier to formulate but generally results in more blowups.

Let $`X`$ be a $`G`$-variety as above. For every point $`xX`$ and any subgroup $`H`$ fixing $`x`$ ($`H`$ may be a proper subgroup of $`\mathrm{Stab}(x)`$) consider the decomposition of the tangent space $`TX_x`$ into irreducible $`H`$-modules

$$TX_x=\underset{i=1}{\overset{k}{}}V_i,$$

where $`dim(V_i)=1`$ for $`ik_1`$ and $`dim(V_i)>1`$ for $`i>k_1`$. Each one-dimensional representation $`V_i`$ gives a character $`H^{}`$. Denote by $`H_1`$ the common kernel of these characters, and let $`Y(x,H)`$ be the subvariety of $`X`$ that is the connected component of $`\mathrm{Fixed}(H_1)`$ passing through $`x`$. Observe that the tangent space to $`Y(x,H)`$ at $`x`$ is precisely $`_{i=1}^{k_1}V_i`$ \[2, Lemma 3\]. Notice that $`Y(x,H)=X`$ if and only if $`H`$ is abelian. Otherwise, the codimension of $`Y(x,H)`$ is at least $`2`$. We will consider the set of the proper subvarieties $`Y_jX`$ that are equal to $`Y(x,H)`$ for at least one pair $`(x,H)`$.

###### Theorem 2.2.

 In the above setup, denote by $`r(X)`$ the maximum number of different subvarieties $`Y_j`$ with non-trivial intersection. Denote by $`Z`$ the set of points $`xX`$ that are contained in $`r(X)`$ different $`Y_j`$. Then $`Z`$ is smooth and $`G`$-equivariant. We let $`X_1`$ be the blowup of $`X`$ along $`Z`$, and iterate the procedure: for $`i1`$ we compute $`r(X_i)`$ and $`Z_i`$, and define $`X_{i+1}`$ to be the blowup of $`X_i`$ along $`Z_i`$. We claim that this process terminates after at most $`r(X)`$ steps, and that the stabilizer of every point of the resulting variety is abelian.

###### Proof.

It is straightforward to see that the intersection of any number of the $`Y_j`$ is smooth. This implies that $`Z`$ is smooth, and it is clearly $`G`$-equivariant. It will now suffice to prove that $`r(X_1)<r(X)`$.

We claim that every $`Y_{j,1}`$ for $`X_1`$ is a proper preimage of some $`Y_j`$ for $`X`$ that is not contained in $`Z`$ and vice versa. Let $`Y_{j,1}=Y(x_1,H)`$ be a proper subvariety of $`X_1`$. Denote the image of $`x_1`$ under the blowdown map by $`x`$, and consider $`Y(x,H)`$. If $`xZ`$, there is nothing to prove. If $`xZ`$, notice that $`Z`$ is a subvariety of $`Y=Y(x,H)`$, so the tangent space $`TZ_x`$ is an $`H`$-submodule of

$$TY_x=\underset{i=1}{\overset{k_1}{}}V_i.$$

If the connected component of $`Z`$ passing through $`x`$ equals $`Y`$, then the fiber of the blowup over $`x`$, which equals $`(_{i>k_1}V_i)`$, has no $`H`$-fixed points. Therefore the tangent space to $`Z`$ is a proper submodule of $`TY_x`$, which we assume to equal $`_{i=1}^{k_2}V_i`$. The point $`x_1`$ then corresponds to a one-dimensional submodule of $`TY_x/TZ_x`$, which we will identify with $`V_{k_2+1}`$. The tangent space $`(TX_1)_{x_1}`$ is isomorphic as an $`H`$-module to

$$\underset{i=1}{\overset{k_2+1}{}}V_i\underset{i=k_2+2}{\overset{k}{}}(V_{k_2+1}V_i).$$

Therefore, the tangent space to $`Y_{j,1}`$ at $`x_1`$ is

$$\underset{i=1}{\overset{k_2+1}{}}V_i\underset{i=k_2+2}{\overset{k_1}{}}(V_{k_2+1}V_i),$$

which is the tangent space of the proper preimage of $`Y`$. The statement in the opposite direction is proved similarly.

To show that $`r(X_1)<r(X)`$, suppose that $`x_1`$ is contained in $`r(X)`$ different proper preimages of the $`Y_j`$. Then $`x_1`$ corresponds to a normal direction to $`Z`$ which is contained in all these $`Y_j`$’s, which leads to a contradiction. ∎

## 3. Batyrev’s abelianization as a maximal wonderful blowup

### 3.1.

The proof of Theorem 2.2 shows that once the $`Y_j`$ are defined for $`X`$, their exact nature is not important for the further blowups. Indeed, at every step one blows up the intersection set of the maximum number of the birational proper preimages of the $`Y_j`$. This allows us to identify the result of Batyrev’s algorithm with the *maximal wonderful blowup* of the MacPherson-Procesi family of blowups  associated to the stratification defined by the $`Y_j`$.

###### Definition 3.2.

The *$`Y`$ stratification* of $`X`$ is the stratification induced by the $`Y_j`$ as follows. For every subset of $`\{Y_j\}`$ one defines a stratum that consists of all points that lie in all $`Y_j`$ from the subset, but do not lie in any other $`Y_j`$. The empty strata are then ignored.

###### Corollary 3.3.

The $`Y`$ stratification of $`X`$ is conical. The maximal wonderful blowup coincides with the result of Batyrev’s algorithm above. In particular, it has only abelian stabilizers.

###### Proof.

First, it is easy to see that the $`Y`$ stratification is conical. Moreover, this can be said about any stratification induced by some connected components of fixed point sets of some subgroups of $`G`$. Indeed, for any $`xX`$ there exists a $`\mathrm{Stab}(x)`$-equivariant isomorphism of a neighborhood $`Ux`$ and a neighborhood of the origin in $`TX_x`$. Under this isomorphism all strata map to linear subspaces, so the stratification is conical.

Second, the maximal wonderful blowup is defined by successively blowing up proper preimages of all strata, starting with the strata of smallest dimension. This is equivalent to blowing up strata in any order that is compatible with the partial ordering on the strata, and Batyrev’s algorithm clearly provides that. ∎

## 4. The stabilizer stratification

### 4.1.

The $`Y`$ stratification is not local, in the sense that its construction does not commute with $`G`$-equivariant open embeddings. A more natural stratification of $`X`$ is induced by the set of all connected components of $`\mathrm{Fixed}(H)`$ for all subgroups $`HG`$. We will call it stabilizer stratification of $`X`$. It is also conical by the argument presented in the proof of Corollary 3.3.

###### Proposition 4.2.

Let $`X^Y`$ be the maximal wonderful blowup associated to the $`Y`$ stratification of $`X`$, and let $`X^{stab}`$ be the maximal wonderful blowup associated to the stabilizer stratification of $`X`$. Then there exists a natural $`G`$-equivariant map

$$X^{stab}X^Y.$$

###### Proof.

The statement is clearly local in $`X`$, so it could be assumed that $`X`$ is a neighborhood of the origin in a $`G`$-vector space $`V`$. Denote by $`S_Y`$ and $`S_{stab}`$ the sets of strata in the $`Y`$ and stabilizer stratifications respectively. According to \[3, Definition p.461, Proposition p.470\] the maximal blowups $`X^Y`$ and $`X^{stab}`$ can be described as closures of the images of an open subset of $`V`$ under the map to

$$V\times \underset{kS}{}(V/V_k).$$

Here $`V_k`$ is the tangent space to the stratum, and $`S`$ is either $`S_Y`$ or $`S_{stab}`$. The functoriality is then obvious. ∎

As a corollary we get the following theorem, which is the main result of this paper.

###### Theorem 4.3.

The stabilizer of every point of the maximal wonderful blowup associated to the stabilizer stratification is abelian.

###### Proof.

Combine Corollary 3.3 and Proposition 4.2. ∎

## 5. Solvable stabilizers

### 5.1.

It is generally much easier to ensure that a blowup of a $`G`$-variety has only solvable stabilizers, due to the following observation.

###### Proposition 5.2.

Let $`\pi :X_1X`$ be a $`G`$-equivariant birational morphism of $`G`$-varieties. Then all stabilizers of $`X_1`$ are solvable if and only if the image of the exceptional divisors of $`\pi `$ contains all points of $`X`$ with non-solvable stabilizers.

###### Proof.

Suppose $`X_1`$ has points with non-solvable stabilizers. Consider the point $`x_1X_1`$ with a non-solvable stabilizer $`H`$ that lies in the minimum number of exceptional divisors of $`\pi `$. If this minimum number is zero, then $`X`$ has a point with non-solvable stabilizer that lies outside the image of the exceptional divisors. Otherwise, let $`E`$ be an exceptional divisor of $`\pi `$ that passes through $`x`$. The tangent space $`TX_{x_1}`$ splits as an $`H`$-module into $`TE_{x_1}`$ and a one-dimensional module. This one-dimensional module corresponds to a character $`\chi :H^{}`$. The kernel of this character is again non-solvable; however its fixed point set contains points in the neighborhood of $`x_1`$ not lying in $`E`$. This contradicts the minimality of $`x_1`$ and proves the “if” part. The “only if” part is obvious. ∎

###### Corollary 5.3.

The minimal wonderful blowup associated to the $`Y`$ stratification has only solvable stabilizers.

###### Proof.

In view of Proposition 5.2, one needs to show that every point $`xX`$ with non-solvable stabilizer is contained in an irreducible stratum of codimension bigger than one. Since all $`Y_j`$ have codimension at least two, it remains to observe that each point $`x`$ with a non-solvable stabilizers is contained in $`Y=Y(x,\mathrm{Stab}(x))X`$. ∎

###### Question 5.4.

Is it true that all stabilizers of the minimal wonderful blowup associated with the stabilizer stratification are solvable?

This is really a question of whether there exists a non-solvable group $`G`$ acting on $`^d`$ with a fixed basis, such that for every $`gG`$ the space $`\mathrm{Ker}(g1)`$ is a coordinate subspace.

## 6. Examples

###### Example 6.1.

As in Section 1, let $`G`$ be the symmetric group $`S_n`$ acting on the product $`X^n`$ of $`n`$ copies of a variety $`X`$. Then the $`Y`$ stratification differs from the stabilizer stratification only in that it does not distinguish general points on the large diagonals (the conjugates of the preimages of the diagonals in $`X\times X`$) from the general points on $`X^n`$. As the result, the maximal wonderful blowup associated to the stabilizer stratification  is obtained from the maximal wonderful blowup for the $`Y`$ stratification by blowing up disjoint proper preimages of the large diagonals. Theorem 4.3 implies the result of Ulyanov  that the stabilizer of the maximal blowup is abelian. It is also easy to see that the discussion of Section 5 implies that the stabilizers of the minimal wonderful blowups are solvable, which was first observed in .

It is natural to ask whether for any finite group action of $`G`$ on $`X`$ there exists a *minimum abelianization*, defined to be the $`G`$-equivariant birational morphism $`\pi :X_1X`$ such that $`X_1`$ has abelian stabilizers and such that every other such map factors through $`\pi `$. Unfortunately, the next example shows that this is not the case.

###### Example 6.2.

Let $`X=^4`$ and let $`G=S_3`$ act on $`X`$ by linear transformations such that $`X`$ is a sum of two irreducible $`2`$-dimensional modules of $`G`$. Then a blowup at the origin or a blowup along any $`2`$-dimensional submodule of $`X`$ has only abelian stabilizers, but there is clearly no minimum abelianization.


text:
# Selfdual backgrounds in 𝑁=2 five-dimensional Chern-Simons Supergravity

## Abstract

We consider five-dimensional $`S(2,2|N)`$ Chern-Simons supergravity on $`M_4\times S_1`$. By fine-tuning the Kaluza-Klein reduction to make the four-dimensional cosmological constant equal zero, it is shown that selfdual curvatures on $`M_4`$ provide exact solutions to the equations of motion if $`N=2`$.

At first sight, the 3-dimensional Chern-Simons field theory defined by the equations of motion

$$g_{ab}F^b=0$$
(1)
seems to be completely trivial. However, it is now well-known that these equations and their generalizations including Wilson lines contain a large amount of information. Most notably, Witten  has shown that all known knots invariants can be described as expectation values of the corresponding Wilson lines. Another application of (1) is three-dimensional gravity. Choosing the Lie algebra to be $`SO(2,2)`$, Eqns. (1) are equivalent to Einstein equations in three dimensions. \[In (1), $`F^a`$ represent the Yang-Mills curvature 2-form for a Lie algebra $`G`$, and $`g_{ab}`$ is a $`G`$invariant tensor.\]

Equations (1) have a known generalization to five dimensions,

$$g_{abc}F^bF^c=0,$$
(2)
where $`g_{abc}`$ is a $`G`$invariant rank three tensor. Contrary to the three-dimensional case, these equations do not seem to have such a wide class of applications. In fact the structure of the space of solutions is largely unknown, and only some special aspects have been studied.

One of the aspects of (2) that has been studied in some detail is its interpretation as a five-dimensional gravitational theory . If the Lie algebra is taken to be $`SO(4,2)`$, this theory has a sensible interpretation as a gravitational theory with a Gauss-Bonnet coupling . Extending $`SO(4,2)`$ to $`SU(2,2|N)`$  we arrive at a supergravity theory in five dimensions . Contrary to the three- dimensional case, the relation between the $`d=5`$ Chern-Simons and standard supergravities is not known (see  for a recent discussion). The interesting feature of Chern- Simons supergravity is that the expected supergroups are incorporated from the very beginning, $`SU(2,2|N)`$ and $`OSp(32|1)`$ in five and eleven dimensions respectively, and the action is invariant by construction. The eleven dimensional case has been conjectured to be related to M-Theory in .

The main difference between (1) and (2) is that the latter has non-flat solutions. Moreover, propagating degrees of freedom can be shown to exist . A “natural” background to study (2) is the flat solution $`F^a=0`$. However, this phase is degenerate because small perturbations around $`F^a=0`$ are completely trivial. In fact, the linearized equations become an identity ***An interesting corollary of this fact is that the pp-wave solutions studied in , satisfying a linear wave equation, become trivial at the Chern-Simons point. In other words, when the coefficients are chosen as in the Chern-Simons action  the wave equation is lost.. Even though a generic Hamiltonian analysis does exhibit the presence of local excitations , no propagating solutions have been found so far.

In this paper we study the $`SU(2,2|N)`$ Chern-Simons theory on $`M_4\times \mathrm{}`$ and point out that selfdual curvatures on $`M_4`$ provide exact solutions to the equations of motion if $`N=2`$.

We start by recalling the construction of the $`SU(2,2|N)`$ Chern-Simons theory. The $`SU(2,2|N)`$ 1-form gauge field $`\mathrm{\Omega }=\mathrm{\Omega }_\mu dx^\mu `$ is a graded $`(N+4)\times (N+4)`$ complex matrix satisfying

$$(\mathrm{\Omega }\mathrm{\Gamma })^{}=\mathrm{\Omega }\mathrm{\Gamma },\text{STr}(\mathrm{\Omega })=0$$
(3)
where

$$\mathrm{\Gamma }=\left(\begin{array}{cc}\gamma _0& 0\\ 0& \text{I}_N\end{array}\right).$$
(4)
IN is the identity operator in $`N`$ dimensions and $`\gamma _0=\text{diag}(1,11,1)`$ We use the Dirac matrices in five dimensions $`\{\gamma _A,\gamma _B\}=2\eta _{AB}`$ with $`\eta _{AB}=\text{diag}(1,1,1,1,1)`$ and $`\gamma _5=\gamma _0\gamma _1\gamma _2\gamma _3`$.. The maximum number of supersymmetries is $`N=4`$ ($`𝒩=8`$) . Conditions (3) can be implemented by expanding $`\mathrm{\Omega }`$ in the form

$$\mathrm{\Omega }=\left(\begin{array}{cc}W& \psi \\ \overline{\psi }& 𝒜\end{array}\right)$$
(5)
where $`WU(2,2)`$ i.e., $`(W\gamma _0)^{}=W\gamma _0`$, $`\psi `$ is a Dirac spinor in five dimensions ($`\overline{\psi }=\psi ^{}\gamma _0`$) and $`𝒜^{}=𝒜`$ is a $`N\times N`$ matrix. Using Dirac matrices, $`W`$ and $`𝒜`$ can be further expanded as,

$`W`$ $`=`$ $`{\displaystyle \frac{i}{4}}A\text{I}_4+{\displaystyle \frac{i}{2}}e^A\gamma _A{\displaystyle \frac{1}{4}}w^{AB}\gamma _{AB},`$ (6)
$`𝒜`$ $`=`$ $`{\displaystyle \frac{i}{N}}A\text{I}_N+𝒜^{},`$ (7)
with $`𝒜^{}SU(N)`$, i.e., Tr$`𝒜^{}=0`$. It is direct to check that conditions (3) are indeed satisfied.

The 1-form $`e^A`$ is identified with the five-dimensional vielbein$`e^A`$ appears divided by the AdS radius which we set equal to one in order to simplify the notation., $`w^{AB}`$ is the spin connection and $`A`$ an Abelian $`U(1)`$ field. The Lagrangian has the usual five-dimensional Chern-Simons form. The explicit expression can be found in .

Our main concern will be the equations of motion which, after setting $`\psi =\overline{\psi }=𝒜^{}=0`$, take the form (see  for a detailed derivation of these equations),

$`ϵ_{ABCDE}\overline{R}^{AB}T^C`$ $`=`$ $`\overline{R}_{DE}F`$ (8)
$`ϵ_{ABCDE}\overline{R}^{AB}\overline{R}^{CD}`$ $`=`$ $`4T_EF`$ (9)
$`{\displaystyle \frac{1}{2}}\overline{R}^{AB}\overline{R}_{AB}T^AT_A`$ $`=`$ $`4(N^24^2)FF`$ (10)
where $`T^A=De^A`$ is the torsion tensor, $`F=dA`$, and

$$\overline{R}^{AB}=R^{AB}+e^Ae^B$$
(11)
with $`R^{AB}=dw^{AB}+w_C^Aw^{CB}`$. Note that all coefficients in (8,9,10) are fixed (up to trivial rescallings) by supersymmetry.

The anti-de Sitter (AdS) background has zero torsion, $`T^A=0`$, and constant spacetime curvature, $`R^{AB}+e^Ae^B=0`$. This configuration solves the equations if, and only if,

$$(N4)FF=0.$$
(12)
Anti-de Sitter space is then a solution only for the theory with the maximum number of supersymmetries $`N=4`$. The solution with $`FF=0`$ is discarded because it corresponds to a degenerate background . On the contrary, the phase for which $`F`$ has maximum rank, which is allowed by the $`N=4`$ equations of motion, carries the maximum number of degrees of freedom .

The $`N=4`$ theory is interesting in many other respects. For example under appropriated boundary conditions it can be shown that the group of asymptotic symmetries is the $`WZW_4`$ algebra . This yields a natural extension to five dimensions of the known relation between 3d Chern-Simons theory and the affine Kac-Moody algebra. In the context of the AdS/CFT correspondence, the $`N=4`$ theory seems to contain the correct ingredients; an AdS background and the maximum number of supercharges. Note also that $`SU(2,2|4)`$ represents both the super $`𝒩=8`$ AdS group in five dimensions and the superconformal $`𝒩=4`$ group in four dimensions (see  and references therein for a detailed study of this group).

Despite the interesting properties of the AdS solution and its associated supergravity theory, it is certainly not the only solution to the equations of motion. In order to further analyze the structure of Chern-Simons gravity we have studied other backgrounds with different number of supersymmetries. In particular, we consider the Chern- Simons theory on $`M_4\times \mathrm{}`$.

We shall perform a Kaluza-Klein reduction of equations (8,9,10). Our procedure will differ from standard Kaluza-Klein theory in that we do not impose the torsion condition. Instead, we shall dimensionally reduce the metric and connection separately. (The standard, second order, Kaluza-Klein reduction of 5d gravity with a Gauss-Bonnet term has been studied, for example, in  and leads to a non-linear electromagnetism.)

The reduction of spacetime indices is trivial. Let $`\mathrm{\Omega }_\mu dx^\mu =\mathrm{\Omega }_idx^i+\mathrm{\Omega }_5dx^5`$. Equations (8,9,10) are 4-forms in 5 dimensions with the structure $`ϵ^{\mu \nu \lambda \rho \sigma }_{\mu \nu \lambda \rho }=0`$. The projection to $`M_4`$ is achieved by setting the free index $`\sigma `$ equal to 5 obtaining $`ϵ^{ijkl}_{ijkl}=0`$. The remaining equations ($`\sigma =i`$) can be solved by the gauge condition $`\mathrm{\Omega }_5=0`$ plus the “zero mode” condition $`_5\mathrm{\Omega }_i=0`$ (see  for more details on the 4+1 decomposition of 5d Chern-Simons theory). In the following, we regard (8,9,10) as 4-forms in 4 dimensions.

The tangent indices are split in a similar way, $`A=\{a,5\}`$ with $`a=\{0,1,2,3\}`$. The vielbein has components $`\{e^a,e^5\}`$ and the spin connection $`\{w^{ab},w^{a5}\}`$. The pair $`e^a`$ and $`w^{ab}`$ are naturally the four dimensional gravitational variables. The 1-form $`e^5`$ will be related by the equations of motion to the gauge field $`A`$, just as in Kaluza-Klein theory. The 1-form $`w^{a5}`$, on the other hand, has a less clear role. Note that both $`e^a`$ and $`w^{a5}`$ transform as vectors under $`SO(3,1)`$ Lorentz rotations acting on the index $`a`$. We then led to consider the ansatz,

$$w^{a5}=\sigma _1e^a$$
(13)
where $`\sigma _1`$ is a constant with dimensions length-1. (Since we have set the AdS radius equal to one, $`\sigma _1`$ is actually dimensionless.)

The parameter $`\sigma _1`$ can be fine-tuned in order to neutralize the four-dimensional cosmological constant. In fact, the cosmological constant enter in the equations only through the combination (11). The four dimensional components of the curvature tensor are $`\overline{R}^{ab}=dw^{ab}+w_C^aw^{Cb}+e^ae^b`$. Using (13) we find,

$$\overline{R}^{ab}=R^{ab}+(\sigma _1^2+1)e^ae^b$$
(14)
where $`R^{ab}=dw^{ab}+w_c^aw^{cb}`$ is the intrinsic 4d curvature 2-form. Thus, if $`\sigma _1=\pm 1`$ there is no cosmological constant on $`M_4`$.

Two useful formulae that follows from (13) are

$$T^5=de^5,\sigma _1\overline{R}^{a5}=\overline{T}^a,$$
(15)
where $`\overline{T}^a=De^a+\sigma _1e^ae^5`$ denotes $`T^A`$ with $`A=a`$.

Inserting this ansatz into the equations of motion we find the reduced set of equations for $`e^a`$, $`w^{ab}`$ and $`A`$,

$`ϵ_{abcd}R^{ab}T^5`$ $`=`$ $`R_{cd}F`$ (16)
$`ϵ_{abcd}R^{ab}R^{cd}`$ $`=`$ $`4T_5F`$ (17)
$`R^{ab}R_{ab}2T^5T_5`$ $`=`$ $`8(N^24^2)FF`$ (18)
$`ϵ_{abcd}R^{ab}\overline{T}^c`$ $`=`$ $`\sigma _1\overline{T}_dF.`$ (19)
The configuration $`R^{ab}=F=T^5=0`$ is a solution to these equations reflecting the absence of a 4d cosmological constant. But we are interested in other backgrounds for which $`F`$ has the maximum rank.

By analogy with Kaluza-Klein theory and motivated by (16) we consider solutions on which the 1-form $`e^5`$ is proportional to the gauge field $`A`$. In this case, the curvature $`R^{ab}`$ will be selfdual. Real selfdual tensors exist only on Euclidean space and thus it would be more natural to analytically continue the above equations to Euclidean signature. However, defining super AdS in Euclidean space is not easy. Since we are interested in the full supergravity theory, in particular the existence of Killing spinors, we will work with Minkowskian signature and consider complex solutions to the equations of motion. We stress that from the purely bosonic point of view most imaginary units ‘$`i`$’ appearing below can be avoided in the Euclidean sector.

A solution to (16) is given by

$$T^5=\frac{i\sigma _2}{2}F,\frac{i\sigma _2}{2}ϵ_{abcd}R^{ab}=R_{cd}$$
(20)
where $`\sigma _2`$ takes the values $`\pm 1`$, and determines whether $`R^{ab}`$ is selfdual or antiselfdual.

Replacing (20) into (17) we find the condition on the gravitational and gauge anomalies,

$$R^{ab}R_{ab}=FF.$$
(21)
Next, we replace all together in (18), and discarding a non-zero factor of $`N`$, we find

$$(N2)FF=0,$$
(22)
which becomes an identity at $`N=2`$. Certainly this condition on the number of extended supersymmetries for the selfdual solutions to exist deserves a deeper explanation. It is interesting and gratifying that $`N`$ is an integer and lies in the range $`1N4`$. An independent derivation for this condition will arise below in the analysis of Killing spinors.

Finally we have to address Eq. (19). Let $`S_a=ϵ_{abcd}R^{ab}e^c\sigma _1Fe_d`$. In terms of $`S^a`$, Eq. (19) reads $`DS^a(i\sigma _1\sigma _2/2)AS^a=0`$. Since the first term in $`S_a`$ is the Einstein tensor, this equation can be regarded as a “generalized” conservation equation for $`S_a`$. Perhaps the most interesting solutions are the “vacuum” ones with $`S_a=0`$ which implies,

$$ϵ_{abcd}R^{ab}e^c=\sigma _1Fe_d.$$
(23)
Sumarising, the four-dimensional fields are $`e^a,w^{ab}`$ and $`A`$. The spin connection is constrained to be selfdual while $`A`$ is constrained by (21). The tetrad $`e^a`$ is determined by (19) or, in the vacuum case, by (23).

A deep analysis of these equations is beyond the scope of this work, we shall come back to this point elsewhere . Here we shall exhibit some solutions in order to prove that they are not self-contradictory.

First, let us point out that the “unbroken” phase $`e^a=0`$ solve the above equations. Let $`w^{ab}`$ be any selfdual connection, and $`A`$ a gauge field such that (21) is satisfied. The only remaining equation is (23) which is trivially solved by $`e^a=0`$.

It is natural to ask whether the Eguchi-Hanson , and other selfdual metrics provide solutions to the above equations. This turns out to be true although they need to be generalized. The reason is that the zero torsion condition $`De^a=0`$ imply $`FA=0`$ and then $`FF=0`$ (this follows directly from (19)). Since for the Eguchi-Hanson spacetime $`R^{ab}R_{ab}0`$, we conclude that this background will not solve (21).

There is, however, a direct modification to the Eguchi-Hanson metric which satisfies the full set of equations. Here we shall only present the general argument, for more details see . Let $`e_{\text{EH}}^a`$ the Eguchi-Hanson tetrad, and $`w^{ab}`$ its associated selfdual spin connection computed via $`De_{\text{EH}}^a=0`$. Let $`A`$ be a gauge field such that (21) is satisfied. A solution to the remaining equation (23) can be found by perturbing the tetrad as $`e^a=e_{\text{EH}}^a+f^a`$. Replacing in (23) we obtain,

$$(ϵ_{bcd}^aR^{cd}+\sigma _1\delta _b^aF)f^b=\sigma _1Fe_{\text{EH}}^a.$$
(24)
This is an algebraic equation of the form $`U_{ab}^{\mu \nu }f_\nu ^b=V_a^\mu `$. It can be checked  that det $`U0`$ and thus the perturbation $`f`$ is fully determined. Note that the torsion is not zero because $`T^a=De^a=D(e_{\text{EH}}^a+f^a)=Df^a0`$. For more details on this construction see .

A more explicit, although less interesting, family of solutions can be found as follows. Let $`w^{ab}=\mathrm{\Lambda }^{ab}B`$ where $`B`$ is an Abelian 1- form, and $`\mathrm{\Lambda }^{ab}`$ a constant selfdual matrix (note the position of the indices)

$$\mathrm{\Lambda }_b^a=\frac{1}{2}\left(\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& i\\ 0& 0& i& 0\end{array}\right).$$
(25)
The curvature becomes $`R^{ab}=\mathrm{\Lambda }^{ab}dB`$ and it is selfdual because $`\mathrm{\Lambda }^{ab}`$ is. We have chosen in this example $`\sigma _2=1`$.

Note that since the curvature is selfdual, Eq. (19) can be arranged as a form-valued eigenvalue equation,

$$R_b^a\overline{T}^b=\frac{i\sigma _1\sigma _2}{2}F\overline{T}^a$$
(26)
A solution to (26) is given by,

$$dB=iF,\mathrm{\Lambda }_b^a\overline{T}^b=\frac{\sigma _1}{2}\overline{T}^a$$
(27)
The first relation imply $`B=iA`$, while the second

$$d(e^0\sigma _1e^1)=0,d(ie^2+\sigma _1e^3)=0.$$
(28)
The local solution to (28) depends on two functions $`f`$ and $`g`$, $`e^0\sigma _1e^1=df`$ and $`ie^2+\sigma _1e^3=dg`$. Note, in particular, that the flat space configuration $`e^a=dx^a`$ is a solution. Finally, we note that (21) is automatically satisfied,

$$R^{ab}R_{ab}=R_b^aR_a^b=\text{Tr}(\mathrm{\Lambda }^2)dBdB=FF.$$
(29)
Let us now analyze the spinor Killing equation on the selfdual background. The supersymmetry transformations in this theory are of the form $`\delta \mathrm{\Omega }=_\mathrm{\Omega }\mathrm{\Sigma }`$ where $`\mathrm{\Sigma }SU(2,2|N)`$ and $`\mathrm{\Omega }`$ is given in (5). These transformations form a closed algebra and differ from the standard $`d=5`$ supergravity formula. The Killing spinor equation is

$$\delta \psi _k=ϵ_k\frac{i}{N}Aϵ_k=0,k=1,2,\mathrm{},N$$
(30)
where $`=d+W`$ is the $`U(2,2)`$ covariant derivative ($`W`$ is given in (6)). The selfdual solution requires $`N=2`$. This condition can actually be re-obtained from the Killing equation.

Applying the derivative $``$ we find the integrability condition,

$$\left(\frac{i}{4}F+\frac{i}{2}T^A\gamma _A\frac{1}{4}\overline{R}^{AB}\gamma _{AB}\right)ϵ_k=\frac{i}{N}Fϵ_k$$
(31)
where the term in brackets is the $`U(2,2)`$ curvature. The term proportional to the torsion $`\overline{T}^a`$ introduces big complications in this equation because it depends on the potential $`A`$ and not only on the curvatures. Fortunately, thanks to (15) this term can be eliminated by imposing the chirality condition,

$$\gamma _5ϵ_k=i\sigma _1ϵ_k.$$
(32)
The integrability condition (31) has now two sectors depending on the relative signs of $`\sigma _1`$ and $`\sigma _2`$ ($`N=2`$)

$$\frac{1}{2}R^{ab}\gamma _{ab}ϵ_k=\{\begin{array}{cc}iFϵ_k& \sigma _1=\sigma _2,\hfill \\ 0& \sigma _1=\sigma _2,\hfill \end{array}$$
(33)
On the exact solution (25,27) with $`\sigma _2=1`$ this equations becomes

$$\frac{1}{2}(\gamma _0\gamma _1+i\gamma _2\gamma _3)ϵ_k=\{\begin{array}{cc}ϵ_k& \sigma _1=1,\hfill \\ 0& \sigma _1=1,\hfill \end{array}$$
(34)
The eigenvalues of $`(\gamma _0\gamma _1+i\gamma _2\gamma _3)/2`$ are $`\{1,1,0,0\}`$. Therefore, the sector $`\sigma _1=1`$ has two solutions, and it can be proved that both of them satisfies the chirality condition (32). This sector preserves 1/2 supersymmetries. On the other hand, the sector $`\sigma _1=1`$ has only one solution which satisfies (32) as well. This sector preserves 1/4 supersymmetries.

We can say something else about the integrability condition (33) which is independent of particular examples. Consider first the case $`\sigma _1=\sigma _2`$. We multiply Eq. (33) by $`F`$. Using the identity $`\frac{1}{2}\{\gamma _{ab},\gamma _{cd}\}=ϵ_{abcd}\gamma _5+\delta _{[ab][cd]}`$ plus the selfduality of $`R^{ab}`$ and chirality of $`ϵ_k`$ one finds,

$$R^{ab}R_{ab}ϵ_k=FFϵ_k$$
(35)
which becomes an identity on the selfdual background satisfying (21). Regarding the sector $`\sigma _1=\sigma _2`$, we note that any selfdual linear combination of Dirac matrices $`R^{ab}\gamma _{ab}`$ has two zero eigenvalues. Thus (33) always have non-zero solutions. This analysis indicates that all selfdual solutions preserve some supersymmetry.

Note that the right hand side of (31) contains a factor $`1/N`$. This factor is present in all the following equations (33), (34) and (35). Non-zero Killing spinors exist only for $`1/N=1/2`$, as expected from (22). In this sense, the Killing spinor equation provides an independent derivation for the condition $`N=2`$.

To conclude, we have explored five-dimensional $`SU(2,2|2)`$ Chern-Simons supergravity compactified to four dimensions, and show that exact, selfdual, solutions exist. These solutions preserve some supersymmetry and we have exhibited some explicit examples. More general solutions will be presented elsewhere .

The author would like to thank A. Chamblin, G. Gibbons, M. Henneaux, A. Ritz, R. Troncoso and J. Zanelli for discussions and useful correspondence. This work was supported in part by FONDECYT (Chile) grant # 1000744 and DICYT (USACH) # 04-0031BL.


text:
# Molecular Gas in Elliptical Galaxies

## 1. Introduction

Not long ago it became clear that many elliptical galaxies do have small amounts of cold interstellar gas and dust (e.g. Knapp et al. 1989; Goudfrooij et al. 1994). The gas and dust have significance far beyond their small masses: their distribution and kinematics provide vital clues to the evolution of ellipticals (van Gorkom & Schiminovich 1997; Sadler et al. 2000). To date, most work on the distribution and kinematics of neutral gas in ellipticals has focused on the atomic gas.

The present contribution describes a survey to make maps of the CO emission in ten elliptical galaxies. The main selection criterion for this sample (after the classification issues, which are admittedly difficult) is the galaxies’ single dish CO line strengths. The resolution of the CO data is about 6<sup>′′</sup>, which is 20 pc (!) for the nearest galaxies and 1.5 kpc for the most distant ones.

The goals of the project are to understand (1) the origin of the molecular gas in ellipticals and its prospects for long-term stability; and (2) how much star formation activity is taking place in ellipticals. The molecular gas distribution and kinematics are being compared with the atomic gas, stellar kinematics, and radio continuum and H$`\alpha `$ emission. The sample galaxies are found in a variety of environments— some are isolated, some in loose groups, and some in the center of Virgo cluster. They also have a wide range of luminosities; some have active nuclei and some show visible signs of recent interactions.

Seven of the ten sample galaxies have already been observed with the BIMA millimeter interferometer, and five are clearly detected. Figure 1 shows the stellar distributions and CO maps for two of them.

## 2. Some Initial Results

Except for NGC 185 and NGC 205, the smallest members of the sample, the molecular gas is distributed in a smooth, regular disk in solid body rotation. There is a strong asymmetric warp in the CO of NGC 3656, but otherwise the disks are not disturbed. In two cases the galaxies’ stars are known to have large velocity dispersions and small rotational velocities; the CO rotation is consistent with that of the stars. Usually the HI kinematics are found to be inconsistent with those of the stars, which is taken as evidence for an external gas origin. Perhaps the situation is more complicated for the molecular gas.

A remarkable result comes from comparing the molecular gas distribution in NGC 5666 to the radio continuum emission in that galaxy (Wrobel & Heeschen 1988). The radio continuum emission is extended; it has exactly the same extent as the resolved CO emission. Thus, the radio continuum and FIR emission can be attributed to star formation (3 M per year) rather than to an active galactic nucleus or to an old stellar population. The other galaxies in the sample will help to show whether this situation is common or uncommon.

Comparisons of the molecular gas properties and the star formation activity in the different galaxies (isolated vs. in clusters, giant vs. dwarf, and so on) will allow us to probe elliptical galaxy evolution as a function of galaxy properties and/or environment.

## References

Goudfrooij, P., Hansen, L., Jorgensen, H. E. and Norgaard-Nielsen, H. U. 1994, A&AS, 105, 341

Knapp, G.R., Guhathakurta, P., Kim, D.-W., & Jura, M. 1989, ApJS, 70, 329

Sadler, E. M., Oosterloo, T. A., Morganti, R., & Karakas, A. 2000, AJ 119, 1180

van Gorkom, J., & Schiminovich, D. 1997, ASP Conf. Ser. 116: The Nature of Elliptical Galaxies; 2nd Stromlo Symposium, 310

Wrobel, J.M. & Heeschen, D.S. 1988 ApJ 355, 677

Young, L.M., & Lo, K.Y. 1996, ApJ 464, L59


text:
# Cosmic Ray Air Shower Characteristics in the Framework of the Parton-Based Gribov-Regge Model NEXUS

## 1 Introduction

Although cosmic rays have been studied for many decades, there remain still many open questions, in particular concerning the high energy cosmic rays above $`10^{14}`$ eV. One does neither know their composition nor the sources and acceleration mechanisms, partly because of the fact that at these high energies direct measurements are impossible due to the weak flux. But since cosmic ray particles initiate cascades of secondaries in the atmosphere, the so-called extensive air showers, one may reconstruct cosmic ray primaries by measuring shower characteristics. This reconstruction requires, however, reliable model predictions for the simulation of EAS initiated by either protons or nuclei from helium to iron. The problem is that the energy of the cosmic rays may exceed by far the energy range accessible by modern colliders, where at most equivalent fixed target energies of roughly $`10^{15}`$ eV can be reached. The projects of new generation EAS arrays are aimed even to the energy region $`10^{20}÷10^{21}`$ eV, and so the gap between the existing energy limit of collider data and the demands of cosmic ray experiments is considerable. Moreover, the real reliable data limits are essentially less than above mentioned, because at present time collider experiments do not register particles going into the extreme forward direction and a number of other drawbacks may be listed. Therefore there exists a real need of “reasonable” models, implementing the correct physics, in order to be able to make extrapolations towards extremely high energies.

Concerning models one has to distinguish between EAS models and hadronic interaction models. The latter ones like venus , qgsjet , and sybill  are modeling hadron-hadron, hadron-nucleus and nucleus-nucleus collisions at high energies, being more or less sophisticated concerning the theoretical input, and relying in any case strongly on data from accelerator experiments. The EAS models like corsika  are actually simulating the full cascade of secondaries, using one of the above-mentioned hadronic models for the hadronic interactions and treating the well known electro-magnetic part of the shower. It turns out that the model predictions of EAS simulations depend substantially on the choice of the hadronic interaction model. In corsika, the average electron number in EAS at primary energy $`10^{15}`$ eV varies from $`1.1110^5`$ to $`1.6210^5`$ (at sea level) depending on the hadronic interaction model . So the right choice of the model and its parameters is extremely important.

Recently a new hadronic interaction model neXus  has been proposed. It is characterized by a consistent treatment for calculating cross sections and particle production, considering energy conservation strictly in both cases (which is not the case in all above-mentioned models!). In addition, one introduces hard processes in a natural way, avoiding any unphysical dependence on artificial cutoff parameters. A single set of parameters is sufficient to fit many basic spectra in proton-proton and lepton-nucleon scattering, as well as in electron-positron annihilation. Briefly: concerning theoretical consistency, neXus is considerably superior to the presently used approaches, and allows a much safer extrapolation to very high energies.

This new approach cures some of the main deficiencies of two of the standard procedures currently used: the Gribov-Regge theory and the eikonalized parton model, the theoretical basis of the above-mentioned interaction models. There, cross section calculations and particle production cannot be treated in a consistent way using a common formalism. In particular, energy conservation is taken care of in case of particle production, but not concerning cross section calculations. In addition, hard contributions depend crucially on some cutoff and diverge for the cutoff being zero.

Having a reliable hadronic interaction model, one may now proceed to do air shower calculations. It may appear that an ideal solution from the user’s point of view is to use direct Monte Carlo technique, where the cascade is traced from the initial energy to the threshold one and the threshold energy corresponds to the minimum energy registered by the array in question. But such an approach takes unreasonably much computer time and sometimes gives no possibility to analyze experimental data in the appropriate way. Already at $`10^{14}`$ eV to $`310^{16}`$ eV, the experimental statistics from the kascade experiment  is ten times bigger than the corresponding number of simulated events. It should be mentioned that at higher energies (above $`10^{17}`$ eV) it is hardly possible to use direct Monte Carlo calculations. Of course, some modifications, like the thinning method , are possible and average values may be computed.

In this paper, we follow an alternative approach, where one treats the air shower development in terms of cascade equations. Here, the cascade evolution is characterized by differential energy spectra of hadrons, which one obtains by solving a system of integro-differential equations. Crucial input for these equations is the inclusive spectra of hadrons produced in hadron-air collisions. These spectra are obtained by performing neXus simulations and finally parameterizing the results. Having solved the cascade equations, we finally analyse the results of the casa-blanka group as a first application of our approach.

## 2 The NEXUS Model

The most sophisticated approach to high energy hadronic interactions is the so-called Gribov-Regge theory . This is an effective field theory, which allows multiple interactions to happen “in parallel”, with phenomenological object called “Pomerons” representing elementary interactions . Using the general rules of field theory, one may express cross sections in terms of a couple of parameters characterizing the Pomeron. Interference terms are crucial, they assure the unitarity of the theory.

A big disadvantage is the fact that cross sections and particle production are not calculated consistently: the fact that energy needs to be shared between many Pomerons in case of multiple scattering is well taken into account when considering particle production (in particular in Monte Carlo applications), but not for cross sections .

Another problem is the fact that at high energies, one also needs a consistent approach to include both soft and hard processes. The latter ones are usually treated in the framework of the parton model, which only allows to calculate inclusive cross sections.

We recently presented a completely new approach  for hadronic interactions and the initial stage of nuclear collisions, which is able to solve several of the above-mentioned problems. We provide a rigorous treatment of the multiple scattering aspect, such that questions of energy conservation are clearly determined by the rules of field theory, both for cross section and particle production calculations. In both (!) cases, energy is properly shared between the different interactions happening in parallel, see fig. 1 for hadron-nucleus collisions. This is the most important and new aspect of our approach, which we consider a first necessary step to construct a consistent model for high energy nuclear scattering.

The elementary diagram, shown as the thick lines in the above diagrams, is the sum of the usual soft Pomeron and the so-called semi-hard Pomeron, where the latter one may be obtained from perturbative QCD calculations (parton ladders). To some extent, our approach provides a link between the Gribov-Regge approach and the parton model, we call it “Parton-based Gribov-Regge Theory”.

## 3 System of Hadronic Cascade Equations

EAS are produced as a result of the hadronic cascade development in the atmosphere. We characterize hadronic cascades by the differential spectra $`h_n(E,X)`$ of hadrons of type $`n`$ with energy $`E`$ at an atmospheric depth $`X`$, the latter one being the integral over the atmospheric density $`\rho `$ along a straight line trajectory (not necessary radial) from some point $`P`$ to infinity,

$$X=_P^{\mathrm{}}\rho (x)𝑑x,$$
(1)
measured usually in $`\mathrm{g}/\mathrm{cm}^2`$. The decrease of the average hadron numbers due to collisions with air nuclei is given as

$$\frac{dh_n}{dX}=\frac{h_n}{\lambda _n},$$
(2)
where the mean inelastic free path $`\lambda _n`$ (in units of mass/area) can be expressed via the average hadron-air cross-section $`\sigma _{\mathrm{inel}}^{(n)}`$ and the average mass of air molecules $`m_{air}`$ :

$$\lambda _n=\frac{m_{\mathrm{air}}}{\sigma _{\mathrm{inel}}^{(n)}}$$
(3)
The second process to be considered is particle decay. The decay rate in the particle c.m. system is $`dh_n/d\tau =h_n/\tau _0`$, with $`\tau _0`$ being the particle life time. For a relativistic particle we find

$$\frac{dh_n}{dX}=\frac{B_n}{EX}h_n,$$
(4)
with the decay constant in energy units

$$B_n=\frac{m_n}{\alpha c\tau _0},$$
(5)
where $`m_n`$ is the hadron mass and $`c`$ the velocity of light. The ratio $`\alpha =\rho /X`$ depends only weakly on $`X`$ and is for the following taken to be constant, which implies constant $`B_n`$. If we take the simple exponential barometric formula for the density $`\rho `$, we get

$$\alpha =\frac{\rho _0g\mathrm{cos}\theta }{P_0},$$
(6)
where $`\rho _0`$ and $`P_0`$ are density and pressure at (for example) sea level, $`g`$ the gravitational acceleration of the earth, and $`\theta `$ the zenith angle of the shower trajectory. In this paper, we only consider the case $`\theta =0`$.

Based on the above discussion, a system of integro-differential equations for the differential energy spectra $`h_n`$ of hadrons may be presented as

$`{\displaystyle \frac{h_n(E,X)}{X}}`$ $`=`$ $`h_n(E,X)\left[{\displaystyle \frac{1}{\lambda _n(E)}}+{\displaystyle \frac{B_n}{EX}}\right]`$
$`+{\displaystyle \underset{m}{}}{\displaystyle _E^{E_{\mathrm{max}}}}h_m(E^{},X)\left[{\displaystyle \frac{W_{mn}(E^{},E)}{\lambda _m(E^{})}}+{\displaystyle \frac{B_mD_{mn}(E^{},E)}{E^{}X}}\right]𝑑E^{},`$
The quantities $`W_{mn}(E^{},E)`$ and $`D_{mn}(E^{},E)`$ are the inclusive spectra of secondaries of type $`n`$ and energy $`E`$ which are produced in interactions ($`W`$) or decays ($`D`$) of primaries of type $`m`$ and energy $`E^{}`$. The energy $`E_{\mathrm{max}}`$ is the maximum energy considered. If the type of primary hadron is $`n_o`$, its energy $`E_0`$ and the cascade originates at $`X_o`$, then one should add

$$h_n(E,X=X_o)=\delta _{nn_o}\delta (EE_0)$$
(8)
as the boundary condition. The more detailed consideration of the problem may be found elsewhere (see ).

It is reasonable to incorporate in the system nucleons (and anti-nucleons), charged pions ($`B_\pi =114`$ GeV), charged kaons ($`B_K=852`$ GeV), and neutral kaons ($`B_{K_L^o}=205`$ GeV). As $`B_{K_S^o}=1.1910^5`$ GeV, there is no sense to account for neutral kaons $`K_o^S`$ at energies $`B_{K_S^o}`$. But these particles should be included if their energy exceeds $`0.01B_{K_S^o}`$. The values of decay constants are given at the height $`11`$ km.

The computational technique to solve the system (3) is based on the same principle as the traditional approaches , but some improvements are introduced, which enable one to avoid too small steps when integrating over the depth . One discretizes the energy as

$$E_i=E_{\mathrm{min}}c^i$$
(9)
with $`E_{\mathrm{min}}1`$ GeV and $`c`$ such that the number of points per order of magnitude is $`10÷20`$. Replacing the integral in the right-hand side of (3) by the corresponding sum, one may write

$$\frac{h_{ni}(X)}{X}=h_{ni}(X)\left[\frac{1}{\lambda _{ni}}+\frac{B_n}{E_iX}\right]+\underset{m}{}\underset{j_{min}}{\overset{j_{\mathrm{max}}}{}}h_{mj}(X)\left[\frac{W_{mn}^{ji}}{\lambda _{mj}}+\frac{B_mD_{mn}^{ji}}{E_jX}\right],$$
(10)
with $`h_{ni}(X)=`$$`h_n(E_i,X)`$, $`\lambda _{ni}=\lambda _n(E_i)`$, and

$`W_{mn}^{ji}`$ $`=`$ $`{\displaystyle _{E_{\mathrm{min}}(i)}^{E_{\mathrm{max}}(i)}}{\displaystyle \frac{E}{E_i}}W_{mn}(E_j,E)𝑑E`$ (11)
$`D_{mn}^{ji}`$ $`=`$ $`{\displaystyle _{E_{\mathrm{min}}(i)}^{E_{\mathrm{max}}(i)}}{\displaystyle \frac{E}{E_i}}D_{mn}(E_j,E)𝑑E`$ (12)
with $`E_{\mathrm{min}}(i)=E_i/\sqrt{c}`$ and $`E_{\mathrm{max}}(i)=E_i\sqrt{c}`$ ( if $`i=j`$ then $`E_{max}(i)=E_j`$ ). The factor $`E/E_i`$ has been added in the integral to ensure exact energy conservation. The spectra $`W_{mn}(E_j,E)`$ and $`D_{mn}(E_j,E)`$ must be provided in order to calculate the above integrals. The standard recipe for calculating $`D_{mn}(E_j,E)`$ may be found elsewhere ( see  or  ) and implies no difficulties. The calculation of $`W_{mn}(E_j,E)`$ must be based on neXus simulations, as discussed below in detail. Since the $`E_j`$ dependence is very smooth, there is no point to calculate the spectra for all energies $`E_j`$. One rather chooses some reference energies $`E_j^{}`$ to calculate the spectra based on neXus simulations, and then interpolates to obtain the spectra for the other energies. This method has proven to be superior concerning the computational time when compared with the direct spectra calculations over all energies. For the calculations in this paper we choose two reference points per order of magnitude, starting from $`10^{11}`$ eV: $`E_j^{}=10^{11}`$ eV, $`10^{11.5}`$ eV, $`10^{12}`$ eV etc. up to $`10^{17}`$ eV which is at present the maximum energy attainable in the neXus model. As the neXus model is not valid at energies below $`10^{11}`$ eV, the data obtained with other codes must be borrowed. In this work we employ results obtained in ,  which are close to predictions of the gheisha code  used in corsika.

The solution of the homogeneous equation

$$\frac{h_{ni}(X)}{X}=h_{ni}(X)\left(\frac{1}{\lambda _{ni}}+\frac{B_n}{E_iX}\right)$$
(13)
has the form

$$h_{ni}(X_o+\mathrm{\Delta }X)=h_{ni}(X_o)\mathrm{exp}\left\{\frac{\mathrm{\Delta }X}{\lambda _{ni}}\right\}\left(\frac{X_o}{X_o+\mathrm{\Delta }X}\right)^{B_n/E_i},$$

if for simplicity we neglect the weak dependence of $`B_n`$ on $`X`$. The solution of the full equation (10) may be written as

$`h_{ni}(X_o+\mathrm{\Delta }X)`$ $`=`$ $`h_{ni}(X_o)\mathrm{exp}\left\{{\displaystyle \frac{\mathrm{\Delta }X}{\lambda _{ni}}}\right\}\left({\displaystyle \frac{X_o}{X_o+\mathrm{\Delta }X}}\right)^{B_n/E_i}`$
$`+{\displaystyle \underset{m}{}}{\displaystyle \underset{j_{min}}{\overset{j_{max}}{}}}{\displaystyle _{X_o}^{X_o+\mathrm{\Delta }X}}h_{mj}(X^{})\left[{\displaystyle \frac{W_{mn}^{ji}}{\lambda _{mj}}}+{\displaystyle \frac{B_mD_{mn}^{ji}}{E_jX^{}}}\right]`$
$`\mathrm{exp}\left\{{\displaystyle \frac{X_o+\mathrm{\Delta }XX^{}}{\lambda _{ni}}}\right\}\left({\displaystyle \frac{X^{}}{X_o+\mathrm{\Delta }X}}\right)^{B_n/E_i}dX^{}`$
which may be verified directly. The above formula allows to calculate $`h_{ni}`$ at depth $`X_o+\mathrm{\Delta }X`$, provided that $`h_{ni}(X)`$ at $`X=X_o`$ is known and all $`h_{mj}(X)`$ are also known for $`X_oXX_o+\mathrm{\Delta }X`$ and $`j>i`$. So, starting from $`h_{ni}(X_o)`$ one may sequentially find $`h_{ni}(X_o+\mathrm{\Delta }X)`$ and so on. Test calculations show that for $`E_i>B_n/3`$ it is quite sufficient to use Simpson’s formula. Here, one needs the values for $`h_{mj}(X_o+\mathrm{\Delta }X/2)`$, which are obtained via interpolation. Whereas for pions this works without problem, for kaons the requirement $`E_i>B_{K^\pm }/3`$ becomes more restrictive and, what is especially essential, errors for one component manifest themselves in other components. In order to retain accuracy without resorting to excessively small $`\mathrm{\Delta }X`$, the integration over $`X^{}`$ in eq. (3) is approximated as

$$\underset{X_o}{\overset{X_o+\mathrm{\Delta }X}{}}f(X^{})\left(\frac{X^{}}{X_o+\mathrm{\Delta }X}\right)^{B_\pi /E_i}𝑑X^{}=A_1f(X_o)+A_2f(X_o+\mathrm{\Delta }X/2)+A_3f(X_o+\mathrm{\Delta }X)$$
(15)
and coefficients $`A_1`$, $`A_2`$, $`A_3`$ are found from the condition that (15) is exact for a second order polynomial. The accuracy of $`1\%`$ may be achieved with $`\mathrm{\Delta }X5`$ g/cm2.

Other EAS characteristics (e.g. electron and muon numbers) are computed in a traditional way as corresponding functionals from functions $`h_n(E,X).`$ Usually one assumes that neutral pions decay immediately at the generation point and do not contribute to the development of the hadronic cascade. This assumption is quite adequate at energies below the corresponding decay constant which for neutral pions is about $`310^{19}`$eV. Moreover, as primary particles are nucleons there is an additional factor of $``$10 in our favour. So the number of neutral pions produced at depth $`X`$ may be obtained as

$$h_{\pi ^0}(E_i,X)=\underset{m}{}\underset{j=i}{\overset{j_{max}}{}}h_{mj}(X)\left[\frac{1}{\lambda _{mj}}W_{m\pi ^0}^{ji}+\frac{B_m}{E_jX}D_{m\pi ^0}^{ji}\right]$$
(16)
where we substitute $`n`$ by $`\pi ^0`$ in the second term of the right-hand side of eq. (10). The electron number $`N_e`$ at depth $`T`$ is given as

$$N_e(T)=\underset{i=1}{\overset{i_{max}}{}}_0^Th_{\pi ^0}(E_i,X_0)N_G(y,t)\frac{2}{1+s}𝑑X_0$$
(17)
where $`y=\mathrm{ln}(E_i/\beta )`$ is the logarithmic energy in units of the critical energy of electrons in air ($`\beta =81.10^6`$ eV) and $`t=(TX_0)/T_0`$ is the depth difference in radiation units ($`T_0=`$37.1 g/cm2). The factor $`2/(1+s)`$ accounts for energy sharing between two photons, with $`s`$ being the shower age parameter $`s=3t/(t+2y)`$. The function $`N_G`$ is referred to as Greisen’s formula ,

$$N_G=\frac{0.31}{\sqrt{y}}\mathrm{exp}\left\{t\left[1\frac{3}{2}\mathrm{ln}(s)\right]\right\},$$
(18)
which predicts the electron number at depth $`T`$ in the shower, produced by a primary photon with energy $`E_i`$ at depth $`X_0`$.

As a rule, experimental EAS arrays can detect muons with energies above a certain threshold $`E_{\mathrm{thr}\mu }`$. The number of such muons may be obtained as follows:

$$N_\mu (E_\mu >E_{\mathrm{thr}\mu },T)=\underset{m}{}\underset{i:E_i>E_{\mathrm{thr}\mu }}{}\underset{ji}{}_0^Th_{mj}(X)\frac{B_m}{E_jX}D_{m\mu }^{ji}W(E_i,X,E_{\mathrm{thr}\mu },T)𝑑X,$$
(19)
where $`D_{m\mu }^{ji}`$ defines the number of muons with energy $`E_i`$ resulting from the decay of hadron $`m`$ with energy $`E_j`$, $`W(E_i,X,E_{\mathrm{thr}\mu },T)`$ is the probability that a muon produced at depth $`X`$ with energy $`E_i`$ will survive between $`X`$ and $`T`$ and its final energy at $`T`$ will be greater than $`E_{\mathrm{thr}\mu }`$.

For the most important channels of muon production ($`\pi ^\pm \mu ^\pm +\nu _\mu `$ and $`K^\pm \mu ^\pm +\nu _\mu `$ are responsible for about 95% of all muons) values of $`D_{m\mu }^{ji}`$ are governed by the simple two body decay kinematics. If we consider the atmosphere to be isotermic then the function $`W(E_i,X,E_{\mathrm{thr}\mu },T)`$ may be written explicitly:

$$W(E_i,X,E_{\mathrm{thr}\mu },T)=\theta \left(E_iE_{\mathrm{thr}\mu }(TX)\frac{dE}{dX}\right)\left\{\frac{X}{T}\frac{E_i(TX)\frac{dE}{dX}}{E_i}\right\}^{B_\mu /(E_i+X\frac{dE}{dX})}$$
(20)
where $`\theta `$ is the step function and $`\frac{dE}{dX}`$ is the ionization loss rate.

## 4 Calculations of the Inclusive Spectra

In this section, we discuss how to obtain the matrices $`W_{mn}^{ji}`$, representing inclusive particle spectra, based on neXus calculations.

The first step evidently amounts to performing a number of neXus simulations to calculate the energy spectrum of produced (secondary) hadrons of type $`h`$,

$$E\frac{dN_{h_{\mathrm{in}}h}}{dE}(E_{\mathrm{in}},E),$$
(21)
in a hadron plus nitrogen reaction $`h_{\mathrm{in}}\mathrm{N}hX`$ for different incident hadrons $`h_{\mathrm{in}}`$, each one at different energies $`E_{\mathrm{in}}`$. For the incident energies, we use

$$E_{\mathrm{in}}^j=10^{2+j/2}\mathrm{GeV},j=0,1,\mathrm{}12,$$
(22)
and for the incident hadron types, we take nucleons, charged pions, charged kaons, and neutral kaons $`K_L^0`$, for the secondary ones in addition neutral pions and $`K_S^0`$ (spectra of $`K_S^0`$ are assumed to be identical to $`K_L^0`$). All other hadrons are assumed to decay immediately.

neXus is based on the Monte Carlo technique, implying automatically statistical fluctuations, which are in particular large, when the production of secondaries with energies approaching the primary one is examined. The influence of the limited statistics obtained via neXus Monte Carlo simulations may be eliminated to some extent if an appropriate smoothing procedure is applied. Such a procedure provides the opportunity to exploit the corresponding ‘smoothened’ spectrum at a fixed primary energy as a continuous function of the energy $`E`$ of secondaries, instead of considering discrete values only. Moreover, it enables one to impose certain restrictions on the shape of inclusive spectra. For example, we assume the $`E`$ dependence of $`W`$ for $`x=E/E_{\mathrm{in}}1`$ to be proportional to $`x^1(1x)^\alpha `$, where $`a`$ is taken from theoretical considerations or just as a fit parameter. We use the Levenberg-Marquardt (LM) method  to fit the Monte Carlo spectra by analytic continuous functions in $`E`$,

$$E\frac{dN_{h_mh_n}}{dE}(E_{\mathrm{in}}^j,E)W_{mn}(E_{\mathrm{in}}^j,E)$$
(23)
for the different incident hadrons $`h_m`$ and secondaries $`h_n`$ for the above-mentioned values $`E_{\mathrm{in}}^j`$ for the incident energy. This method finds the minimum $`\chi ^2`$ and its algorithm consists in the combination of the inverse-Hessian method and the steepest descent method. It is one of the standard non-linear least-square approaches and its detailed description may be found elsewhere . The statistics of $`10^5`$ events for each of the reference energies $`E_{\mathrm{in}}^j`$ is sufficient to calculate $`W_{mn}(E_{\mathrm{in}}^j,E)`$ within $`(0.1÷0.2)\%`$ accuracy.

As an example, we show in figs. 2 and 3

the Monte Carlo results for incident charged pions with energies of $`10^{11}`$ eV and $`10^{17}`$ eV, together with the corresponding fit functions $`W_{mn}(E_{\mathrm{in}}^j,E)`$. Similarly excellent fits are obtained for all the other spectra. For purposes of more efficient computing, two sets of Monte Carlo spectra are used, which exploit different $`E`$-scales: a linear one for large energies of secondaries (right figures) and a logarithmic scale for small ones (left figures).

Using analytic expressions for the spectra, we can now proceed to calculate integrated spectra,

$$W_{mn}^i(E_{\mathrm{in}}^j)=_{E_i/\sqrt{c}}^{E_i\sqrt{c}}\frac{E}{E_i}W_{mn}(E_{\mathrm{in}}^j,E)𝑑E$$
(24)
with $`E_i=E_{\mathrm{min}}\mathrm{c}^\mathrm{i}`$ for $`1ii_{\mathrm{max}}`$. For the actual calculation, we use $`c=10^{0.1}`$ and $`i_{\mathrm{max}}=81`$. We obtain the integrated spectra for arbitrary incident energies $`E_{\mathrm{in}}`$ via an interpolation formula $`\stackrel{~}{W}_{mn}^i(E_{\mathrm{in}})`$, and we can thus calculate these quantities in particular for all the energies $`E_j`$:

$$W_{mn}^{ji}=\stackrel{~}{W}_{mn}^i(E_j).$$
(25)
One could have calculated the integrated spectra directly for all the energies $`E_j`$, but due to the weak energy dependence of the spectra it is much more efficient to proceed as discussed above.

## 5 Solving the Cascade Equations

Having all the ingredients, we are now able to solve the cascade equations as discussed above to obtain the differential hadron spectra $`h_n(E_i,X)`$. The method adopted gives the possibility to calculate average characteristics of EAS within $`1`$% accuracy. Based on the inclusive spectra, we may calculate numbers of different hadrons $`N_{h_i}`$ as well as numbers of electrons $`N_e`$ and $`N_\mu `$ for a number of observation levels $`X`$ in the atmosphere.

In fig. 4, we show the dependencies of the number of electrons $`N_e`$, muons $`N_\mu `$ ($`E_\mu >1`$ GeV), and all hadrons $`N_h`$ ($`E_h>50`$ GeV) on the depth $`X`$ for different incident energies. One observes the expected increase of the particle numbers with energy and as well the shift of the shower maximum towards larger $`X`$ with increasing energy. In fig. 5, we show the corresponding characteristics for individual hadrons. The pions are by far most dominant, followed by nucleons, then charged kaons and $`K_L^0`$, whereas $`K_S^0`$ are the least frequent due to the short life time. In fig. 6, we show the total hadronic energy as a function of the atmospheric depth $`X`$ for different incident energies. Obviously the hadronic energy is highest for the highest incident energy. With increasing atmospheric depth, the hadronic energy drops exponentially, due to its conversion into the energy of electro-magnetic cascade (and to some extent into muon and neutrino energy).

## 6 Some Results of Calculations

It is certainly the main purpose of this paper to introduce a new, sophisticated hadronic interaction model, particularly suited for high energy hadronic interactions, and, at the same time, to explain the use of cascade equations rather than explicit simulations for the air shower calculation. So we do not want to present extensive applications of this approach, but rather discuss one instructive example.

The casa-blanca group  published recently very interesting results concerning the composition of cosmic rays in the energy range $`10^{14}`$-$`10^{17}`$ eV. From the results shown in the previous section, we can easily calculate the shower maximum $`X_{\mathrm{max}}`$ as a function of the incident energy $`E_{\mathrm{in}}`$. We calculate as well the shower maximum for incident nuclei (with mass number $`A`$) assuming that it is given by the result for nucleon at a reduced energy $`E_{\mathrm{in}}/A`$. As it had been shown (see, for example, ) this simple superposition prescription works well for average characteristics of nucleus-induced EAS. In fig. 7

we show the shower maximum $`X_{\mathrm{max}}`$ for incident nucleons (upper curve) and iron (lower curve) as a function of the incident energy $`E_{\mathrm{in}}`$ together with the data. Comparing our calculations with the data we find an excellent agreement with the results of analysis carried out in  on the basis of qgsjet and venus models, i.e. we confirm the change to a heavier primary composition at energies above the knee of the primary cosmic ray spectrum, which has been earlier reported by the group of Moscow State University , as well as the novel feature discovered by the casa-blanca group - “lightening” of the composition at energies just before the knee.

The above statement can be made more quantitative by studying the so-called mean nuclear mass $`A`$, which is defined to be the nuclear mass which would fit best the experimental data. In fig. 8, we plot the mean logarithmic mass ln$`A`$ as a function of the incident energy $`E_{\mathrm{in}}`$ together with qgsjet and venus results obtained in . As already mentioned above, the mass number has clearly a minimum at around $`10^{15.5}`$ eV, for higher energies the mass is increasing again. The neXus results are quite similar to the qgsjet ones, whereas venus has a tendency towards higher masses. One should keep in mind, however, that venus is strictly speaking already outside its energy range of validity, which makes its prediction somewhat uncertain.

It is worth to note that the observed qualitative behaviour for the primary composition may be obtained in the framework of the diffusion model for cosmic ray propagation if one assumes a large magnetic halo for the Galaxy with antisymmetric magnetic field . As it has been shown in  such a configuration of the galactic magnetic field agrees with the measurements and allows to explain naturally the observed “sharpness” of the knee in the primary energy spectrum.

## 7 Summary and Outlook

We introduced a new type of hadronic interaction model (neXus), which has a much more solid theoretical basis when compared, for example, to presently used models such as qgsjet and venus, and so provides much more reliable predictions at superhigh energies where there are no collider data yet. A particular feature of the model is a consistent treatment for calculating cross sections and particle production considering energy conservation strictly in both cases. In addition, one introduces hard processes avoiding any unnatural dependence on artificial cutoff parameters. Using a single set of parameters, neXus is able to fit many basic spectra in proton-proton and lepton-nucleon scattering, as well as in electron-positron annihilation. So it is worth to point out once more that concerning theoretical consistency, neXus is considerably superior to the models in current use, and allows a much safer extrapolation to superhigh energies which are very important in cosmic ray studies.

We explained in detail an air shower calculation scheme, based on cascade equations, which is quite accurate in calculating main characteristics of air showers, being extremely fast concerning computing time. We perform the calculation using the neXus model to provide the necessary tables concerning particle production in hadron-air collisions.

As an application, we calculated the shower maximum as a function of the incident energy for incident protons and iron, to compare with corresponding data. Based on these data, we calculated as well the so-called mean nuclear mass. We are thus able to confirm that the average mass as a function of the incident energy shows a minimum around $`10^{15.5}`$ eV.

Problems where it is possible to ignore fluctuations are not too numerous. So an unavoidable question arises how to retain calculation efficiency at reasonable level and at the same time to account for fluctuations. But the answer to this really crucial question is well known and proves to be rather simple. One needs to employ a combination of the Monte Carlo techhnique and numerical solutions of hadronic cascade equations in the atmosphere. The explicit simulation should be carried out from the primary (initial) energy $`E_{\mathrm{in}}`$ to some threshold $`E_{\mathrm{thr}}=kE_0`$, where $`k10^210^3`$. Numerous calculations showed (see, for example, ) that there is no sense in going below this threshold as it does not increase the accuracy of EAS fluctuation determination. So contributions from hadrons with energies below $`E_{\mathrm{thr}}`$ may be accounted for in average. Numerical solutions of hadronic cascade equations can produce corresponding tables with sufficient accuracy and in a very short time.

Here we do not consider calculations of lateral distributions of shower particles as this problem will be discussed in our next paper. But we would like to point out that there is a rather important class of problems connected with giant EAS arrays aimed at primary energies $`10^{20}`$eV (see ). For these arrays it is sufficient to calculate lateral distributions only at large distances (above 100 m) from the shower axis. Such a situation makes it possible to treat the problem as a combination of the one-dimensional approach for hadrons with energy $`10^{11}`$eV and the rigorous three-dimensional technique for low energy region only.

## 8 Acknowledgements

This work has been funded in part by the IN2P3/CNRS (PICS 580) and the Russian Foundation of Basic Researches (RFBR-98-02-22024).