text:
# On the stability of classical chaotic motion under system’s perturbations

## Abstract

We study in detail the time behavior of classical fidelity for chaotic systems. We show in particular that the asymptotic decay, depending on system dynamical properties, can be either exponential, with a rate determined by the gap in the discretized Perron-Frobenius operator, or algebraic, with the same power as for correlation functions decay. Therefore the decay of fidelity is strictly connected to correlations decay.

As it is known, exponential separation of orbits starting from slightly different initial conditions has been associated to classical chaos. It has been remarked that the situation in quantum mechanics is drastically different. Indeed the scalar product of two states $`\psi _1|\psi _2`$ is time independent. This has led to the introduction of fidelity as a measure of stability of quantum motion . More precisely one considers the overlap of two states which, starting from the same initial conditions, evolve under two slightly different Hamiltonians $`H_0`$ and $`H_ϵ=H_0+ϵV`$. The fidelity is then given by $`f(t)=|\psi |\mathrm{exp}(iH_ϵt/\mathrm{})\mathrm{exp}(iH_0t/\mathrm{})|\psi |^2`$. The quantity $`f(t)`$ can be seen as a measure of the so-called Loschmidt echo: a state $`|\psi `$ evolves for a time $`t`$ under the Hamiltonian $`H_0`$, then the motion is reversed and evolves back for the same time $`t`$ under the Hamiltonian $`H_ϵ`$ and the overlap with the initial state $`|\psi `$ is considered.

However, we would like to stress that, in principle, such difference between classical and quantum mechanics actually does not exist. The Liouville equation, which describes classical evolution, is unitary and reversible as the Schrödinger equation. However, as stressed in several occasions (see, e.g., Ref. ), there exist time scales up to which quantum motion can share the properties of classical chaotic motion including local exponential instability. Due to the existence of such time scales, what may be different, and indeed it is, is the degree of stability of dynamical motion. Indeed, as clearly illustrated in the analysis of Loschmidt echo in , quantum motion turns out to be more stable than the classical motion.

The growing interest in quantum computers has attracted recent interest in this quantity as a measure of the stability of quantum computation in the presence of hardware imperfections or noisy gates operations. Confining ourselves to classically chaotic systems, the emerging picture which results from analytical and numerical investigations  is that both exponential and Gaussian decay are present in the time behavior of fidelity. The strength of the perturbation determines which of the two regimes prevails. The decay rate in the exponential regime appears to be dominated either by the classical Lyapunov exponent or, according to Fermi golden rule, by the spreading width of the local density of states.

In addition, at least for short times, the decaying behavior depends on the initial state (coherent state, mixture, etc.). While it can be true that, for practical purposes, the short time behavior of fidelity may be the most interesting one, it is also true, without any doubt, that in order to have a clear theoretical understanding and identify a possible universal type of quantum decay one needs to consider the asymptotic behavior of fidelity.

This arises the problem of understanding the corresponding classical decay of fidelity and later on inquiring about the times scales at which quantum decay mimics the classical one. In the present paper we concentrate our attention on the classical behavior.

What do we know about the decay of classical fidelity for chaotic systems? What is the relation with correlation functions? Can we derive the decay of fidelity from the behavior of correlations or is fidelity a completely independent function? In a recent paper  it has been found numerically that, after an initial transient, classical fidelity decays exponentially and the rate is given by the Lyapunov exponent (see also ). This is in also in agreement with previous papers  indicating that quantum fidelity, for strong enough quantum perturbation (which, for a fixed classical perturbation strength, corresponds to semiclassical region), decays exponentially with a rate given by the Lyapunov exponent of the corresponding classical system. On the other hand we know that the decay of correlation functions is not ruled by Lyapunov exponent. In the first place there is the general phenomenon of long time tails which means power law decay. In addition, for the special cases in which one can prove exponential decay, the rate is determined by the gap in the discretized Perron-Frobenius operator and not by the Lyapunov exponent.

In this paper, we show that the asymptotic decay of classical fidelity for chaotic systems is not related to the Lyapunov exponent: Similarly to correlation functions, this decay can be either exponential or power law. In the first case, the decay rate is determined by the gap in the discretized Perron-Frobenius operator, in the latter case the power law has the same exponent as for correlation functions.

The classical fidelity $`f(t)`$ is defined as follows:

$$f(t)=_\mathrm{\Omega }𝑑𝐱\rho _ϵ(𝐱,t)\rho _0(𝐱,t),$$
(1)
where the integral is extended over the phase space, and

$$\rho _0(𝐱,t)=U_0^t\rho (𝐱,0),\rho _ϵ(𝐱,t)=U_ϵ^t\rho (𝐱,0)$$
(2)
give the evolution after $`t`$ steps of the initial density $`\rho (𝐱,0)`$ (assumed to be normalized, i.e. $`𝑑𝐱\rho ^2(𝐱,0)=1`$) as determined by the $`t`$-th iteration of the Frobenius-Perron operators $`U_0`$ and $`U_ϵ`$, corresponding to the Hamiltonians $`H_0`$ and $`H_ϵ`$, respectively. The above definition can be shown to correspond to the classical limit of quantum fidelity . In the ideal case of perfect echo ($`ϵ=0`$), the fidelity does not decay, $`f(t)=1`$. However, due to chaotic dynamics, when $`ϵ0`$ the classical fidelity decay sets in after a time scale

$$t_\nu \frac{1}{\lambda }\mathrm{ln}\left(\frac{\nu }{ϵ}\right),$$
(3)
required to amplify the perturbation up to the size $`\nu `$ of the initial distribution, with $`\lambda `$ Lyapunov exponent. Thus, for $`tt_\nu `$ the recovery of initial distribution via the imperfect time-reversal procedure fails, and the fidelity decay is determined by the decay of correlations for a system which evolves forward in time according to the Hamiltonians $`H_0`$ (up to time $`t`$) and $`H_ϵ`$ (from time $`t`$ to time $`2t`$). This is conceptually similar to the “practical” irreversibility of chaotic dynamics: due to the exponential instability, any amount of numerical error in computer simulations rapidly effaces the memory of the initial distribution . In the present case, the coarse-graining which leads to irreversibility is due not to round-off errors but to a perturbation in the Hamiltonian.

In the following, we illustrate this general phenomenon in standard models of classical chaos, characterized by uniform exponential instability (the sawtooth map), marginal stability (the stadium billard), or mixed phase space dynamics (the kicked rotator).

The sawtooth map is defined by

$$\overline{p}=p+F_0(\theta ),\overline{\theta }=\theta +\overline{p},$$
(4)
where $`(p,\theta )`$ are conjugated action-angle variables, $`F_0=K_0(\theta \pi )`$, and the overbars denote the variables after one map iteration. We consider this map on the torus $`0\theta <2\pi `$, $`\pi Lp<\pi L`$, where $`L`$ is an integer. For $`K_0>0`$ the motion is completely chaotic and diffusive, with Lyapunov exponent given by $`\lambda =\mathrm{ln}\{(2+K_0+[(2+K_0)^24]^{1/2})/2\}`$. For $`K_0>1`$ one can estimate the diffusion coefficient $`D`$ by means of the random phase approximation, obtaining $`D(\pi ^2/3)K_0^2`$. In order to compute the fidelity (1), we choose to perturb the kicking strength $`K=K_0+ϵ`$, with $`ϵK_0`$. In practice, we follow the evolution of $`10^8`$ trajectories, which are uniformly distributed inside a given phase space region of area $`A_0`$ at time $`t=0`$. The fidelity $`f(t)`$ is given by the percentage of trajectories that return back to that region after $`t`$ iterations of the map (4) forward, followed by the backward evolution, now with the perturbed strength $`K`$, in the same time interval $`t`$. In order to study the approach to equilibrium for fidelity, we consider the quantity $`g(t)=(f(t)f(\mathrm{}))/(f(0)f(\mathrm{}))`$; in this way $`g(t)`$ drops from $`1`$ to $`0`$ when $`t`$ goes from $`0`$ to $`\mathrm{}`$. We note that $`f(0)=1`$ while, for a chaotic system, $`f(\mathrm{})`$ is given by the ratio $`A_0/A_c`$, with $`A_c`$ the area of the chaotic component to which the initial distribution belongs.

The behavior of $`g(t)`$ is shown in Fig. 1, for $`K_0=(\sqrt{5}+1)/2`$ and different $`L`$ values. One can see that only the short time decay is determined by the Lyapunov exponent. It takes place for $`t_\nu <t<t_ϵ`$, with $`t_\nu `$ defined in Eq. (3) and $`t_ϵ=(1/\lambda )\mathrm{ln}(2\pi /ϵ)`$ time scale required to amplify the effect of the Hamiltonian perturbation up to the maximum extension in the angle $`\theta `$. The Lyapunov regime is followed by a power law decay  $`1/\sqrt{Dt}`$ until the diffusion time $`t_DL^2/D`$ and then the asymptotic relaxation to equilibrium takes place exponentially, with a decay rate $`\gamma `$ (shown in Fig. 2), which, as discussed below, is ruled not by the Lyapunov exponent but by the largest Ruelle-Pollicott resonance . In particular, it is $`ϵ`$-independent.

We determine numerically these resonances for the sawtooth map using the following method :
(i) the phase space torus ($`0\theta <2\pi `$, $`\pi Lp<\pi L`$) is divided in $`N\times NL`$ square cells;
(ii) the transition matrix elements between cells are determined numerically by iterating for one map step the phase space distributions given by the characteristic functions of each cell: in this way we build a finite dimensional approximation of the one-period evolution operator $`U_0`$;
(iii) this truncated evolution matrix $`U_0^{(N)}`$ (of size $`LN^2\times LN^2`$) is diagonalized: it is no longer unitary, and its eigenvalues $`z_i^{(N)}`$ are inside the unit circle in the complex plane (an example is shown in Fig. 3). The non-unitarity of the coarse-grained evolution is due to the fact that the transfer of probability to finer scale structures in the phase space is cut-off, and this results in an effective dissipation ;
(iv) resonances correspond to “frozen” non-unimodular eigenvalues, namely $`z_i^{(N)}z_i`$ when $`N\mathrm{}`$, with $`|z_i|<1`$. Convergence of eigenvalues to values inside the unit circles comes from the asymptotic self-similarity of chaotic dynamics .
As it is known, the asymptotic ($`t\mathrm{}`$) relaxation of correlations is determined by the resonance with largest modulus, $`|\stackrel{~}{z}|=\mathrm{max}_i|z_i|<1`$, giving a decay rate $`\gamma _0=\mathrm{ln}|\stackrel{~}{z}|`$. In Fig. 2 we illustrate the good agreement between the asymptotic decay rate of fidelity (extracted from the data of Fig. 1) and the decay rate as predicted by the gap in the discretized Perron-Frobenius spectrum. It should be stressed that, since fidelity involves forward and backward evolution, the fidelity decay at time $`t`$ has to be compared with the correlations decay at time $`2t`$. For this reason in Fig. 2 the circles correspond to $`\gamma =2\gamma _0`$.

We would like to stress that the same qualitative behavior of Fig. 1 is obtained in the presence of stochastic noise, e.g. the backward evolution is driven by a kicking strength $`K(t)=K_0+ϵ(t)`$, with $`\{ϵ(t)\}_{t=1,2,\mathrm{}}`$ uniformly and randomly distributed inside the interval $`[ϵ,ϵ]`$. In particular, we observed the initial Lyapunov decay and the asymptotic exponential relaxation with the same rate $`\gamma `$. This means that the effect of a noisy environment on the decay of fidelity for a classical chaotic system is similar to that of a generic static Hamiltonian perturbation.

Further confirmation for the validity of the above illustrated scenario has been obtained by analyzing systems in which the asymptotic decay of correlations is algebraic. This happens in two cases:
(i) when the system possesses marginally stable orbits. A typical example is the stadium billiard in which, as it is known , correlations decay as $`1/t`$;
(ii) when there is mixed phase space . A typical example is the kicked rotator model (described by Eq. (4) with $`F_0=K_0\mathrm{sin}\theta `$).
Since in the long time limit the fidelity decay at time $`t`$ is still related to the decay of correlations at time $`2t`$, in the case of power law decay of correlations as $`t^\alpha `$, we expect a power law decay of fidelity with the same exponent $`\alpha `$. This is indeed confirmed by our numerical results. In Fig. 4 it is shown that, for the stadium billiard, fidelity decays asymptotically $`1/t`$, as expected. In Fig. 5 we compare the fidelity decay (at time $`t`$) and the correlations decay (at time $`2t`$) for the kicked rotator with kicking parameter $`K_0=2.5`$, for which the phase space contains chaotic components and stable islands. The correlator is given by $`D(t)=(C(t)C(\mathrm{}))/(C(0)C(\mathrm{}))`$, with $`C(t)=_\mathrm{\Omega }𝑑𝐱\rho _0(𝐱,t)\rho (𝐱,0)`$. It is seen that, after an initial Lyapunov decay, fidelity approaches the same asymptotic power law decay of correlations .

In summary, we have shown that in chaotic systems the asymptotic decay of classical fidelity, which describes the structural stability of motion under system’s perturbations, is analogous to the asymptotic decay of correlation functions. This asymptotic decay can be either exponential or algebraic, depending on the dynamical properties of the system. In any instance, it is not related to the local exponential instability ruled by the Lyapunov exponent and it is $`ϵ`$-independent. It would be interesting to understand what are the implications of these findings for the decay of quantum fidelity.

We thank Roberto Artuso for useful discussions. This work was supported in part by the EC RTN contract HPRN-CT-2000-0156, the NSA and ARDA under ARO contract No. DAAD19-02-1-0086, the PA INFM “Weak chaos: Theory and applications”, and the PRIN 2000 “Chaos and localization in classical and quantum mechanics”.


format: text

text:
# 1 The diffusion coefficient vs the tilting force for 𝑘=0.2. Curve 1 - 𝑇=0.01, curve 2 - 𝑇=0.1, curve 3 - 𝑇=0.5

## Acknowledgement

The authors are grateful to Prof. Romi Mankin and Dmitri Martila, M.Sci., for useful discussions and assistance.


text:
# Radiative Neutralino Decay in Supersymmetric Models

## 1 Introduction

One of the major goals of collider experiments is to verify or disprove the existence of weak scale supersymmetric matter. Given sufficient energy and integrated luminosity, sparticles should be produced at large rates, if indeed they exist. Once produced, sparticles will decay via a cascade into other sparticles until the state containing the lightest SUSY particle (LSP) is reached. If kinematically accessible, tree level two-body decays will comprise the dominant sparticle branching fractions. If two-body modes are suppressed or closed, then tree level three-body decays can dominate. It is also possible for sparticles to decay into two-body final states where the decay is mediated by loop diagrams. Important examples include $`\stackrel{~}{g}g\stackrel{~}{Z}_i`$, $`\stackrel{~}{Z}_i\stackrel{~}{Z}_j\gamma `$ ($`i,j=14`$ with $`i>j`$) and $`\stackrel{~}{t}_1c\stackrel{~}{Z}_1`$ for sparticles, and $`h\gamma \gamma `$ for the lightest SUSY Higgs scalar.

The radiative decay $`\stackrel{~}{Z}_i\stackrel{~}{Z}_j\gamma `$ has been examined previously in a number of papers, since it can give rise to a distinctive energetic, isolated photon plus missing energy signature. The decay proceeds through loops containing intermediate charged scalars and charginos. After several approximate calculations, the complete calculation was presented in both the non-linear and linear $`R_\xi `$ gauge by Haber and Wyler. These authors note that the radiative neutralino decay branching fraction may be large in MSSM parameter space regions where $`M_1,M_2,|\mu |M_Z`$ (light neutralino radiative decay), or when $`|\mu |,M_ZM_1,M_2`$ TeV (higgsino to higgsino radiative decay).

Detailed numerical studies of the $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ decay rate have been performed by Ambrosanio and Mele. These authors found that the $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ decay rate can in fact be the dominant $`\stackrel{~}{Z}_2`$ decay mode in regions of MSSM model parameter space where tree level two-body $`\stackrel{~}{Z}_2`$ decay modes are not allowed, and the competing $`\stackrel{~}{Z}_2`$ three-body decay modes are dynamically and/or kinematically suppressed.

At low $`\mathrm{tan}\beta `$, the three-body neutralino decay $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1f\overline{f}`$ occurs dominantly via five Feynman diagrams involving virtual sfermion and $`Z`$ boson exchange. The dynamical suppression of the three-body decay modes occurs when one of $`\stackrel{~}{Z}_2`$ and $`\stackrel{~}{Z}_1`$ is gaugino-like, while the other is higgsino-like. Since the $`\stackrel{~}{Z}_if\stackrel{~}{f}`$ vertex contains the gaugino component of the neutralinos, the decay via a virtual sfermion $`\stackrel{~}{f}^{}`$ is suppressed if one (or both) of $`\stackrel{~}{Z}_2`$ or $`\stackrel{~}{Z}_1`$ is higgsino-like. In addition, the $`Z\stackrel{~}{Z}_2\stackrel{~}{Z}_1`$ coupling contains the higgsino components of the neutralinos so that if one of $`\stackrel{~}{Z}_1`$ and/or $`\stackrel{~}{Z}_2`$ is gaugino-like, then decay via the $`Z^{}`$ diagram is suppressed. The region of MSSM parameter space where dynamical suppression of three-body decay modes occurs is when the gaugino masses $`M_1M_2`$, $`\mathrm{tan}\beta 1`$ and $`\mu <0`$ with $`|\mu |`$ being sufficiently small. In light of recent limits from LEP2 that $`m_{\stackrel{~}{W}_1}>103.5`$ GeV and $`m_h>114.1`$ GeV (for a SM-like $`h`$ boson), the region of dynamical suppression is largely excluded, since $`m_h`$ is quite light at very low values of $`\mathrm{tan}\beta `$.

In addition, there can exist a region of reduced dynamical suppression. In this case, three-body decays via $`Z^{}`$ can be suppressed if one (or both) of the neutralinos is gaugino-like (thus suppressing $`\stackrel{~}{Z}_2`$ decay via $`Z^{}`$), while the sfermions are quite heavy (thus suppressing decay through $`\stackrel{~}{f}^{}`$). This interesting case occurs naturally as we shall see in the “focus-point” scenario of the minimal supergravity (mSUGRA or CMSSM) model.

Finally, there can exist regions of kinematical suppression of three-body $`\stackrel{~}{Z}_2`$ decays. In this case, $`m_{\stackrel{~}{Z}_2}m_{\stackrel{~}{Z}_1}`$ so that $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1f\overline{f}`$ is suppressed kinematically by a factor $`ϵ^5`$, where $`ϵ=(1m_{\stackrel{~}{Z}_1}/m_{\stackrel{~}{Z}_2})`$. The radiative decay $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ is suppressed by $`ϵ^3`$, so that as the $`m_{\stackrel{~}{Z}_2}m_{\stackrel{~}{Z}_1}`$ mass gap approaches zero, the radiative decay increasingly dominates. The kinematic suppression region occurs when $`\mathrm{tan}\beta 1`$ or when $`M_1M_2`$. In the region of parameter space where kinematical suppression occurs, however, the small mass gap between neutralinos also means that the final state photon energy will be small, unless enhanced by a Lorentz boost if $`\stackrel{~}{Z}_2`$ is produced either directly or via cascade decay with a high velocity.

The radiative neutralino decays have received some recent attention by Kane et al. in attempts to explain an $`ee\gamma \gamma `$ event seen by the CDF experiment in run 1 of the Fermilab Tevatron. These “higgsino-world” scenarios seem to need $`\mathrm{tan}\beta 1`$, $`\mu <0`$ and a light top squark. These scenarios seem largely excluded now by LEP2 limits on the masses of the Higgs bosons, which require $`\mathrm{tan}\beta >23`$.

In this paper, we re-examine the radiative neutralino decay branching fraction as a function of supersymmetric model parameter space. Our calculations include a number of improvements over previous analyses.

1. We implement the complete $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ decay calculation into the event generator ISAJET v7.64. Our calculation includes full third generation mixing effects amongst the squarks and sleptons.
2. Our treatment of neutralino three-body decays includes the effect of all third generation Yukawa couplings and sfermion mixings in the decay calculations. These are especially important at large values of the parameter $`\mathrm{tan}\beta `$, where the Yukawa couplings $`f_b`$ and $`f_\tau `$ can become large. In addition, along with the five $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1f\overline{f}`$ decay diagrams proceeding via sfermion and $`Z`$ boson exchange, we include decays through $`h`$, $`H`$ and $`A`$ bosons. These diagrams can be important, especially at large $`\mathrm{tan}\beta `$.
3. We include the latest constraints from LEP2. As noted before, these constraints eliminate many of the regions of parameter space where $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ can be large, especially around $`\mathrm{tan}\beta 1`$.
4. We present results for the paradigm mSUGRA model. We also present a scan of MSSM parameter space which includes negative as well as positive gaugino masses. We find the radiative decay branching fraction can be large in a particular SUGRA model with non-univeral gaugino masses.
5. Finally, in our numerical scans, we take note of parameter space regions where the radiative neutralino branching fraction is not dominant, but may be nonetheless observable. For instance, at the CERN LHC, many millions of sparticle production events may be recorded each year, provided sparticles are light enough. In this case, branching fractions of a few per cent or even less may be interesting, and allow one to measure the neutralino radiative decay branching fraction. Since the $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ branching fraction is very sensitive to parameter space, its measurement may help determine or constrain some of the fundamental SUSY model parameters.

The rest of this paper is organized as follows. In Sec. 2, we examine the radiative neutralino branching fraction in the mSUGRA model. We find it can reach levels approaching 1%, especially in the well-motivated focus-point region. In Sec. 3, we examine the radiative decay branching fraction in the MSSM. In this case, we find regions of parameter space where the branching fraction can approach 50-100%, even for large values of $`\mathrm{tan}\beta `$. These regions occur dominantly where $`M_1M_2`$ or when $`M_1M_2`$. We note in Sec. 4 that the latter MSSM regions of enhancement naturally occur in a particular SUGRA model with non-universal gaugino masses. We map out the interesting regions of parameter space where the radiative decay branching fraction can be large. Finally, in an appendix, we present our formulae for radiative neutralino production, cast in a form consistent with the ISAJET event generator program.

## 2 $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ branching fraction in mSUGRA model

Our first goal is to examine the radiative neutralino branching fraction in the paradigm mSUGRA (or CMSSM) model. Most analyses of supersymmetric phenomena take place within this model, so it is convenient to know the rates for radiative neutralino decay in relation to other expected phenomena. In the mSUGRA model, the MSSM is assumed to be valid between the grand unified scale ($`M_{GUT}2\times 10^{16}`$ GeV) and the weak scale ($`1`$ TeV). At $`M_{GUT}`$, it is assumed all scalar masses are unified to the parameter $`m_0`$, while all gaugino masses unify to $`m_{1/2}`$. In addition, trilinear soft SUSY breaking masses unify to $`A_0`$. The gauge and Yukawa couplings and soft SUSY breaking parameters evolve from $`M_{GUT}`$ to $`M_{weak}`$ according to the well-known MSSM renormalization group equations (RGEs). At $`M_{weak}`$ it is assumed electroweak symmetry is broken radiatively, which fixes the magnitude of the superpotential $`\mu `$ parameter, and allows one to effectively trade the bilinear soft breaking parameter $`B`$ for the ratio of Higgs vevs $`\mathrm{tan}\beta `$. Thus the parameter space is characterized by

$$m_0,m_{1/2},A_0,\mathrm{tan}\beta \mathrm{and}sign(\mu ),$$
(1)
where we take the top quark pole mass $`m_t=175`$ GeV. In our mSUGRA model calculations, we utilize the iterative RGE solution embedded within the program ISAJET v7.64. This version also includes our radiative neutralino decay calculation.

In Fig. 1, we show the $`m_0vs.m_{1/2}`$ plane for $`\mathrm{tan}\beta =10`$, $`A_0=0`$ and $`\mu >0`$. The red-shaded regions are excluded either by the presence of a stau LSP (upper left) or a lack of REWSB (lower and right). The pink shaded region is excluded by LEP2 searches for charginos and a light Higgs scalar. The upper and left regions indicate where $`\stackrel{~}{Z}_2`$ two-body decays occur: either to $`\stackrel{~}{\mathrm{}}_1\mathrm{}`$ ($`\mathrm{}`$ denotes any of the charged or neutral leptons) on the left, or to $`\stackrel{~}{Z}_1Z`$ or $`\stackrel{~}{Z}_1h`$ in the upper plane. The remaining regions allowed by LEP2 have dominant three-body $`\stackrel{~}{Z}_2`$ decays. The green shaded region has $`0.001<BF(\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma )<0.005`$ while the yellow has $`0.005<BF(\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma )<0.01`$. In these regions, where $`m_{1/2}\stackrel{<}{}250`$ GeV, we also have $`m_{\stackrel{~}{g}}\stackrel{<}{}650`$ GeV. Then $`\stackrel{~}{g}\stackrel{~}{g}`$ cross sections exceed $`2000`$ fb at the CERN LHC, leading to at least 20,000 SUSY events per 10 fb<sup>-1</sup> of integrated luminosity. After suitable cuts are made on number of jets, jet $`E_T`$ and $`\overline{)}E_T`$, a nearly pure sample of SUSY events will remain. These SUSY candidates can then be scrutinized for the presence of hard, isolated photons, and may allow a determination at some level of the neutralino radiative decay branching fraction.

We note that the regions of mSUGRA parameter space with the largest $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ branching fraction occur at rather large $`m_0`$ and low $`m_{1/2}`$ values. The very largest allowed $`m_0`$ areas correspond to the focus-point region of Feng, Matchev and Moroi, where fine-tuning may be low, but squark and slepton masses are at the multi-TeV level. These heavy sfermion masses may be sufficient to suppress many flavor-changing or $`CP`$ violating processes that could arise due to the presence of off-diagonal flavor changing soft SUSY masses, or imaginary components of soft SUSY breaking terms.

We also show the regions of $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ branching fraction for large $`\mathrm{tan}\beta =45`$ in Fig. 2. Again, we see a significant part of the three-body $`\stackrel{~}{Z}_2`$ decay region inhabited with possibly measurable branching fractions for radiative neutralino decay. The region with the largest branching fraction is adjacent to the right-side region forbidden by lack of REWSB. In that region, the $`\mu `$ parameter is becoming small, and the light neutralinos are becoming increasingly higgsino-like. In addition, sfermion masses are heavy, so we are in a region of reduced dynamical suppression of three-body decay modes, also with some kinematic suppression.

In both the low and the high $`\mathrm{tan}\beta `$ cases, it is important to note that the $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ branching fraction reaches no greater than 1%, so is never dominant. However, regions do occur where the radiative decay should be detectable, and the branching fraction should be measureable. Detailed simulations including detector effects for identifying isolated photons will be necessary to ascertain the precision with which such measurements can be made.

## 3 $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ branching fraction in the MSSM

It is apparent from the numerical studies of Ambrosanio and Mele that the greater freedom of parameter choices in the unconstrained MSSM can allow for regions of parameter space with larger $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ branching fractions than in the mSUGRA model. In this section, we present a scan over unconstrained MSSM parameters. To make the analysis tractable, we do assume all MSSM soft SUSY breaking parameters and Yukawa couplings to be purely real, and we take all sfermion masses to be degenerate, with no intergenerational mixing. For simplicity, we take $`A_t=A_b=A_\tau `$, and fix $`m_{\stackrel{~}{g}}=1500`$ GeV. Our scan of MSSM parameter space is then over the following region:

$`1500\mathrm{GeV}<`$ $`M_1`$ $`<1500\mathrm{GeV},`$
$`0<`$ $`M_2`$ $`<1500\mathrm{GeV},`$
$`1500\mathrm{GeV}<`$ $`\mu `$ $`<1500\mathrm{GeV},`$
$`2000\mathrm{GeV}<`$ $`A_t`$ $`<2000\mathrm{GeV},`$
$`150\mathrm{GeV}<`$ $`m_A`$ $`<1500\mathrm{GeV},`$
$`2<`$ $`\mathrm{tan}\beta `$ $`<50,`$
$`100\mathrm{GeV}<`$ $`m_{\stackrel{~}{f}}`$ $`<1500\mathrm{GeV},`$
where we have the freedom to choose one of the gaugino masses, in this case $`M_2`$, to be real and positive. We require $`m_{\stackrel{~}{W}_1}>103.5`$ GeV and $`m_h>114.1`$ GeV, according to LEP2. The latter Higgs mass requirement is overly harsh in that the limit on the lightest SUSY Higgs boson is somewhat lower, but this doesn’t alter our conclusions. We also require a mass gap $`m_{\stackrel{~}{Z}_2}m_{\stackrel{~}{Z}_1}>5`$ GeV, so that the photon that appears in scattering events from any of these models has a detectable energy.

Our results are presented in Figs. 3 and 4, as plots of branching fraction ($`BF`$) versus model parameter. Fig. 3a and 3b show the branching fractions versus the gaugino masses $`M_1`$ and $`M_2`$. We see immediately that in fact there do exist models with substantial radiative decay branching fractions in the range 10%-100%. For both gaugino masses, these large branching fraction models populate the region with $`|M_1|,M_2<1000`$ GeV, with the distribution of points being somewhat asymmetric for $`M_1`$ being positive or negative. If $`M_1`$ becomes large negative or $`M_2`$ becomes large, the branching fractions can remain large, but the mass gap $`m_{\stackrel{~}{Z}_2}m_{\stackrel{~}{Z}_1}`$ decreases below 5 GeV. However, for large positive $`M_1`$ the mass gap stays large, but in this case we do not get large branching fractions because either two-body decay modes of $`\stackrel{~}{Z}_2`$ turn on, or else $`\stackrel{~}{Z}_2\stackrel{~}{W}_1f\overline{f}^{}`$ becomes large, thus supressing the radiative decay branching fraction. The models with $`|M_1|,M_2`$ near zero are usually excluded by constraints from LEP on $`Z\stackrel{~}{Z}_i\stackrel{~}{Z}_j`$ decay or from LEP2 on the chargino mass.

Fig. 3c shows the $`BF`$ versus the $`\mu `$ parameter. Again, small values of $`|\mu |`$ are excluded by LEP and LEP2 constraints. However, we see that really large radiative decay branching fractions only occur for very large values of $`|\mu |\stackrel{>}{}300400`$ GeV. The suppresion of lower values of $`|\mu |`$ comes from the mass gap requirement. Fig. 3d shows the branching fraction versus the parameter $`A_t`$. Here, large values of branching fraction are possible for all values of $`A_t`$, although the very largest seem to favor large values of $`|A_t|`$.

In Fig. 4 we show the $`BF`$ versus the parameters $`m_A`$, $`\mathrm{tan}\beta `$ and $`m_{\stackrel{~}{f}}`$ in frames a, b and c. We see that large branching fractions are possible for all values of pseudoscalar mass $`m_A`$, and for all values of $`\mathrm{tan}\beta `$. Large values of BF have some preference for $`m_{\stackrel{~}{f}}\stackrel{>}{}300400`$ GeV, so that two-body decays $`\stackrel{~}{Z}_2f\stackrel{~}{f}`$ are closed, and three-body decays $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1f\overline{f}`$ through an intermediate $`\stackrel{~}{f}`$ have some suppression.

In the last frame, Fig. 4d, we show the branching fraction not against an input parameter, but against the neutralino mass gap, $`m_{\stackrel{~}{Z}_2}m_{\stackrel{~}{Z}_1}`$. Here it is clear that the largest values of $`BF`$ prefer the smallest mass gaps, where there exists at least some amount of kinematical suppression of three-body $`\stackrel{~}{Z}_2`$ decays. The large branching fractions have a cut-off around 90-100 GeV, where $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1Z`$ decays become possible.

To gain more understanding of the parameter space where radiative decays are large, we plot in Fig. 5 only models with $`BF>10\%`$ (black dots) or $`>1\%`$ (red dots), in correlated parameter planes. Again, we require $`m_{\stackrel{~}{Z}_2}m_{\stackrel{~}{Z}_1}>5`$ GeV, and all LEP and LEP2 constraints to be satisfied. From Fig. 5a in the $`M_1vs.M_2`$ plane, we see that models with the largest $`BF`$ require $`|M_1|M_2`$, a feature noted by Ambrosanio and Mele in their analysis. The region with $`M_1M_2`$ turns out to be quite narrow, while the region with $`M_1M_2`$ has slightly more breadth to it. This latter point is even more apparent when we relax to the condition $`BF>1\%`$.

In Fig. 5b and 5c, we show models in the $`M_1vs.\mu `$ and $`M_2vs.\mu `$ planes. For models with $`BF>10\%`$, evidently $`|M_1|<|\mu |`$ and also $`M_2<|\mu |`$. If we relax to models with $`BF>1\%`$, then these requirements no longer hold. In frame 5d, we show models in the $`M_2/M_1vs.\mu `$ plane. Models with $`BF>10\%`$ clearly favor $`M_2/M_11`$, with either sign of $`\mu `$ being equally preferable. If we allow models with $`BF>1\%`$, then the $`M_1\pm M_2`$ restriction can be relaxed, and rather large ratios of $`M_2/M_1`$ are allowed, especially if $`|\mu |`$ is small.

An important element of our analysis is the size of the $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1`$ mass gap, so there is sufficient photon energy for observations. In Fig. 6, we show the mass gap $`m_{\stackrel{~}{Z}_2}m_{\stackrel{~}{Z}_1}`$ versus the parameters $`\mu `$, $`M_1`$, $`M_2`$ and $`\mathrm{tan}\beta `$ in frames a-d, respectively. From Fig. 6a, we see that models with $`BF>10\%`$ prefer large absolute values of the $`\mu `$ parameter, but for these models the mass gap is bounded typically by 40-60 GeV. By relaxing to models with $`BF>1\%`$, we pick up models with mass gaps ranging to 100 GeV. However, for models with small values of $`\mu `$, where higgsino$``$ higgsino transitions are possible, the mass gap is very small ($`\stackrel{<}{}20`$ GeV), and these cases will be more difficult to observe. Figs. 6b and 6c show that the mass gap for models with $`BF>10\%`$ is typically less than 30-40 GeV for $`|M_1|,M_2`$ ranging up to 1000 GeV. The models with $`BF>1\%`$ and $`|M_1|,M_2>1000`$ GeV typically have a small mass gap, less than 10-15 GeV. Finally, in Fig. 6d, the mass gap versus $`\mathrm{tan}\beta `$ shows no major correlations, and significant mass gaps can be achieved for any $`\mathrm{tan}\beta `$ values.

To conclude this section, we summarize our findings by noting that the largest neutralino radiative decay branching fractions are found in regions of parameter space with

$`|M_1|,M_21001000\mathrm{GeV},`$
$`|\mu |>|M_1|,M_2\mathrm{and}`$
$`|M_1|M_2.`$
## 4 $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ branching fraction in a SUGRA model with non-universal gaugino masses

So far, our findings are that in the mSUGRA model the radiative decay $`BF`$ never exceeds 1%, and is largest in regions with large $`m_0`$ such as the focus point region. For the MSSM case, we found distinct regions of parameter space where the radiative branching fraction can range up to 10%-100%, although these conditions are never realized in the mSUGRA model. In this section, we note that there do exist models where the above MSSM parameter space restrictions are realized, namely in SUGRA models with non-universal gaugino masses.

Since supergravity is not a renormalizable theory, in general we may expect a non-trivial gauge kinetic function. Expanding the gauge kinetic function as

$$f_{AB}=\delta _{AB}+\widehat{h}_{AB}/M_P+\mathrm{},$$
(2)
where the fields $`\widehat{h}_{AB}`$ transform as left handed chiral superfields under supersymmetry transformations, and as the symmetric product of two adjoints under gauge symmetries, we parametrize the lowest order contribution to gaugino masses by,

$$\frac{1}{4}d^2\theta _L\frac{\widehat{h}_{AB}}{M_P}\overline{\widehat{W}_A^c}\widehat{W}_B\frac{F_h_{AB}}{M_P}\lambda _A\lambda _B+\mathrm{}$$
(3)
where the $`\widehat{W}_A`$ are curl superfields, the $`\lambda _A`$ are the gaugino fields, and $`F_{hAB}`$ are the auxillary field components of the $`\widehat{h}_{AB}`$ that acquire a SUSY breaking $`vev`$.

If the fields $`F_h`$ which break supersymmetry are gauge singlets, universal gaugino masses result. However, in principle, the chiral superfield which communicates supersymmetry breaking to the gaugino fields can lie in any representation in the symmetric product of two adjoints, and so can lead to gaugino mass terms that (spontaneously) break the underlying gauge symmetry.

In the context of $`SU(5)`$ grand unification, $`F_h`$ belongs to an $`SU(5)`$ irreducible representation which appears in the symmetric product of two adjoints:

$$(\mathrm{𝟐𝟒}\times \mathrm{𝟐𝟒})_{\mathrm{symmetric}}=\mathrm{𝟏}\mathrm{𝟐𝟒}\mathrm{𝟕𝟓}\mathrm{𝟐𝟎𝟎},$$
(4)
where only $`\mathrm{𝟏}`$ yields universal masses. The relations amongst the various $`GUT`$ scale gaugino masses have been worked out by Anderson. The relative $`GUT`$ scale $`SU(3)`$, $`SU(2)`$ and $`U(1)`$ gaugino masses $`M_3`$, $`M_2`$ and $`M_1`$ are listed in Table 1 along with the approximate masses after RGE evolution to $`QM_Z`$. These scenarios represent the predictive subset of the more general (and less predictive) case of an arbitrary superposition of these representations. The model parameters may be chosen to be,

$$m_0,M_3^0,A_0,\mathrm{tan}\beta \mathrm{and}sign(\mu ),$$
(5)
where $`M_i^0`$ is the $`SU(i)`$ gaugino mass at scale $`Q=M_{GUT}`$. $`M_2^0`$ and $`M_1^0`$ can then be calculated in terms of $`M_3^0`$ according to Table 1, and the weak scale SUSY particle spectrum can again be calculated using ISAJET v7.64.

We note that the 75 model leads to weak scale gaugino masses in the ratio $`M_2:M_16:5`$, which is approaching the conditions required in Sec. 3 for models with large radiative branching fractions. We adopt the GUT scale gaugino mass boundary conditions of the 75 model, and explore the model parameter space for large radiative decay branching fractions.

Our first results are shown in Fig. 7 in the $`M_3^0vs.m_0`$ plane, for $`\mathrm{tan}\beta =5`$, $`A_0=0`$ and $`\mu >0`$. Again, the red shaded regions are excluded theoretically, while the pink regions are excluded by constraints from LEP2. The region with mass gap $`m_{\stackrel{~}{Z}_2}m_{\stackrel{~}{Z}_1}>5`$ GeV is indicated between the dotted contours. The substantial green region has branching fractions $`BF:1\%10\%`$, while the yellow region has branching fractions in the range $`BF:10\%50\%`$. Almost all the large branching fraction region has mass gap $`m_{\stackrel{~}{Z}_2}m_{\stackrel{~}{Z}_1}<5`$ GeV. The value of $`|\mu |`$ in the yellow region is generally much smaller than $`|M_1|`$ and $`M_2`$, so that $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ consists of a higgsino $``$ higgsino transition, with a rather small mass gap between $`\stackrel{~}{Z}_2`$ and $`\stackrel{~}{Z}_1`$; there is significant kinematic suppression of three-body decay rates. In the lower right unshaded regions, the value of $`\mu `$ is comparable or larger than $`|M_1|`$ and $`M_2`$, and the $`\stackrel{~}{Z}_2`$-$`\stackrel{~}{Z}_1`$ mass gap is usually larger, so three-body decay rates increase, suppressing the radiative decay branching fraction. An exception occurs near $`m_03000`$ GeV and $`M_3^0300`$ GeV where the light neutralinos are transitioning between higgsino-like and gaugino-like. In this area, the mass gap drops to very small values of order 1 GeV, and the radiative decay is momentarily enhanced.

In Fig. 8, we show again the 75 model, but this time with $`\mathrm{tan}\beta =10`$. As $`\mathrm{tan}\beta `$ is increased beyond 5, the allowed parameter space rapidly decreases. In this case, only a small region in the center of the plot gives a viable SUSY spectrum satsifying all constraints. In this case, the yellow region represents $`1\%<BF<5\%`$ (lower values than in Fig. 7). The larger $`\mathrm{tan}\beta `$ value allows for a larger mass gap between $`\stackrel{~}{Z}_2`$ and $`\stackrel{~}{Z}_1`$ for small $`\mu `$ values compared to Fig. 7 ($`m_{\stackrel{~}{Z}_2}m_{\stackrel{~}{Z}_1}>5`$ GeV throughout the $`\mathrm{tan}\beta =10`$ plane). Consequently, three-body decays suffer less kinematical suppression, and the radiative branching fraction is lower.

We have also examined briefly two cases of “more minimal” supersymmetric models. In these models, also known as inverted mass hierarchy (IMH) models, there can exist non-universal scalar masses, leading to multi-TeV first and second generation scalars (to suppress flavor and CP violating SUSY processes), while third generation scalars remain in the sub-TeV regime (as required by naturalness constraints). The two models examined were the radiative inverted hierarchy model (RIMH) and the GUT scale IMH model (GSIMH). In both these models, usually two-body $`\stackrel{~}{Z}_2`$ decays are allowed. In the RIMH model, for some parameter regions the $`\mu `$ parameter can become very small. In that case, results are similar to the SUGRA model focus point region already presented with branching fractions near the $`1\%`$ level.

## 5 Conclusions

We have re-calculated the decay rate for radiative neutralino decay within the MSSM. We have also made detailed scans of parameter space in the mSUGRA model, the MSSM and a SUGRA model with non-universal gaugino masses. In the mSUGRA model, radiative neutralino branching fractions for $`\stackrel{~}{Z}_2`$ can reach up to 1%. While this rate comprises only a small fraction of the $`\stackrel{~}{Z}_2`$ decay modes, it may be sufficiently large to be measureable at collider experiments at the CERN LHC, where many $`\stackrel{~}{g}`$ and $`\stackrel{~}{q}`$ events can be produced, and $`\stackrel{~}{Z}_2`$ will occur in cascade decays. The rate is largest in the focus point region of parameter space. If the branching fraction is observed, then it may aid in the determination of underlying model parameters.

In the MSSM, the $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ decay rate can reach nearly 100%. In this case, the gaugino masses $`|M_1|,M_2\stackrel{<}{}1000`$ GeV, with $`|M_1|,M_2<|\mu |`$ and $`|M_1|M_2`$. The branching fraction can be large for all values of $`\mathrm{tan}\beta `$.

The conditions on parameters of the MSSM for a large radiative neutralino branching fraction turn out to be realized in a SUGRA model with non-universal gaugino masses, where the hidden sector field $`\widehat{h}_{AB}`$ transforms as a 75 dimensional symmetric product of two adjoints of $`SU(5)`$. At low $`\mathrm{tan}\beta 5`$, in this case the branching fraction can reach up to 50%, although the mass gap $`m_{\stackrel{~}{Z}_2}m_{\stackrel{~}{Z}_1}`$ tends to be rather small - of order a few GeV. For higher $`\mathrm{tan}\beta `$, the parameter space in the 75 model rapidly diminishes. However, scans made for the $`\mathrm{tan}\beta =10`$ case show the radiative branching fraction reaching up to $`5\%`$.

Our final conclusion is that once a significant sample of SUSY particle candidate events is discovered at collider experiments, then that sample should be scrutinized for the presence of additional isolated photons, beyond the level expected from bremsstrahlung, misidentified jets, etc. If the events contain in addition isolated photons, then it may be possible to test whether they come from radiative neutralino decay, and to in fact measure that particular branching fraction. If the $`\stackrel{~}{Z}_2\stackrel{~}{Z}_1\gamma `$ branching fraction is in fact measured, then its value should aid in the attempt to extract the values of underlying model parameters.

## Appendix

In this appendix, we present our formulae for the radiative neutralino decay rate. The complete radiative decay width of $`\stackrel{~}{Z}_j`$ was calculated by Haber and Wyler. We organize our calculation in a similar fashion to theirs, working in the non-linear $`R_\xi `$ gauge. We use the Lagrangian conventions of Tata, which is useful for inclusion into ISAJET. The relevant Feynman diagrams labelled a-p are presented in Ref. , and will not be repeated here. We have checked that our formula reproduces analytically that of Haber and Wyler after converting conventions.

In addition, we have compared our numerical results for various diagram subsets with a calculation by Ambrosanio and Mele which is embedded in the SUSYGEN 3.00/36 program. Again, we find excellent agreement, except for the case of sfermion loop diagrams, where mixing effects cause our result to differ from those of SUSYGEN. If we neglect mixing effects, then we agree with the SUSYGEN result.

The radiative decay width is given by

$$\mathrm{\Gamma }(\stackrel{~}{Z}_j\stackrel{~}{Z}_i\gamma )=\frac{g_{\stackrel{~}{Z}_j\stackrel{~}{Z}_i\gamma }^2}{8\pi }\left(\frac{m_{\stackrel{~}{Z}_j}^2m_{\stackrel{~}{Z}_i}^2}{m_{\stackrel{~}{Z}_i}}\right)^3.$$
(6)
Here

$$g_{\stackrel{~}{Z}_j\stackrel{~}{Z}_i\gamma }=g_{abcd}+g_{efgh}+g_{ijkl}+g_{mnop},$$
(7)
with the separate contributions given by the following:

$$g_{abcd}=\underset{sfermions}{}\frac{eq_f}{(4\pi )^2}C_f\left(m_{\stackrel{~}{Z}_i}C_i^{abcd}K+m_{\stackrel{~}{Z}_j}C_j^{abcd}(I^2K)+m_fC_0^{abcd}I\right).$$
(8)
Here $`C_f=1`$ for leptons and $`C_f=3`$ for quarks.

$`C_i^{abcd}`$ $`=`$ $`(i)^{\theta _i}(i)^{\theta _j}\left(\alpha _k^{(i)}\alpha _k^{(j)}\beta _k^{(i)}\beta _k^{(j)}\right),`$
$`C_j^{abcd}`$ $`=`$ $`(i)^{\theta _i}(i)^{\theta _j}\left(\alpha _k^{(i)}\alpha _k^{(j)}\beta _k^{(i)}\beta _k^{(j)}\right),`$
$`C_0^{abcd}`$ $`=`$ $`(i)^{\theta _i+\theta _j}\left(\beta _k^{(i)}\alpha _k^{(j)}\alpha _k^{(i)}\beta _k^{(j)}\right),`$
$`\alpha _{\stackrel{~}{f}_1}^{(i)}`$ $`=`$ $`A_f^{(i)}\mathrm{cos}\theta _ff_fv_a^{(i)}\mathrm{sin}\theta _f,`$
$`\alpha _{\stackrel{~}{f}_2}^{(i)}`$ $`=`$ $`A_f^{(i)}\mathrm{sin}\theta _f+f_fv_a^{(i)}\mathrm{cos}\theta _f,`$
$`\beta _{\stackrel{~}{f}_1}^{(i)}`$ $`=`$ $`B_f^{(i)}\mathrm{sin}\theta _f+f_fv_a^{(i)}\mathrm{cos}\theta _f,`$
$`\beta _{\stackrel{~}{f}_2}^{(i)}`$ $`=`$ $`B_f^{(i)}\mathrm{cos}\theta _f+f_fv_a^{(i)}\mathrm{sin}\theta _f,`$
$`A_u^{(i)}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(gv_3^{(i)}+{\displaystyle \frac{g^{}}{3}}v_4^{(i)}\right),`$
$`A_d^{(i)}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(gv_3^{(i)}+{\displaystyle \frac{g^{}}{3}}v_4^{(i)}\right),`$
$`A_l^{(i)}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(gv_3^{(i)}+g^{}v_4^{(i)}\right),`$
$`B_u^{(i)}`$ $`=`$ $`{\displaystyle \frac{4}{3\sqrt{2}}}g^{}v_4^{(i)},`$
$`B_d^{(i)}`$ $`=`$ $`{\displaystyle \frac{2}{3\sqrt{2}}}g^{}v_4^{(i)},`$
$`B_l^{(i)}`$ $`=`$ $`\sqrt{2}g^{}v_4^{(i)}.`$
For the neutralino mixing element $`v_a^{(i)}`$, $`a=1`$ if $`f`$ is an up type quark while $`a=2`$ for leptons and down type quarks.

The $`g`$ and $`g^{}`$ are the $`SU(2)`$ and $`U(1)_Y`$ gauge couplings, while the $`v_i^{(j)}`$ are neutralino mixing elements, as defined by Tata. The $`f_f`$ is the Yukawa coupling for fermion $`f`$, and $`\theta _f`$ is the sfermion mixing angle. The integrals $`I`$, $`I^2`$, $`J`$ and $`K`$ are defined in Haber and Wyler.

Also,

$$g_{efgh}=\frac{e}{8\pi ^2}\left(m_{\stackrel{~}{Z}_i}C_i^{efgh}(JK)m_{\stackrel{~}{Z}_j}C_j^{efgh}(I^2JK)2m_{\stackrel{~}{W}_k}C_0^{efgh}J\right),$$
(9)
where

$`C_i^{efgh}`$ $`=`$ $`(i)^{\theta _i}(i)^{\theta _j}\left(D_k^{(i)}D_k^{(j)}E_k^{(i)}E_k^{(j)}\right),`$
$`C_j^{efgh}`$ $`=`$ $`(i)^{\theta _i}(i)^{\theta _j}\left(D_k^{(i)}D_k^{(j)}E_k^{(i)}E_k^{(j)}\right),`$
$`C_0^{efgh}`$ $`=`$ $`(i)^{\theta _i+\theta _j}\left(E_k^{(i)}D_k^{(j)}D_k^{(i)}E_k^{(j)}\right),`$
$`D_1^{(i)}`$ $`=`$ $`g(1)^{\theta _{\stackrel{~}{W}_1}}\left({\displaystyle \frac{\mathrm{cos}\gamma _R}{\sqrt{2}}}v_1^{(i)}+\mathrm{sin}\gamma _Rv_3^{(i)}\right),`$
$`D_2^{(i)}`$ $`=`$ $`g(1)^{\theta _{\stackrel{~}{W}_2}}\theta _y\left({\displaystyle \frac{\mathrm{sin}\gamma _R}{\sqrt{2}}}v_1^{(i)}+\mathrm{cos}\gamma _Rv_3^{(i)}\right),`$
$`E_1^{(i)}`$ $`=`$ $`g\left({\displaystyle \frac{\mathrm{cos}\gamma _L}{\sqrt{2}}}v_2^{(i)}+\mathrm{sin}\gamma _Lv_3^{(i)}\right),`$
$`E_2^{(i)}`$ $`=`$ $`g\theta _x\left({\displaystyle \frac{\mathrm{sin}\gamma _L}{\sqrt{2}}}v_2^{(i)}+\mathrm{cos}\gamma _Lv_3^{(i)}\right).`$
Furthermore,

$$g_{ijkl}=\frac{e}{(4\pi )^2}\left(m_{\stackrel{~}{Z}_i}C_i^{ijkl}K+m_{\stackrel{~}{Z}_j}C_j^{ijkl}(I^2K)+m_{\stackrel{~}{W}_k}C_0^{ijkl}I\right).$$
(10)
$`C_i^{ijkl}`$ $`=`$ $`(i)^{\theta _i}(i)^{\theta _j}\left(\mathrm{sin}^2\beta A_{k2}^{(i)}A_{k2}^{(j)}\mathrm{cos}^2\beta A_{k1}^{(i)}A_{k1}^{(j)}\right),`$
$`C_j^{ijkl}`$ $`=`$ $`(i)^{\theta _i}(i)^{\theta _j}\left(\mathrm{sin}^2\beta A_{k2}^{(i)}A_{k2}^{(j)}\mathrm{cos}^2\beta A_{k1}^{(i)}A_{k1}^{(j)}\right),`$
$`C_0^{ijkl}`$ $`=`$ $`(i)^{\theta _i+\theta _j}(1)^{\theta _{\stackrel{~}{W}_k}}\mathrm{sin}\beta \mathrm{cos}\beta \left(A_{k1}^{(i)}A_{k2}^{(j)}A_{k2}^{(i)}A_{k1}^{(j)}\right),`$
$`A_{11}^{(i)}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(gv_3^{(i)}+g^{}v_4^{(i)}\right)\mathrm{cos}\gamma _Rgv_1^{(i)}\mathrm{sin}\gamma _R,`$
$`A_{21}^{(i)}`$ $`=`$ $`\theta _y\left({\displaystyle \frac{1}{\sqrt{2}}}\left(gv_3^{(i)}+g^{}v_4^{(i)}\right)\mathrm{sin}\gamma _Rgv_1^{(i)}\mathrm{cos}\gamma _R\right),`$
$`A_{12}^{(i)}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(gv_3^{(i)}+g^{}v_4^{(i)}\right)\mathrm{cos}\gamma _L+gv_2^{(i)}\mathrm{sin}\gamma _L,`$
$`A_{22}^{(i)}`$ $`=`$ $`\theta _x\left({\displaystyle \frac{1}{\sqrt{2}}}\left(gv_3^{(i)}+g^{}v_4^{(i)}\right)\mathrm{sin}\gamma _L+gv_2^{(i)}\mathrm{cos}\gamma _L\right).`$
Finally,

$$g_{mnop}=\frac{e}{(4\pi )^2}\left(m_{\stackrel{~}{Z}_i}C_i^{mnop}K+m_{\stackrel{~}{Z}_j}C_j^{mnop}(I^2K)+m_{\stackrel{~}{W}_k}C_0^{mnop}I\right).$$
(11)
Here $`C_I^{mnop}`$ ($`I=i,j,0`$) are obtained from $`C_I^{ijkl}`$ by making the replacement $`\mathrm{sin}\beta \mathrm{cos}\beta ,\mathrm{cos}\beta \mathrm{sin}\beta `$. In the above, $`\gamma _L`$ and $`\gamma _R`$ are chargino mixing angles as defined by Tata; in that work, the $`\theta _x`$, $`\theta _y`$, $`\theta _i`$ and $`\theta _{\stackrel{~}{W}_i}`$ are also defined.

## Acknowledgments

We thank Xerxes Tata for discussions. This research was supported in part by the U.S. Department of Energy under contract number DE-FG02-97ER41022.


text:
# Multiple solutions of coupled-cluster equations for PPP model of [10]annulene

## I Introduction

The coupled-cluster (CC) method cizek:66  is one of the basic tools for studying the electronic-correlation effects in many-electron systems. Due to their nonlinearity, the CC equations may have multiple solutions; this problem was first analyzed by Živković and Monkhorst zivkovic:78 . Recently, Kowalski and Jankowski kowalski:98a  applied to the CC equations a new powerful mathematical tool—the homotopy (continuation) method drexler:78 . By using this tool they were able to find the complete sets of solutions for the CCD and the CCSD equations corresponding to some four-electron four-orbital systems, known as the H4 and P4 models. Very recently, Podeszwa et al. podeszwa:02b  used the homotopy method in a study of multiple solutions to the CCD equations for the Pariser-Parr-Pople (PPP) model pariser:53a ; pople:53  of benzene (6 electrons, 6 orbitals). The available data on the CC multiple-solution problem is still very limited. The purpose of this paper is to get some insight into the multiple solutions of the CCD equations (as well as the equations of some modified CCD methods) for the PPP model of annulene (10 electrons, 10 orbitals), a system which is already too big to be treated by the homotopy method.

Benzene and annulene are the first two members of the family of cyclic polyenes (annulenes) $`\mathrm{C}_M\mathrm{H}_M`$, where $`M=4m_0+2`$, $`m_0=1,2,\mathrm{}`$, and the assumed point-symmetry group is DMh. The problem of electronic correlations in annulenes, which may serve as prototypes of one-dimensional metal, have been thoroughly studied within the PPP model by Paldus and co-workers paldus:84a ; paldus:84b ; paldus:84c . Due to the high symmetry of annulenes, the contributions from the single excitations are absent in the PPP model, and the CCSD method reduces to the CCD one. Paldus et al. discovered and documented the breakdown of the CCD method for these systems. They also showed that this breakdown may be attributed to the neglect in the CCD method the terms providing the coupling between the D and Q excitations. Paldus et al. devised an approximate coupled-pair method, called ACPQ paldus:84a ; paldus:84b , in which these coupling terms (approximately) canceled certain quadratic terms in the CCD equations. Unlike the CCD method, the ACPQ method, and a related ACP-D45 method (ACP in short) introduced earlier jankowski:80 , were shown to perform well for the annulenes, being convergent and giving the correlation energies close to the full configuration-interaction (FCI) results. However, a recent CC study podeszwa:02a  of annulenes, taking into account the double (D), triple (T), and quadruple (Q) excitations in the CC operator, showed that even the CCDTQ method breaks down for these systems when the correlation effects become sufficiently strong. Moreover, the ACP and ACPQ methods were shown to fail for very large annulenes podeszwa:02c . Thus, the annulenes remain a challenging problem for the many-electron theories.

## II Finding multiple solutions of CCD, ACP, and ACPQ equations for annulene

In the CCD method, the CC operator corresponding to the double excitations depends on some linear parameters, hereafter referred to as the $`t_2`$ amplitudes. The CCD equations comprise a set of quadratic equations for the unknown $`t_2`$ amplitudes. In the nonorthogonally spin-adapted formalism (see, e.g., Ref. stolarczyk:84, ), the number of $`t_2`$ amplitudes is equal to $`K=n_\mathrm{o}n_\mathrm{u}(n_\mathrm{o}n_\mathrm{u}+1)/2`$, where $`n_\mathrm{o}`$ and $`n_\mathrm{u}`$ is the number of occupied and unoccupied orbitals, respectively (no symmetry other then the spin one is assumed). According to the Bézout theorem, such a set of quadratic equations may have up to $`2^K`$ solutions, a truly astronomical number in most applications. Therefore, finding the complete set of CCD solutions may be attempted only when the number of $`t_2`$ amplitudes is very small (e.g., due to the spatial symmetry of the many-electron system). In the homotopy study of the H4 model kowalski:98a  (corresponding to 6 spin- and spatial-symmetry–adapted $`t_2`$ amplitudes), the complete set of CCD solutions numbered 12, which is distinctly smaller than the Bézout limit (64). On the other hand, the application of the homotopy method to the PPP model of benzene podeszwa:02b  (having 8 $`t_2`$ amplitudes adapted to spin, spatial, and time-reversal symmetry) has brought in some quite disturbing findings: not only the total number of solutions (230) approached the theoretical limit (256), but some solutions led to the correlation energies much more negative than that of the “standard” CCD solution (which approximate the ground-state FCI solution); in addition all the “nonstandard” solutions were hardly similar to the FCI ones.

In the PPP model of annulene, the CCD equations adapted to spin, spatial, and time-reversal symmetry correspond to 29 $`t_2`$ amplitudes. The ACP and ACPQ equations are obtained by neglecting (ACP and ACPQ) and modifying (ACPQ) some of the quadratic terms in the CCD equations, see, e.g., Ref. podeszwa:02c, . In the present paper we use the PPP model of annulenes described in Ref. podeszwa:02a,  and follow the notation used there: $`\beta `$ $`(0)`$ eV is the resonance integral of the PPP model, and $`\gamma (R)`$ is a potential describing the long-range electrostatic interactions. We shall report the results obtained for two variants of the PPP model: the PPP-MN model, employing the Mataga-Nishimoto potential mataga:57  ($`\gamma (R)=e^2[R+e^2(\gamma ^0)^1]^1`$), with $`\gamma ^0=10.84`$ eV, and the PPP-P model, using Pople’s potential pople:53  ($`\gamma (R)=e^2/R`$ for $`R0`$), with $`\gamma ^0=14.25878`$ eV ($`e`$ is the elementary charge).

If the CCD results for benzene podeszwa:02b  may be extrapolated, the CCD (and ACP and ACPQ) equations for annulene may have more than $`10^8`$ solutions; that makes the application of the homotopy method hopeless in this case. On the other hand, it was found in Ref. podeszwa:02b,  that some of the real CCD solutions for benzene were stable within their close neighborhoods: such a solution (subject to the numerical rounding errors) led to a convergent iteration process when the direct-inversion-of-the-iterative-subspace (DIIS) method pulay:82  was used. Thus, one may expect that the DIIS method is able to trace various solutions, if only the starting points are suitably chosen.

A standard method of solving the CCD equations is based on the iterative procedure represented in Eq. (14) of Ref. podeszwa:02b, . By starting from the $`t_2`$ amplitudes equal to zero, one recovers in the first iteration the $`t_2^{(1)}`$ amplitudes corresponding to the first order of the Møller-Plesset (MP) perturbation theory; in the subsequent iterations higher-order MP contributions are added. Such an iteration procedure, if convergent, furnishes the standard CCD solution corresponding (as it is often a priori assumed) to the ground state of the system. The CCD electronic-correlation energy is then calculated as a certain linear function of the converged $`t_2`$ amplitudes (the $`t_2^{(1)}`$ amplitudes correspond to the second-order MP correlation energy and will be, after Ref. podeszwa:02a, , referred to as the MP2 amplitudes). The above treatment applies also to the ACP and ACPQ equations. In annulenes one may also use a kind of analytical-continuation procedure to find the standard solution for those values of parameter $`\beta `$ for which the standard procedure (starting from the MP2 amplitudes) does not converge: one uses as the starting point the converged $`t_2`$ amplitudes found for the sufficiently close $`\beta `$ value.

The multiple solutions of the CCD, ACP, and ACPQ equations presented in this paper were found in a quasi-random manner, by generating a pool of the initial $`t_2`$ amplitudes and then applying the DIIS procedure. The starting pool included the results of several hundred converged CCD, ACP, and ACPQ calculations obtained for different $`\beta `$ values and for different forms of function $`\gamma (R)`$ (these sets of $`t_2`$ amplitudes indeed looked random). Most of the calculations starting from this quasi-random amplitudes either diverge, or converge to the standard solution, but in some cases a new solution is obtained. Certainly, this is not a very efficient method of finding the multiple solutions, but since a single CCD calculation for annulene usually takes only several seconds, a lot of data may be generated that way. We have found that the strongly correlated regime ($`\beta >1.5`$ eV) is rich in multiple solutions, while in the weakly correlated regime ($`\beta <3.0`$ eV) practically only the standard solution could be reached that way. Therefore, we performed calculations for several values of $`\beta `$ from interval $`(1.0\text{ eV},0\text{ eV})`$. These results were then extended as far as possible by a careful application of our analytical-continuation procedure. The obtained results are not supposed to provide the complete set of real solutions for annulene. Nevertheless, they provide a rather intriguing glimpse at the “tip of the iceberg” of the multiple solution problem for this system. Let us note that the multiple solutions of the ACP and ACPQ equations have never been studied before.

## III Results

In Figs. 13 we present the electronic-correlation energies corresponding to the multiple solutions of the CCD, ACP, and ACPQ equations, respectively. These results are obtained within the PPP-MN model for $`5.0\text{ eV}\beta <0\text{ eV}`$, each solution is represented by a continuous line. The line endpoints (not touching the graph’s boundaries) indicate the limiting $`\beta `$ values corresponding to our analytical-continuation procedure. It is seen that some of the new solutions exist for a very broad range of $`\beta `$, while some can be continued only within a very narrow interval. The ground-state FCI correlation energies, shown as the broken line, are provided for a reference. The solution line that approaches the FCI line in the weakly correlated regime corresponds to the standard solution discussed in the previous papers paldus:84a ; podeszwa:02a .

The CCD results in Fig. 1 are rather simple: The standard-solution line is the only extending for all the $`\beta `$ values. (When the standard iteration procedure is combined with the DIIS technique, the standard CCD solution for annulene can be obtained for all negative values of $`\beta `$, as shown previously by Paldus et al. paldus:84a  by using the Newton-Raphson method.) In the strongly correlated regime ($`\beta >1.5`$ eV), the standard CCD solution provides a rather poor approximation to the ground state, and the corresponding correlation energy falls well below the FCI line. Other solutions have the energies higher than the FCI result, some of them may represent excited states. Most of them can be found only in the proximity of $`\beta =0`$, where the quasidegeneracy effects are the strongest.

The ACP and ACPQ results (Figs. 2 and 3 respectively) display a different and more complicated behavior. Unlike the CCD standard solution, the ACP and ACPQ standard solutions stay close to the FCI line also in the strongly correlated regime (although the corresponding $`t_2`$ amplitudes show increasing deviations from the FCI $`t_2`$ amplitudes as $`\beta 0`$ podeszwa:02a ). The ACP and ACPQ standard solutions can be found for $`\beta <0`$ by applying the standard iteration procedure starting from the MP2 amplitudes podeszwa:02a . Quite a few nonstandard solutions have been found, some of them extending for the whole range of $`\beta `$. Quite unusual is the presence of the solution lines lying below the standard-solution line (the “underground” lines). Such exotic solutions have been observed in the study of benzene podeszwa:02b ; however, they were unstable in the DIIS iterative procedure (could be found only by the homotopy method). The “underground” solutions obviously have no physical meaning, and it is a little disturbing that they may appear while iteratively solving the ACP or ACPQ equations. One of the ACPQ “underground-solution” lines is especially peculiar: it starts with the lowest energy for $`\beta =0`$ and rises quite steeply with the decreasing $`\beta `$, crossing several solution lines, including the standard-solution line (apparently there is no “non-crossing rule” for different solution lines).

Quite intriguing are the solutions of the CCD equations corresponding to the PPP-P model, displayed in Fig. 4 together with the ground-state FCI results. We were able to find only two solutions, hereafter referred to as CCD1 (the standard solution) and CCD2 (the nonstandard one). Despite the use of the extensive quasi-random search, the CCD iterations invariably led to CCD1 or CCD2 (these solutions were so strong attractors that the iterations rarely diverged). In contrast to the CCD results for the PPP-MN model, the iterations starting directly from the MP2 amplitudes converge to the CCD1 solution in the weakly correlated regime, while in the strongly correlated regime the CCD2 solution is obtained. In Fig. 4 we mark by the squares the endpoints of the respective direct-convergence domains. One may thus obtain the CCD1 and CCD2 solutions in a deterministic way, by starting the calculations within these regions. Moreover, by applying our analytical-continuation procedure, we are able to extend a given solution line outside its domain of (direct) convergence.

The CCD1-solution line resembles closely the line corresponding to the standard CCD solution for the PPP-MN model, see Fig. 1, while the CCD2 line has no counterpart in that figure. The CCD1 and CCD2 lines cross for $`\beta =1.17479`$ eV, in the vicinity of the FCI-solution line. In Table 1 we present the CCD1, CCD2, MP2 and FCI correlation energies calculated for three representative $`\beta `$ values, including the intersection coordinate. In addition, we show there the values of parameters $`\theta `$ and $`\eta `$ which provide a measure of similarity between the vector of $`t_2`$ amplitudes corresponding to a given solution and the vector of $`t_2`$ amplitudes extracted from the linear coefficients of the ground-state FCI wavefunction: $`\theta `$ is the angle between the vectors, and $`\eta `$ is the ratio of the vector lengths (for details, see Ref. podeszwa:02a, ).

From Fig. 4 and Table 1 a consistent picture emerges: The CCD1 solution provides an acceptable approximation to the ground-state FCI solution up to the vicinity of the CCD1/CCD2 crossing; then, in the strongly-correlated region, it deteriorates and hardly deserves the name of the standard solution. The CCD2 solution, on the other hand, starts as the “underground” in the weak correlation regime (bearing no resemblance to the FCI solution), but it improves in the strongly-correlated regime, and, finally, becomes a slightly better (but still poor) approximation to the FCI solution (the CCD2 correlation energies are surprisingly similar to the FCI ones in the strongly correlated regime). Thus, it looks as CCD1 and CCD2 swap their roles as the standard and “underground” solutions. The behavior of the MP2 results agrees with the observed convergence pattern for iterations starting from the MP2 amplitudes. Interestingly, if one performed a CCD study of the PPP-P model of annulene by applying a standard iteration procedure (starting with the MP2 amplitudes) at some representative $`\beta `$ values (say, $`5.0`$, $`2.5`$, $`1.5`$, and $`0.5`$ eV), the conspicuous switch between the CCD1 and CCD2 solutions would have been passed unnoticed. By looking at the CCD correlation energies alone, one would then proclaimed a very good performance of the CCD method up to very small $`\beta `$ values, contrary to the results found by Paldus et al. paldus:84a  for the PPP-MN model.

The study of the ACP and ACPQ equations in the PPP-P model brings in some additional unexpected findings: the standard ACP solution cannot be found for $`\beta `$ between $`0.80`$ and $`0.4`$ eV, while the standard ACPQ solution encounters problems already for $`\beta >2.53`$ eV. Apparently, in the strongly correlated regime of the PPP-P model both methods become unstable.

Our search of the multiple real solutions of the CCD, ACP, and ACPQ equations has brought several unexpected findings: (i) while solving the CC equations by iterations, one may arrive at an “underground” solution, (ii) in the PPP model, the character of multiple solutions is sensitive to the form of the long-range potential (the $`\gamma `$ function), (iii) different solutions may play the role of the standard solution when the strength of the electronic-correlation effects is varied (this may apply to the changes in molecular geometry, e.g., to the bond breaking). In conclusion: it seems that the problem of multiple solutions of the CC equations deserves still more attention. Some variants of our random-search approach may be tried also on larger many-electron systems.

This work was supported in part by the Committee for Scientific Research (KBN) through grant No. 7 T09A 019 20.