text:
# Statistics of Lyapunov exponent in one-dimensional layered systems

## I Introduction

The fact that the electronic localization in disordered systems is of wave nature has led to suggestion that classical waves could be similarly localized in random systems. The effort in searching for localization of classical waves such as acoustic and electro-magnetic waves is tremendous. It has drawn intensive attentions from both theorists  and experimentalists. Snice the pioneering work of Anderson et.al., the concepts of universality and scaling have become important. These ideas stem from the insensitivity of the macroscopic laws to microscopic details; that is, systems or models which differ from each other on a microscopic level can show identical macroscopic behavior. According to the hypothesis of single-parameter scaling (SPS), if the localization behaviors of a 1D disordered system obey SPS, the Lyapunov exponent or inverse localization that characterizes the degree of localization will be proportional to its variance.

In a recent study of acoustic (AC) waves propagation in a one dimensional randomly layered system we have found that the statistics relation between Lyapunov exponent (LE) $`\gamma `$ and its variance $`\text{var}(\gamma )`$ do not follow the preditions of single parameter scaling (SPS). However, in an earlier study on 1D localization behaviors of electromagnetic (EM) waves, the author claimed that the non-universal bahaviors of LE will disappear and SPS will be restorted while the randomness of the system exceeds a critical value. Since for 1D propagation the AC and EM waves are in fact mathematically equivalent, i.e., there exists one to one correspondence between these two kinds of waves, it looks quite impossible that they can have different localization behaviors. In order to understand where do the main differences between these two kinds of models come from, in this paper we study EM and AC systems simultaneously. Two EM models and one AC model are studied in this paper. We find that though the statements made by Deych et.al.  are correct in their chosen case, however, the applicability of SPS is more or less based on the fact that the impedance contrast between the constituents of the wave media is closed to 1. Without this restriction then even in the EM systems SPS will not be restored in the high randomness limit.

This paper is organized as follows. In the next section we explain the correspondence between EM and AC waves and define the three models employed in this paper. In Sec. III we we first review the previous results and then discuss the numerical results of the three chosen models. A possible explanation of the origin of the novel properties of the AC localization are also proposed. Concluding remarks are given in section IV.

## II Theory and models

To begin with, we explain the one to one correspondence of the 1D propagation between AC and EM waves. For simplicity while without destroying the generality we restrict our discussion to the monochromatic waves with time dependence $`e^{i\omega t}`$. We also assume that waves are nomally incident on the left boundary of the media and propagate along the x-axis. The 1D propagation of AC waves under these assumptions is govened by

$$\frac{d}{dx}\left(\frac{1}{\rho }\frac{dp}{dx}\right)=\frac{\omega ^2p}{\rho c^2}$$
(1)
where $`\rho =\rho (x)`$, $`p=p(x)`$ and $`c=c(x)`$ represent the mass density, the pressure, and the phase velocity of the wave in the media, respectively. In the special case of layered media considered in , $`\rho `$ and $`c`$ are all constants in a single layer. Across an interface that separate two layers either $`\rho `$ or $`c`$ jump to a diferent value but pressure $`p`$ and media vibration velocity

$$u=\frac{1}{i\omega \rho }\frac{dp}{dx}$$
(2)
must be continuously connected. The continuity conditions for every interfaces together with the wave equation (1) itself determine the dynamics of the whole system.

Similarly, the equation governing 1D propagation of monochromatic EM waves can be deduced from Maxwell equation and is written as

$$\frac{d}{dx}\left(\frac{1}{\mu }\frac{dE}{dx}\right)=\frac{\omega ^2ϵE}{c_0^2}=\frac{\omega ^2E}{\mu c^2}$$
(3)
or

$$\frac{d}{dx}\left(\frac{1}{ϵ}\frac{dH}{dx}\right)=\frac{\omega ^2\mu H}{c_0^2}=\frac{\omega ^2H}{ϵc^2},$$
(4)
where $`E=E_y`$ and $`H=H_z`$ are the electric and magnetic fields , $`ϵ`$ and $`\mu `$ stand for the permittivity and permeability, $`c_0`$ is the speed of light in vacuum and $`c=c_0/\sqrt{ϵ\mu }`$ is the speed of light in the media. Comparing Eq. (3) and (4) with Eq. (1) one can easily recognize the equivalence between AC and EM models via the substitutions

$$Ep,Hu,ϵ1/\rho c^2,\mu \rho $$
(5)
in Eq. (3) or

$$Hp,Eu,ϵ\rho ,\mu 1/\rho c^2$$
(6)
in Eq. (4).

In spite of these similarities, the two models discussed in  and  indeed have some different features. First, in the optical case that described in Eq. (3), one usually assumes $`\mu =1`$. However, in our AC model the corresponding quantity is $`\rho `$. The mass density ratio btween water and air is about 775, a very large value. Second, the phase velocity ratio in our model is 4.455, whereas in  the ratio is 1.0945. Third, the thickness ratio in our model is 9999/1, much larger than 1/1 that considered in . Fourth, we randomize the thickness of the water layers (medium with high phase velocity) and keep the thickness of air layers (medium with low phase velocity) constant. In  they randomize the thickness of the layers with low phase velocity.

To clarify where do the major differeces come from, we define three models and study them numerically. The system for each model is a composite made of two kinds of material $`A`$ and $`B`$ with corresponding thickness $`a_j`$ and $`b_j`$ in the $`j`$th $`A/B`$ layer. For simplicity hereafter we assume $`a_j=a`$, i.e., all $`A`$ type layers have the same size. Any quantity $`Q`$ in $`A`$ and $`B`$ type layers are denoted as $`Q_a`$ and $`Q_b`$ respectively.

* Model 1 is an optical model with $`ϵ_b/ϵ_a=2`$ and $`\mu _b/\mu _a=1`$. The thickness ratio is $`b/a=1`$, where $`bb[1\mathrm{\Delta },1+\mathrm{\Delta }]`$ and $`\mathrm{\Delta }[0,1]`$.
* Model 2 is also an optical model with $`ϵ_b/ϵ_a=20`$ and $`\mu _b/\mu _a=1`$. The thickness ratio is $`b/a=1`$, where $`bb[1\mathrm{\Delta },1+\mathrm{\Delta }]`$ and $`\mathrm{\Delta }[0,1]`$.
* Model 3 is our previously considered acoustic model with $`c_b/c_a=4.455`$ and $`\rho _b/\rho _a=755.2`$. The thickness ratio is given by $`b/a=9999`$, where $`bb[1\mathrm{\Delta },1+\mathrm{\Delta }]`$ and $`\mathrm{\Delta }[0,1]`$.

In a single layer the impedance of AC waves is given by $`\rho c`$ and the impedance of EM waves is given by $`ϵc`$ or $`\mu c`$. Thus we see that in Model 3 the impedance contrast is about $`3365`$ and impedance contrast in Model 1 and 2 are about $`1.414`$ and $`4.472`$ respectively.

## III Numerical results

Before the discussion of the numerical results we first summarize the relevant results of our previous study. There both LE and its variance as functions of frequency were studied. At low disorders, the variance of LE inside the gaps is small. Contrast to the optical case, there are no double maxima inside the gap. With increasing disorder, double peaks appears inside the allow bands. When exceeding a certain critical value, however, the double peaks emerge. The higher frequency, the lower is the critical value. The increasing disorder reduces the band gap effect and smears LE. We also ploted LE-variance relations. However, with increasing disorder, we did not observe linear dependence between LE and its variance, as expected from the single parameter scaling theory.

Now we turn to the discussion of the three chosen models. When randomness $`\mathrm{\Delta }=0`$, the layered systems become periodic. Eigenfunctions of wave equation in a periodic environment are Bloch waves. As is well known, band structure appears in this situation. The understanding of band structure is very important and helpful in the following discussions. The dispersion relation of the Bloch waves in the underlying periodic system is given by

$$\mathrm{cos}Kd=\mathrm{cos}k_aa\mathrm{cos}k_bb\mathrm{cosh}2\eta \mathrm{sin}k_aa\mathrm{sin}k_bb$$
(7)
with $`d=a+b`$ representing the thickness of the space period and $`K`$ the Bloch wave number. Here for model 1 and 2 function $`\mathrm{cosh}2\eta `$ is defined as

$$\mathrm{cosh}2\eta =\frac{1}{2}\left(\sqrt{\frac{ϵ_b}{ϵ_a}}+\sqrt{\frac{ϵ_a}{ϵ_b}}\right)$$
(8)
and for model 3 it is given by

$$\mathrm{cosh}2\eta =\frac{1}{2}\left(gh+\frac{1}{gh}\right),g=\rho _a/\rho _b,h=c_a/c_b.$$
(9)
In the frequency ranges that $`Kd`$ are real, the waves are freely propagating in the media and by definition the frequency ranges correspond to the pass bands. Beyond the pass bands $`Kd`$ are not purely real numbers, solutions of wave equation that satisfy appropriate boundary conditions in both the left and right infinity do not exist, and thus the frequency ranges are refer to as the band gaps. If the media is semi-infinite and the periodicity is ended by a boundary, say, the left boundary, then the waves are localized in the vicinity of the boundary. The penetration depth is equal to $`1/|\mathrm{Im}(K)|`$. Fig. 1(a1)-(c1) plot the band structures of model 1-3. Solid curves represent $`\mathrm{Re}(Kd)`$ and cover the pass bands. Broken lines cover the band gaps and represent the inverse penetration depths $`\mathrm{Im}(Kd)`$. Model 1 has very wide pass bands and very narrow band gaps. Band gaps are wider than pass bands in model 2. Model 3 has very wide band gaps and very narrow pass bands. An important feature of model 3 is that in the band gaps the the penetration depth is very small. Even in the first gap (which has the longest penetration depth) the penetration depth is smaller than one period.

To see how randomness influences the transmission properties we select the frequency range that around the second gap and plot the transmission curves on (a2)-(c2). In this calculation 1600 layers (800 periods) are used in model 1 and 200 layers (100 periods) are used in model 2 and 3. In the gaps the transmission rate is almost zero. When randomness is increased we observe that the transmission rate in the pass bands have been reduced much in the model 2 and model 3. On the other hand, the band gaps in model 2 and 3 seem to be more robust than in model 1.

To further explore the influence of randomness we study both LE and its variance. The results are shown in Fig. 2. (A1)-(A5) are results of model 1 and (B1)-(B5) are results of model 2. LE is denoted as $`\gamma `$ and its variance is denoted as $`\text{var}(\gamma )`$. Here $`\gamma `$ and $`\text{var}(\gamma )`$ are defined as

$$\gamma =\underset{N\mathrm{}}{lim}\gamma _N$$
(10)
with

$$\gamma _N=\frac{1}{2N}\mathrm{ln}\left(\frac{1}{T_N}\right)$$
(11)
and

$$\mathrm{var}(\gamma )=\underset{N\mathrm{}}{lim}(\gamma _N^2\gamma _N^2),$$
(12)
here $`T_N`$ is the transmission rate for system with $`2N`$ layers ($`N`$ periods in the corresponding periodic system), and notation $`\mathrm{}`$ represents the ensemble average. The sample size is chosen in such a way that it is much larger than the localization length and the ensemble average is carried out over 200 random configurations. As expected, when the randomness is small the LE can be approximated by the inverse penertration depth for wave propagating in the underlying semi-infinite periodic system. Like in the case discussed in Ref., the double peaks of $`\text{var}(\gamma )`$ first appear near the gap edges. The peaks of $`\text{var}(\gamma )`$ imply the fluctuations of transmission. Further increasing $`\mathrm{\Delta }`$, the $`\text{var}(\gamma )`$ peaks become fatter and flatter. We observe that if in the vicinity of a peak there is another peak, then increasing the randomness will cause them to merge. For model 1, as $`\mathrm{\Delta }`$ is increased, a pair of $`\text{var}(\gamma )`$ peaks merge with each other inside the gap, following the scenario of Ref. However, for model 2 when the randomness is increased the pairs of $`\text{var}(\gamma )`$ peaks tend to merge with each other in the pass bands and finally destroy the pass bands. It seems that in model 2 the merging of the double peaks of $`\text{var}(\gamma )`$ in a pass band will never be completed. The merging tendency of a pair of double peaks merely increase the LE and destroy the pass band they belong to.

We also plot the LE versus its variance in Fig. 3. There we observe that although for model 1 (with small dielectric contrast between two layers like that studied in ) SPS seems still a good approximation in large disorder limit, deviation from SPS is clearly observed in model 2. Similar quantities have also been calculated for model 3 and ploted in Fig. 4 and Fig. 5. There we also observe novel behaviors of $`\gamma `$ and $`\text{var}(\gamma )`$ and find even larger deviation from SPS as reported in .

From these observations, we find:

* How does the $`\text{var}(\gamma )`$ vary with the randomness depends on many parmeters. For example, the impedance constrast, the thickness ratio and which kind of layers (for example, high phase velocity, low phase velocity etc.) is randomized.
* When impedance constrast is large, it seems that the deviation from SPS is a usual feature. The results reported in  is based more or less on the fact that the impedance contrast between two layers is closed to 1.

In order to understand why the large deviation from SPS in model 3 is established we study its transmission properties in more detail. In Fig. 6(a) we plot the transmission rate of the AC waves through a pair of air blocks of thickness $`a`$. The two air blocks are seperated by a water layer of thickness $`b=9999a`$. Similarly in Fig. 6(b) we plot the transmission for the system with 100 air blocks. Comparing these two diagrams we find that the air blocks are very strong scatterers and thus as few as only two air blocks are enough to determine the ranges of band gaps and pass bands. If the localization effect of AC waves in model 3 is mainly determined by the multiple sacttering of AC waves between pairs of air blocks, then one would expect that for a system with 100 air blocks the transmission $`T_{N=100}`$ in the band gaps can be approximated by $`T_{N=2}^{50}`$, here $`T_{N=2}`$ refers to the transmission rate for two air blocks illustrated in 6(a). We indeed observed this result in Fig. 6(c). When randomness is not very large, this local effect explains why the phase averaging process in model 3 is so inefficient that the randomness cannot modify the LE much in the band gaps.

## IV Concluding Remarks

In this paper we studied the statistics of localization properties in one-dimensional layered systems. Luapunov exponent and its variance are compared for three chosen models. We find that the band structures of the corresponding periodic systems influence the localization properties much if the impedance contrast between neighboring layers is not close to 1. In general the single parameter scaling is not very accurate and more model-dependent parametrs should be included in the detailed discriptions of localization behaviors.

## ACKNOWLEDGEMENT

The work received support from National Science Council (No. NSC89-2611-M008-002 and NSC89-2112-M008-008).


format: text

text:
# Missing Link Found? — The “runaway” path to supermassive black holes

## 1 Introduction

There is rapidly growing evidence for SMBHs in the centers of many galaxies (for a review see Kormendy and Richstone (1995)). There are too many examples to list here; indeed, there are only a few galaxies for which observations indicate that a central SMBH does not exist(Kormendy and McClure, 1993). Many authors have pointed out that the mass of the central black hole $`m_{BH}`$ correlates linearly with the mass of bulge $`M_b`$, i.e., the ratio of $`m_{BH}`$ to $`M_b`$ is almost constant (0.002 (Kormendy and Richstone, 1995) to 0.006(Magorrian et al., 1998)). This suggests that the formation of the central BH is somehow related to the formation of the bulge.

The formation mechanism of SMBHs is not well understood. Our theoretical understanding has not advanced much beyond the scenarios described by Rees(Rees, 1978, 1984) in the early 1980s. In the famous diagram by Rees, there were basically two paths from gas clouds to massive black holes. The first is direct monolithic collapse, the second is via the formation of a star cluster, with subsequent runaway collisions leading to BH formation. Previous numerical studies, however, have demonstrated that neither path is likely. In the first, a massive gas cloud is more likely to fragment into many small clumps, in which stars then form, so direct formation of a massive BH from a gas cloud seems unlikely. In the second, stellar dynamics in star clusters does not easily lead to the formation of massive black holes. A number of low-mass black holes (masses around 10M) are formed via the evolution of massive stars, and these black holes do indeed sink to the center of the cluster through dynamical friction and form binaries by three-body encounters. However, recent $`N`$-body simulations(Portegies Zwart and McMillan, 2000) have demonstrated that practically all of these black-hole binaries are ejected from the cluster by recoil following interactions with other black holes (or BH binaries) before they can merge through gravitational radiation.

## 2 IMBHs in M82

Matsumoto et al. (2001) have identified nine bright compact X-ray sources in the central region of M82 using recent Chandra data. The brightest source (No. 7 in their Table 1) had a luminosity of $`9\times 10^{40}\mathrm{ergs}^1`$ in Jan 2000, corresponding to a black hole with a minimum mass of $`700\mathrm{M}_{}`$ (assuming emission at the Eddington limit). It probably consists of a single compact object, as its X-ray flux shows rapid time variation(Matsumoto et al., 2001). This is the first detection of a BH with a mass much greater than $`100\mathrm{M}_{}`$ but much less than $`10^6\mathrm{M}_{}`$. Among the eight other sources, at least three (5, 8 and 9) have Eddington masses greater than $`30\mathrm{M}_{}`$.

Matsushita et al. (2000) observed the same region in using Nobeyama Millimeter Array and found a huge expanding shell of the molecular gas. They estimated the age and kinetic energy of the shell to be around 10 Myr and $`10^{55}`$ erg, which suggest a strong starburst took place a few Myrs ago.

Harashima et al. (2001) observed the same region in the IR (J, H, and K-band) using the CISCO instrument on the SUBARU telescope(Harashima et al., 2001). They identified a number of young compact star clusters, at least four of them coinciding with the X-ray sources within the position uncertainty of Chandra and SUBARU. The SUBARU field of view is smaller than that of Chandra, and three of the X-ray sources are outside of the IR field. Thus, four out of six bright X-ray sources in the IR field of view have been identified with bright star clusters. The logical conclusion from these observations is that most of Chandra X-ray sources, including the brightest one with an Eddington mass of $`700\mathrm{M}_{}`$, are formed in star clusters.

Therefore, we now have two important observational results. The first is that a BH with intermediate mass ($`100<M_{bh}/M_{}<10^6`$) has been found. The second is that it coincides with a young compact star cluster. In the following, we discuss how these findings change our understanding of the formation of supermassive BHs. We first discuss how IMBHs can be formed in young compact star clusters, then how IMBHs might grow into SMBHs.

## 3 IMBH formation through runaway growth

Here we consider some formation mechanisms for IMBHs in star clusters (see figure 1). An obvious way to form an IMBH is the collapse of a supermassive star. However, this possibility seems to have little observational or theoretical support. Observationally, no stars with masses larger than $`200\mathrm{M}_{}`$ are known. The most luminous star currently known is the Pistol Star(Figer et al., 1998), with an estimated present mass of less than $`200\mathrm{M}_{}`$. Recent radio observations of this star have detected signs of a strong stellar wind(Lang, et al., 1999), which is consistent with the theoretical prediction that such massive stars are dynamically unstable(Figer et al., 1998). We therefore regard it as unlikely that an IMBH with mass exceeding $`700\mathrm{M}_{}`$ could have formed by the collapse of a single massive star.

An alternative possibility is the formation and growth of massive stars (and IMBHs) through successive merging. Massive stars in star clusters have higher merger rates than less massive cluster members (or field stars), for the following three reasons. First, they have larger geometrical cross sections, and these are further enhanced by gravitational focusing. Second, because of mass segregation, these stars tend to be found in the high-density region near the cluster center. Third, the most massive stars at the center tend to form binaries through three-body interactions, and binary membership greatly increases their chances of merging. A significant fraction of binary-single star encounters lead to complex resonances during which very close encounters can occur. Many such encounters lead to physical collisions between stars (Hut and Inagaki, 1985; McMillan, 1986). If these effects are strong enough, we expect that a “merging instability”(Lee, 1987), or a runaway growth of the most massive star, will occur in the cluster core.

Lee (1987) performed Fokker-Planck simulations of globular clusters with merging of stars taken into account. He included the first two of the three effects mentioned above, and concluded that the merger instability is unlikely unless the mass of the globular cluster is very large ($`N>10^7`$, where $`N`$ is the number of stars). Recently, however, $`N`$-body simulations carried out by Portegies Zwart et al. (1999) have demonstrated that runaway merging can take place in much smaller systems containing only $`12,000`$ stars. The primary reason for this discrepancy is that, in $`N`$-body simulations, all three of the effects described above are automatically included, whereas in the Fokker-Plank calculation by Lee the third was not modeled. A second reason is that Portegies Zwart et al. considered a young, compact cluster, similar to R136 in the LMC or the Arches and Quintuplet systems in the galactic center, while Lee considered normal globular clusters. Thus, the dynamical timescale was much shorter in the models considered by Portegies Zwart et al., allowing merging to occur before stellar evolution eliminated the most massive stars.

Portegies Zwart et al. found that, in one case, the most massive star experienced more than ten collisions and reached a mass of around $`200\mathrm{M}_{}`$ before experiencing a supernova. There is considerable uncertainty as to how much mass would remain as a BH after the supernova explosion of such a massive star, but it is quite likely that the remnant black hole would still be one of the most massive objects in the cluster, and that the runaway merging process would continue. Although the geometrical cross section of a BH is small, “merging” would take place when a star approached within its tidal radius, leading to a relatively large merger cross section.

In order for runaway merging to occur, the dynamical friction timescale for the most massive stars must be short enough that they can sink to the center during their lifetime of several Myr. The dynamical friction timescale may be expressed asBinney, and Tremaine (1987):

$$t_{fric}=\frac{1.17}{\mathrm{log}\mathrm{\Lambda }}\frac{r^2v_c}{Gm}2.7\times 10^9\left(\frac{r}{10\mathrm{p}\mathrm{c}}\right)^2\left(\frac{r_h}{10\mathrm{p}\mathrm{c}}\right)^{1/2}\left(\frac{M}{10^6M_{}}\right)^{1/2}\left(\frac{20M_{}}{m}\right)\mathrm{yr},$$
(1)
where $`\mathrm{log}\mathrm{\Lambda }`$ is the Coulomb logarithm, $`G`$ is the gravitational constant, $`v_c`$ is the local velocity dispersion, $`r`$ is the distance from the center of the cluster, $`r_h`$ and $`M`$ are the half-mass radius and the total mass of the cluster, and $`m`$ is the mass of the star.

In the following, we consider how the dynamical friction works in the cluster found in M82. From the infrared luminosity, Harashima et al. (2001) estimate that the total mass of the cluster is $`5\times 10^6\mathrm{M}_{}`$. They also estimated the seeing-corrected radius of the cluster as 5 pc. The dynamical friction timescale is around 800 Myr for stars at the half mass radius of the cluster. This looks too long, but in fact it is not. As can be seen from equation (1), the timescale depends rather strongly on $`r`$. For $`r<r_h/10`$, a volume which still contains about 5% of the total cluster mass, the dynamical friction timescale is less than 10 Myr. We therefore conclude that there is time for a fair fraction of the most massive stars in the cluster to sink to the cluster center and undergo runaway merging before exploding as supernovae.

After the BH has become much more massive than other cluster members, it forms a cusp near the cluster center (Bahcall, and Wolf, 1976), and continues to consume other stars. Unfortunately, no realistic simulations of this phase of the evolution are available. Marchant and Shapiro (1980) performed Monte-Carlo simulations of this stage for a simplified cluster containing $`3\times 10^5`$ solar-mass stars and one 50 solar-mass seed BH. They found that the BH mass jumped to over $`10^3M_{}`$ (0.3% of the cluster mass) almost immediately after they put the BH into the system. After this initial rapid growth, a slower phase ensued, with a doubling timescale comparable to the relaxation time of the cluster. Their result should be regarded as a lower limit on the BH growth rate, since realistic effects, in particular the presence of a mass spectrum, would greatly enhance the accretion rate. Taking these effects into account, it seems safe (even conservative) to suppose that 0.1% of the total cluster mass accretes to form a $`5000M_{}`$ central BH in a few Myr.

As stated above, there are more than 10 bright star clusters in the vicinity of the IMBH host cluster in M82, some of them apparently hosting small BHs. Their age is around 10 MyrHarashima et al. (2001). Also, the starburst in M82 is a long-duration event, having started at least 200 Myr agode Grijs et al. (2001). If we assume that the clusters are formed in a constant rate, we would conclude that around 200 clusters have been formed. The star formation rate certainly has not been constant, so there is considerable uncertainty in our estimate. Even so, we believe it is safe to assume that around one hundred clusters similar to our host cluster have formed in total, and that a considerable fraction of them host IMBHs.

## 4 Building up the central SMBH

The growth rate of the IMBH in a star cluster slows once massive stars disappear (after $`100`$ Myr). What happens to the cluster and the IMBH after the initial rapid growth phase? In the discussion below, we assume a cluster mass of $`5\times 10^6`$ M and a half-mass radius of 5pc.

The cluster is subject to two evolutionary processes: evaporation through two-body relaxation and orbital decay (sinking) via dynamical friction. Evaporation is driven partly by thermal relaxation and partly by stellar mass loss. Portegies Zwart, et al. (2000) estimated that the evaporation timescale for a Roche-lobe-filling compact star cluster is around 2-3 half-mass relaxation times, which is of the order of a few Gyr for our star clusters. Rewriting equation (1) using appropriate scaling for this case, we find that timescale on which the cluster sinks to the galactic center via dynamical friction is

$$t_{fric}6\times 10^8\left(\frac{r}{1\mathrm{k}\mathrm{p}\mathrm{c}}\right)^2\left(\frac{v_c}{100\mathrm{k}\mathrm{m}\mathrm{s}^1}\right)\left(\frac{5\times 10^6M_{}}{m}\right).$$
(2)
Clusters initially within 1 kpc of the galactic center can therefore reach the center within one Gyr. Note that an IMBH can reach the galactic center only if its host cluster can sink to the center before it evaporates. If the cluster dissolves before significant orbital decay occurs, the timescale for the IMBH to fall to the center increases greatly.

According to our estimate in the previous section, around 100 compact clusters have formed close to the center of M82 in the last 200 Myrs. If we assume that half of these clusters contain $`5000\mathrm{M}_{}`$ IMBHs, and that these IMBHs actually merge, then the total BH mass at the center of the galaxy would be at least $`2.5\times 10^5\mathrm{M}_{}`$, high enough to be called a supermassive BH, if we take into account the relatively small total mass of M82.

Having demonstrated that $`5000\mathrm{M}_{}`$ IMBHs can form and reach the galactic center in a reasonable timescale, we now turn to the question of whether the multiple IMBHs at the center can merge. Begelman, Blandford, and Rees (1980) discussed the evolution of a supermassive BH binary at the center of a galaxy, taking into account dynamical friction from field stars and energy loss via gravitational radiation. They found that the merging timescale depends strongly on mass, and for the very massive BHs in which they were interested, merging took much longer than a Hubble time.

There are two reasons why they obtained such a long timescale. The first is that they considered SMBHs with masses exceeding $`10^8M_{}`$, and the timescale for orbital decay via gravitational radiation is proportional to the mass, for given orbital velocity. The second reason is that they assumed that “loss-cone depletion” would occur. In other words, they assumed that stars approaching close enough to the BH binary to interact were supplied only though two-body relaxation, which is a very slow process in large elliptical galaxies.

However, for the IMBHs we consider here, the timescale for merging through gravitational radiation is many orders of magnitude shorter than that for the SMBHs considered by Begelman et al. One reason is simply that the IMBH mass is many orders of magnitude smaller. An IMBH of mass $`10^4M_{}`$ has a timescale $`10^4`$ times shorter than a $`10^8M_{}`$ SMBH when compared at the same orbital velocity. A second reason is that loss-cone depletion is not as effective as originally assumed, at least for relatively small BH mass. Therefore IMBHs can reach high orbital velocities in short timescales. Recent extensive numerical simulations (Makino et al., 1993; Makino, 1997; Merritt and Quinlan, 1997) have shown that the hardening of the BH binary through dynamical friction is in fact several orders of magnitude faster than the prediction from loss-cone arguments. A third reason for the long SMBH merger time is the fact that Begelman et al. ignored the effects of eccentricity, which can significantly reduce the merging time scale. Thus, we can safely conclude that the merging time scale for IMBHs with masses less than $`10^4\mathrm{M}_{}`$ is no more than a few Myr.

Since we expect many IMBHs to fall to the center of the galaxy on a timescale of 100 Myr, in principle, slingshot ejection might remove some of them from the region. However, once one BH has become more massive than typical infalling BHs, it becomes extremely unlikely that it will be ejected, since the kick velocity is inversely proportional to the mass (because of momentum conservation). Thus, even though some of the infalling BHs might be ejected, the central BH would continue to grow. Figure 2 summarizes the discussion in this section.

## 5 Discussion

In this paper, we have discussed the implications for our understanding of the SMBH formation mechanism of the recent discovery of an IMBH in M82. Our conclusion is that the IMBH found in M82 plays the role of “missing link” between stellar-mass BHs and SMBHs.

Since we now know that IMBHs exist, it seems natural to expect that SMBHs might be formed from them. We propose that IMBHs are formed in the cores of young compact star clusters through merging of massive stars and BHs formed from them. These compact young clusters sink to the galactic center by dynamical friction. At the same time, they evaporate via thermal relaxation, stellar mass loss and the effect of the parent galaxy’s tidal field. Thus, IMBHs are created and transported to the center of the galaxy, where they eventually merge to form SMBHs.

In the following, we discuss how we might seek to confirm our new scenario. The most direct evidence would be the observation of gravitational radiation from close binary IMBHs or merging IMBHs. LISA (Jafry, Cornelisse, and Reinhard, 1994), when completed, will be able to detect IMBH merging events even at cosmological distances. The event rate for merging of SMBH is estimated to be one per 1–10 years. In our scenario, each SMBH is a product of $`100`$ mergings of IMBHs or IMBH and growing SMBH. Therefore we predict a much higher event rate for IMBH-IMBH and IMBH-SMBH merging, of the order of 1 per month or even 1 per week.

To test our hypothesis, searches for IMBHs in other galaxies are clearly necessary. In our view, IMBHs are likely to form in young compact star clusters created in nuclear starbursts. We predict that coordinated observations of nearby starburst galaxies at IR, X-ray and radio wavelengths, like those performed for M82, will reveal many more candidate IMBHs.

It is also vital to determine internal and external kinematics of the host star clusters of IMBHs. High-dispersion spectroscopy in the IR with large ground-based telescopes such as SUBARU should be able to determine the velocity dispersion of such a star cluster. Observations by HST would easily resolve the cluster and give us detailed information of its structure. Comparison of these results with theoretical models will then determine whether or not runaway merging can actually take place there.

The Ultra Luminous Compact X-ray Sources(Makishima et al., 2000; Colbert, and Mushotzky, 1999) may be directly related to IMBHs. These are luminous compact X-ray sources with Eddington masses exceeding 100 $`\mathrm{M}_{}`$. Most are found in the arms of spiral galaxies. Thus, they are probably related to young star clusters, and might well also have been formed through merging of massive stars. Their relatively low masses may reflect the differences sizes and masses of their parent star clusters.

## Acknowledgments

This work was supported in part by NASA through Hubble Fellowship grant HF-01112.01-98A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA under contract NAS 5-26555 and by the Research for the Future Program of Japan Society for the Promotion of Science (JSPS-RFTP97P01102). SLWM and SPZ acknowledge support of NASA ATP grants NAG5-6964 and NAG5-9264.


text:
# Anomalous Spectral Statistics of the Asymmetric Rotor Model

## 1 Introduction

As is well known, in the semiclassical limit there is a clear connection between the behavior of classical systems (regular or chaotic) and the corresponding quantal ones. In particular for quantal systems, corresponding to classical regular systems, the spectral statistics ($`P(s)`$ and $`\mathrm{\Delta }_3(L)`$) follow the Poisson ensemble, while for systems corresponding to chaotic ones the Wigner ensemble is followed (see, for example, Ref. 3 and references therein).

Nevertheless, some exceptions are known. The most famous case is perhaps the harmonic oscillator one, discussed in great detail in Ref. 4 and 5.

The aim of this paper is to discuss another pathological case: the classically integrable triaxial rotator model (see, for instance, Ref. 6). Incidentally, this model has been used very often in the description of the low-lying states of the even-even atomic nuclei.

The asymmetric top described by the rotor model is a classically integrable system, but an analytical formula for its quantum energy spectrum it not known. Nevertheless, numerical results can be obtained. By following the Landau approach, the Hamiltonian operator is split into 4 submatrices, corresponding to different symmetry classes. Each truncated submatrix is numerically diagonalized. Finally, the nearest-neighbor spacing distribution $`P(s)`$ and the spectral rigidity $`\mathrm{\Delta }_3(L)`$ are calculated. Surprisingly, the spectral statistics of energy levels do not follow the predictions of Poisson statistics. This result may suggest the presence of a hidden symmetry in the system but, to our knowledge, such a symmetry has not yet been identified.

## 2 The Asymmetric Rotor Model

Let us consider a system of coordinates with axes along the three principal axes of intertia of the top, and rotating with it. The classical Hamiltonian $`H`$ of the top is given by

$$\widehat{H}=\frac{1}{2}\left(aJ_1^2+bJ_2^2+cJ_3^2\right),$$
(1)
where $`\stackrel{}{J}=(J_1,J_2,J_3)`$ is the angular momentum of the rotation and $`a=1/I_1`$, $`b=1/I_2`$, $`c=1/I_3`$ are three parameters such that $`I_1`$, $`I_2`$ and $`I_3`$ are the principal momenta of intertia of the top. The Hamiltonian is classically integrable and its action variables are precisely the three components $`J_s`$, $`s=1,2,3`$, of the angular momentum.

The quantum Hamiltonian $`\widehat{H}`$ is obtained by replacing the components of the angular momentum, in the classical expression of the energy, by the corresponding quantum operators $`\widehat{J}_1`$, $`\widehat{J}_2`$ and $`\widehat{J}_3`$. The commutation rules for the operators of the angular momentum components in the rotating system of coordinates are given by

$$\widehat{J}_r\widehat{J}_s\widehat{J}_s\widehat{J}_r=i\mathrm{}ϵ_{rst}\widehat{J}_t,$$
(2)
where $`ϵ_{rst}`$ is the Ricci tensor and $`r,s,t=1,2,3`$. Note that these commutation rules differ from those in the fixed system in the sign on the right-hand side.

As usual, the two operators $`\widehat{J}^2=\widehat{J}_1^2+\widehat{J}_2^2+\widehat{J}_3^2`$ and $`\widehat{J}_3`$ are simultaneously diagonalized on the basis of eigenstates $`|J,k`$ with integer eigenvalues $`J`$ and $`k`$ ($`k=J,J+1,\mathrm{},J1,J`$), respectively. The non-zero matrix elements of the quantum Hamiltonian $`\widehat{H}`$ in the basis $`|J,k`$ are given by

$$J,k|\widehat{H}|J,k=\frac{\mathrm{}^2}{4}(a+b)(J(J+1)k^2)+\frac{\mathrm{}^2}{2}ck^2,$$
(3)
$$J,k|\widehat{H}|J,k+2=J,k+2|H|J,k=$$

$$=\frac{\mathrm{}^2}{8}(ab)\sqrt{(Jk)(Jk1)(J+k+1)(J+k+2)}.$$
(4)
The quantum Hamiltonian $`\widehat{H}`$ has matrix elements only for transitions with $`kk`$ or $`k\pm 2`$. The absence of matrix elements for transitions between states with even and odd $`k`$ has the result that the matrix of degree $`2J+1`$ is the direct product of two matrices of degrees $`J`$ and $`J+1`$. One of these contains matrix elements for transitions between states with even $`k`$, and the other contains those for transitions between states with odd $`k`$.

It is useful to introduce a new basis, given by

$$|J,k;S=\frac{1}{\sqrt{2}}\left(|J,k+|J,k\right),|J,0,S=|J,0,$$

$$|J,k;A=\frac{1}{\sqrt{2}}\left(|J,k|J,k\right),k0.$$
(5)
By using this new basis, the total Hamiltonian matrix is decomposed in the direct product of $`4`$ submatrices by considering the parity of the quantum number $`k`$: even (E) or odd (O), and the symmetry of the state: symmetric (S) or anti-symmetric (A). So the submatrices are labelled as follow: (E,S), (E,A), (O,S), (O,A). These are the classes of symmetry of the system. In Table 1 we show the dimension of each submatrix for a fixed $`J`$.

The matrix elements of the Hamiltonian $`\widehat{H}`$ in the new basis, with respect to the old basis, are given by

$$J,k,S|\widehat{H}|J,k,S=J,k,A|\widehat{H}|J,k,A=J,k|\widehat{H}|J,k,k1$$
(6)
$$J,1,S|\widehat{H}|J,1,S=J,1|\widehat{H}|J,1+J,1|\widehat{H}|J,1$$
(7)
$$J,1,A|\widehat{H}|J,1,A=J,1|\widehat{H}|J,1J,1|\widehat{H}|J,1$$
(8)
$$J,k,S|\widehat{H}|J,k+2,S=J,k,A|\widehat{H}|J,k+2,A=J,k|\widehat{H}|J,k+2,k0$$
(9)
$$J,0,S|\widehat{H}|J,2,S=\sqrt{2}J,0|\widehat{H}|J,2,k0$$
(10)
We calculate the eigenvalues of each submatrix for different values of $`J`$ using a fast implementation, in double precision, of the Lanczos algorithm with a LAPAC code.

In Figure 1 we plot the density of levels $`\rho (E)`$ of each submatrix of $`\widehat{H}`$ and $`J=1000`$. The results show that the density of levels is practically the same for the four classes. $`\rho (E)`$ displays a high peak at the left-center of the energy interval and a long tail for large energy values.

## 3 Spectral Statistics

As previously discussed, according to the Berry-Tabor theory, given a classical integrable Hamiltonian that, written in action variables $`J_r`$, satisfies the condition $`\left|\frac{^2H}{J_rJ_s}\right|0`$, then, in the semiclassical limit, its spectral statistics should follow the Poisson statistics. Note that a system of linear harmonic oscillators, whose Hamiltonian is given by $`H=\stackrel{}{\omega }\stackrel{}{I}`$, does not satisfy the previous condition. In fact, a system of linear harmonic oscillators is integrable but it does not follow Poissonian statistics.

The triaxial rigid rotator is integrable and satisfies the previous Berry-Tabor condition. Thus, one expects that the spectral statistics of the quantized rigid rotator should be Poissonian. We shall show that is not the case.

In general, various statistics may be used to show the local correlations of the energy levels but the most used spectral statistics are $`P(s)`$ and $`\mathrm{\Delta }_3(L)`$. $`P(s)`$ is the distribution of nearest-neighbor spacings $`s_i=(\stackrel{~}{E}_{i+1}\stackrel{~}{E}_i)`$ of the unfolded levels $`\stackrel{~}{E}_i`$. It is obtained by accumulating the number of spacings that lie within the bin $`(s,s+\mathrm{\Delta }s)`$ and then normalizing $`P(s)`$ to unit. As shown by Berry and Tabor, for quantum systems whose classical analogs are integrable, $`P(s)`$ is expected to follow the Poisson distribution

$$P(s)=\mathrm{exp}(s).$$
(11)
The statistic $`\mathrm{\Delta }_3(L)`$ is defined, for a fixed interval $`(L/2,L/2)`$, as the least-square deviation of the staircase function $`N(E)`$ from the best straight line fitting it:

$$\mathrm{\Delta }_3(L)=\frac{1}{L}\underset{A,B}{\mathrm{min}}_{L/2}^{L/2}[N(E)AEB]^2𝑑E,$$

where $`N(E)`$ is the number of levels between E and zero for positive energy, between $`E`$ and zero for negative energy. The $`\mathrm{\Delta }_3(L)`$ statistic provides a measure of the degree of rigidity of the spectrum: for a given interval L, the smaller $`\mathrm{\Delta }_3(L)`$ is, the stronger is the rigidity, signifying the long-range correlations between levels. For this statistic the Poissonian prediction is

$$\mathrm{\Delta }_3(L)=\frac{L}{15}.$$
(12)
It is useful to remember that Berry, on the basis of the Gutwiller semiclassical formula for the density of states, has shown that $`\mathrm{\Delta }_3(L)`$ deviates from the universal Poissonian predictions for large $`L`$: $`\mathrm{\Delta }_3(L)`$ should saturate to an asymptotic value performing damped oscillations.

In Figure 2 the spectral statistic $`P(s)`$ is plotted for the four submatrices of $`\widehat{H}`$ and $`J=1000`$. Note that the level spectrum is mapped into unfolded levels with quasi-uniform level density by using a standard procedure described in Ref. 13. As expected from the previous analysis of density of levels, $`P(s)`$ is practically the same for the four classes of symmetry. Moreover, $`P(s)`$ has a pathological behavior: a peak near $`s=1`$ and nothing elsewhere.

Compared to $`P(s)`$, the spectral rigidity $`\mathrm{\Delta }_3(L)`$ is less pathological. As shown in Figure 3, $`\mathrm{\Delta }_3(L)`$ follows quite well the Poisson prediction $`\mathrm{\Delta }_3(L)=L/15`$ for small $`L`$ but for larger values of $`L`$ it gets a constant mean value with fluctuations around this mean value. These fluctuations becomes very large by increasing $`L`$, in contrast with the Berry prediction. Note that we have repeated the caculations with other values of $`a`$, $`b`$ and $`c`$ but the results do not change.

The behavior of the density of levels $`\rho (E)`$ and of the spectral statistics $`P(s)`$ and $`\mathrm{\Delta }_3(L)`$ does not change by changing the matrix dimension, namely the quantum number $`J`$. In Figure 4 we plot the density of levels and the spectral statistics for $`J=2000`$ and $`J=4000`$.

## Conclusions

The main conclusion of this paper is that the asymmetric rotor is, like the harmonic oscillator, another pathological case with respect to the classical-quantum correspondence between integrability and Poisson statistics. In our opinion, the pathology of the asymmetric rotor model is more interesting because, unlike the harmonic oscillator, the asymmetric rotor satisfies the conditions of the Berry-Tabor theory. The presence of hidden symmetries could explain the pathological behavior of spectral statistics but such symmetries have not yet been identified.

V.R.M. is greately indebted to Y. Alhassid, M.V. Berry, S. Graffi and V. Zelevinsky, for enlightening discussions. L.S. acknowledges M. Robnik for fruitful conversations.


text:
# Transition Between Ground State and Metastable States in Classical 2D Atoms

## I INTRODUCTION

Wigner suggested in 1934 that a liquid to solid phase transition should occur in a three-dimensional (3D) Fermi system at low densities . The quest for the observation of such a Wigner crystal has been the object of very intense and not stagnant work. After the first discovery of Wigner crystallization of electrons on the surface of liquid helium , there has been considerable theoretical and experimental progress in the study of the mesoscopic system consisting of a finite number of charged particles, which are laterally confined by a parabolic potential and repel each other through a Coulomb potential. This system is the classical analog of the well-known quantum dot problem. These quantum dots are atomic-like structures which have interesting optical properties and may be of interest for single electron devices. These systems and their configurations have been observed experimentally, and are important in solid-state physics, plasma physics as well as in atomic physics. The classical approach is valid for quantum dots in high magnetic fields where the kinetic energy of the electrons is quenched, or for other classical systems, such as laser cooled ions in a trap  which are realized by electric and magnetic fields, trapped ions cooled by laser techniques , ions in a radio-frequency (RF) trap (Paul trap)  or a Penning trap  which can also serve as an illustration of 3D Coulomb clusters . Very large Coulomb clusters have been created recently in strongly coupled RF dusty plasmas  which are like a two-dimensional (2D) layered system. Examples of 2D Coulomb clusters are electrons on the surface of liquid helium  and electrons in quantum dots . The vortex clusters in an isotropic superfluid , vortices in superfluid He4, vortices in a Bose-Einstein condensate stirred with a laser beam  and in superconducting grains  have many common features with those of 2D charged particles . Colloidal particles dissolved in water  and placed between two glass plates are another example of an experimental system where classical particles exhibit Wigner crystallization . Very recently, macroscopic 2D Wigner islands, consisting of charged metallic balls above a plane conductor were studied and ground state, metastable states and saddle point configurations were found experimentally .

In a finite system there is a competition between the bulk triangular lattice and the circular confinement potential which tries to force the particles into a ring like configuration. Those configurations were systematically investigated in Ref.  and a Mendeleev-type of table for these classical atomic-like structures was constructed. The spectral properties of the ground state configurations were presented in Refs.  and generalized to screened Coulomb  and logarithmic  interparticle interactions.

In the present paper we want to go one major step further and calculate not only all the different metastable states but also the saddle points between those local energy minima and the path followed by the particles to transit between those energy minima. The present work is motivated by recent experimental work  where it was found that: i) some of the configurations did not agree with the previous theoretical published one, and ii) they were able to observe some of the saddle points which are the key configurations for transition between different stable (ground or metastable) states. Therefore, we also investigated the stability of the ground state configurations against the functional form of the confinement potential and the exact form of the inter-particle interaction potential.

The present paper is organized as follows. In Sec. 2, we describe the model system. In Sec. 3, our numerical technique, used to obtain the ground and metastable states, is outlined. The technique we used to find the saddle point is similar to the Cerjan-Miller algorithm . After the saddle points are found, we connect the saddle point to the global minimum or a local minimum by the ‘walking downhill’ method. Sec. 4 is devoted to the structural and static properties of the ground and metastable states for $`N=140`$. The configurations are analyzed and compared with available experimental data and the results of previous theoretical approaches. The phase diagram for 9 and 16 particles in the ground state with different functional forms of confinement potential and interparticle interaction is also calculated. The discussion on the saddle point is presented in Sec. 5, and the connecting path from the ground state to the metastable states is found, and we investigate the completely geometric properties of the energy landscape. Our conclusions are presented in Sec. 6.

## II Model system

The model system consists of identical charged particles interacting through a Coulomb repulsive interacting and moving in a 2D plane where they are confined by a parabolic potential

$$H=\frac{q^2}{\epsilon }\underset{i>j}{}\frac{1}{\left|\stackrel{}{r_i}\stackrel{}{r}_j\right|}+\underset{i}{}V(\stackrel{}{r}_i).$$
(1)
The confinement potential $`V(\stackrel{}{r})=\frac{1}{2}m^{}\omega _0^2r^2`$ is taken circular symmetric and parabolic, where $`m^{}`$ is the effective mass of the particles, $`q`$ is the particle charge, $`\omega _0`$ is the radial confinement frequency and $`\epsilon `$ is the dielectric constant of the medium the particles are moving in. Note that for the quantum dot problem an additional term appears in Eq. (1) which is the kinetic energy of the particles which is absent in our statical classical problem. Here the motion of the particles is restricted to the $`(x,y)`$ plane. To exhibit the scaling of the system, we introduce the characteristic scales in the problem: $`r_0=(2q^2/mϵ\omega _0^2)^{1/3}`$ for the length and $`E_0=(m\omega _0^2q^4/2ϵ^2)^{1/3}`$ for the energy. After the scaling transformations ($`rr/r_0,EE/E_0`$), the Hamiltonian can be rewritten in a simple dimensionless form as

$$H=\underset{i>j}{}\frac{1}{\left|\stackrel{}{r_i}\stackrel{}{r}_j\right|}+\underset{i}{}V(\stackrel{}{r_i}),$$
(2)
with $`V(\stackrel{}{r})=x^2+y^2`$ and which only depends on the number of particles $`N`$. The numerical values for the parameters $`\omega _0,r_0,E_0`$ for some typical experimental systems were given in Ref. .

## III Numerical approach

Most of the previous works have treated the quantum mechanical problem of a small number of electrons. In the present paper, we consider only the classical system. Although a classical approach for the description of the behavior of electrons in quantum dots, in principle, is not applicable, it is possible that certain features of the classical system may survive in a quantum system.

Due to the presence of confinement energy and electron-electron Coulombic interaction, a complete description of the cluster system is complicated and can’t be obtained analytically. The Monte Carlo simulation technique  is relatively simple and rapidly convergent and it provides a reliable estimation of the total energy of the system in cases when relatively small number of Metropolis steps is sufficient. However, the accuracy of this method in calculating the explicit states is poor in certain cases. It becomes more difficult for clusters with a large number of particles, which have significantly more metastable states. To circumvent this problem we employ the numerical technique of Newton optimization which was outlined and compared with the standard Monte Carlo technique in Ref. . In this way, we are able to obtain not only the ground state but also metastable states. It also yields the eigenfrequencies and the eigenmodes of the ground state configuration. Now only a small number of calculation steps is needed to obtain the same accuracy. Moreover, using the modified Newton approach, we can explore the stability of the system in its ground-state configuration through its spectrum.

By studying the characteristics of the energy landscape and the energy barrier between the different local minima, we are able to find the saddle point configurations which are very important and are the key configurations for transition between different stable states. The technique we used to find the saddle point is explained in more detail in Ref. , and is similar to the Cerjan-Miller algorithm . After the saddle points are found, we connect the saddle point to the global minimum or a local minimum by the ‘walking downhill’ method. In this algorithm the direction of the steepest gradient is followed to force the system to transit from the saddle point state to the local minimum state. Which minimum is finally reached depends on the initial step, therefore we repeat this procedure several times to determine both minima which the saddle point state connects. Thus the connecting path followed by the particles to transit between those energy minima is found, from which the geometric properties of the energy landscape are obtained.

## IV Ground state and metastable state

In Table I we list for $`N`$=1, 2,…,40 the energy per particle $`E/N`$ in the ground state and in the metastable states, where we also list the energy difference with the ground state $`\mathrm{\Delta }E/N.`$ The configuration is indicated by the number of particles in the different rings, the position of the center of the ring and the radius of the different rings, the width of the ring which is defined as the difference of the maximum radius and minimum radius in the same ring, and the energy of the lowest three normal mode frequencies of the ground state are also given in Table I. This table is rather exhaustive and should be compared with a similar one published in Ref.  for a logarithmic interacting system.

For different values of $`N`$ there exist different possible values for $`E/N`$ which are nothing else than the metastable states. The difference in energy between the metastable and the ground state is given in the third column and the corresponding configuration in the fourth column. Note that with increasing $`N`$ the number of metastable configurations increases and in general (but not always) the widths of the rings for metastable configurations are larger and the central ring/particle is not exactly located in the center of the parabolic potential well. For sufficiently large $`N`$, the simple ring structure gradually disappears in the center and the triangular Wigner lattice appears. There is a competition between two types of ordering: ordering into a triangular-lattice structure (Wigner lattice) and ordering into a shell structure, which leads to clusters with interesting self-organized patterns which show concentric shells at small $`N`$ and hexagonal cores surrounded by circular outer shells at large $`N`$.

The lowest non-zero normal mode frequency is a measure for the stability of the ground state, it tells us how easy or difficult it is to deform this state. Therefore, intuitively we would expect that the value of this frequency would be correlated to $`\mathrm{\Delta }E/N`$, the energy difference between the first metastable and the ground state. Those values are plotted in Fig. 1 as function of $`N`$. Notice that there exist such a correlation in general, but that this is not true for all $`N`$-values, e.g. for $`N=12,18,19,20,21,30`$ there is no correlation.

The rings have sometimes a finite width which are shown in Fig. 2 as a function of $`N`$. Notice that the widths fall into three bands; i) width $``$ 0.003 which is practically a perfect ring, ii) width $``$ 0.02, and iii) width $``$ 0.5. Usually, but not always, the outer ring has the largest width. The width of the rings increases with increasing $`N`$ and at the same time the widest ring becomes often the next to outer ring.

We compare our ground state configuration with available experimental data  and the results of previous theoretical approaches . For very small number of particles ($`N<16`$), all theoretical and experimental results for the ground state configurations are the same expect for $`N`$=9 and 15 whatever kind of interparticle interaction. The experimental observation  for the ground state of 9 particles is (1,8) and for 15 particles it is (4,11), which compares to our result (2,7) and (5,10), respectively. For $`17<N<30,`$ the experimental result and all the calculated patterns present three shells. Our result differs with the experimental data of Ref.  for $`N`$= 17, 20, 22, 24, 25, 27-30. Because of the discrepancy between some of the experimental configurations and the ‘numerical exact’ theoretical ground state configurations it is possible that experimentally the inter-particle interaction is not exactly a Coulombic potential and the confinement potential is not purely quadratic. Therefore, we investigated the effect of such deviations of these potentials on the ground state configuration. As an example we took $`N`$=9 and use confinement potentials $`Vr^n`$ and for the inter-particle interaction $`Vr^m`$. The resulting phase diagram is shown in Fig. 3. Notice that, depending on the values of $`n`$ and $`m,`$ the system can be either in the (1,8) or the (2,7) configuration. For the harmonic confined Coulomb interacting system i.e. ($`n,m`$) = (2,1) the system is in the (2,7) configuration but from the phase diagram it is clear that if we change the confinement potential slightly and make it more steep up to $`n`$2.2 the configuration (1,8) becomes the ground state. The experimentally determined ground state configuration for 9 particles was (1,8) .

There is also a difference with the experimental data and our results for $`N`$=16 particles. Therefore, we did the same investigation and present the phase diagram in Fig. 4. Notice that the harmonic confined Coulomb interacting system, i.e. $`(n,m)`$ = (2,1) is again close to the phase boundary between the configuration (1,5,10) and (5,11). This is probably the explanation why the experimental configuration differs from our simulation results, since it is hard to guarantee that ($`n,m`$) is exactly (2,1) during the experiment.

Notice that for both $`N`$=9 and $`N`$=16 the metastable configuration has an energy very close to the one of the ground state, the difference is less than 0.2$`\%`$. These metastable configurations correspond indeed with the experimentally observed ones. Consequently, an alternative explanation for the difference with the experiment is that the experimental configuration got stuck in the metastable configuration.

## V Saddle points

Between metastable states and the ground state there are potential barriers. The system will prefer to transfer over the lowest potential barrier, which is the saddle point configuration between these energy minima, in order to transit from one stable configuration to the other. We plot in Fig. 5 the trajectories of the particles for the $`N`$=5 system making a transition from the ground state (5) to the metastable state (1,4) and the saddle point connecting them. The trajectories of the particles can also be obtained by moving one of the particles to the center of the system.

For 6 particles, the ground state (1,5) and the metastable state (6), corresponding to the hexagonal configuration, are obtained. Moreover, the unstable equilibria associated to saddle point configurations are also obtained, and the energy landscape is shown schematically in Fig. 6. There are two saddle points for this case, one of them is very close to the metastable state in both energy and configuration, and will therefore be hard to see experimentally . In Fig. 6, the insets show the arrangement of the particles for the different states. Using the ‘walking downhill method’, we found the central particle slowly moving to the periphery of the cluster. We would like to stress that the configuration with 6 particles on a perfect ring is a saddle point state in contrast to the claim made in Ref. . This can be understood from the following simple model calculation: if 3 particles are placed on a circle with radius $`A`$, on the corners of an equilateral triangle, and the other 3 particles on another equilateral triangle’s corners with radius $`B`$ rotated over $`60^0`$, the energy is

$$E(c)=\frac{9}{2}\left(\frac{1+c^2}{36}\right)^{1/3}(\frac{1}{1+c}+\frac{1+c}{\sqrt{3}c}+\frac{2}{\sqrt{1c+c^2}})^{2/3}$$
(3)
where $`c=\frac{A}{B}`$. This function is shown in Fig. 7. It is clear that the perfect circle configuration i.e. $`c=\frac{A}{B}=1`$ is a saddle point, and that the minimum is obtained if 3 particles move a bit to the center, and the other 3 particles move away from the center (see the insets in Fig. 7). Both shown metastable states are just connected by a rotation over $`120^0`$. The two minima in Fig. 7 correspond to the same configuration in which inner and outer ring are interchanged. Comparing our results with the Fig. 2 (‘$`N=6`$: ground state, saddle point configuration and the hexagonal metastable state’) of Ref. , we see that the other saddle point is observed experimentally.

A list of the saddle point energies up to 20 particles is given in Table II. From this table, we notice that there is only one saddle point state for $`N=3,4,5`$ particles. But, on the other hand it is well-known that there are $`(k1)`$ saddle points when there are $`k`$ minima. For $`N`$=3 and 4 one saddle point is found, although there is no metastable configuration. The reason is that the saddle point state connects two equilateral ground state configurations which can be obtained from each other by a simple rotation. For the simple case of 3 particles, we show the energy surface and the corresponding configurations schematically in Fig. 8. Notice that there are always more saddle points than minima for $`N`$$`>`$6. With increasing the number of particles, more saddle point states are obtained and the energy landscape gets more complicated. For example for 9 particles, we obtain three saddle points and one metastable state. The results for the trajectories and energy landscape are shown in Fig. 9. Again, the ground state configurations corresponding with the black and the white dot are connected by a simple rotation, i. e. a symmetry operation.

## VI Conclusion

We presented the results of a numerical calculation of the configurations of the ground and all metastable states and their energies, the system consisting of classical 2D charged particles that are confined in a parabolic confinement potential for $`N`$=1,…,40. These artificial atoms undergo configurational changes when the system transits from the ground state to the different metastable states, or between the different metastable states. Such transitions move through the lowest energy barrier connecting those states, i. e. through a saddle point. The connecting path from the ground state to all metastable states is found and the geometric properties of the energy landscape were discussed.

Sensitivity of the configuration on the form of the confinement potential and the interparticle interaction is investigated and a phase diagram was obtained. This sensitivity on e.g. the form of the confinement potential is probably the explanation why the experimental configuration differs from our simulation results.

## VII Acknowledgments

B. Partoens is a post-doctoral researcher of the Flemish Science Foundation (FWO-Vlaanderen). Stimulating discussions with Dr. J. Shi and M. Milošević are gratefully acknowledged. This work is supported by the Flemish Science Foundation (FWO-VI), the Belgian Inter-University Attraction Poles (IUAP-VI), the “Onderzoeksraad van de Universiteit Antwerpen”(GOA), the EU Research Training Network on “Surface Electrons on Mesoscopic Structures”, and INTAS.

Tables

Figure captions


text:
# 1 Leading chargino – up-squark contribution to 𝐾-𝐾̄ mixing.

SUSX-TH/01-022

Chargino Contributions to $`\epsilon `$ and $`\epsilon ^{}`$
Shaaban Khalil1,2 and Oleg Lebedev1

1Centre for Theoretical Physics, University of Sussex, Brighton BN1 9QJ, UK
2Ain Shams University, Faculty of Science, Cairo, 11566, Egypt

Abstract
We analyze the chargino contributions to the $`K\overline{K}`$ mixing and $`\epsilon ^{}`$ in the mass insertion approximation and derive the corresponding bounds on the mass insertion parameters. We find that the chargino contributions can significantly enlarge the regions of the parameter space where CP violation can be fully supersymmetric. In principle, the observed values of $`\epsilon `$ and $`\epsilon ^{}`$ may be entirely due to the chargino – up-squark loops.

A convenient way to parameterize SUSY contributions to the flavor changing processes is to employ the so called mass insertion approximation . The advantage of this approach is that it allows to treat such contributions in a model independent way without resorting to specific assumptions about the SUSY flavor structures (a technical definition of this approximation will be given below).

The gluino contributions to the kaon observables in the mass insertion approximation have been studied in detail -, but the chargino contributions have not received similar attention. The latter have been considered either in the context of minimal flavor violation, that is in SUSY models with the flavor mixing given by the CKM matrix ,  or as contributing to the $`K\overline{K}`$ mixing only . In general, the flavor structure in the squark sector may be very complicated. In particular, flavor patterns in the up and down sectors can be entirely different, which may result in the dominance of the chargino contributions to the K and B observables. On the other hand, the neutralino contributions involve the same mass insertions as the gluino ones (i.e. down type mass insertions) and thus cannot qualitatively change the picture. In this Letter, we study the chargino contributions to the $`K\overline{K}`$ mixing and $`\epsilon ^{}`$ using the mass insertion approximation and derive the corresponding bounds on the mass insertion parameters.

Let us first consider the $`K\overline{K}`$ mixing. The two observables of primary interest are the $`K_L`$-$`K_S`$ mass difference and indirect CP violation in $`K\pi \pi `$ decays:

$`\mathrm{\Delta }M_K=M_{K_L}M_{K_S},`$
$`\epsilon ={\displaystyle \frac{A(K_L\pi \pi )}{A(K_S\pi \pi )}}.`$ (1)
The experimental values for these parameters are $`\mathrm{\Delta }M_K3.489\times 10^{15}`$ GeV and $`\epsilon 2.28\times 10^3`$. The Standard Model predictions for them lie in the ballpark of the measured values, however a precise prediction cannot be made due to the hadronic and CKM uncertainties.

Generally, $`\mathrm{\Delta }M_K`$ and $`\epsilon `$ can be calculated via

$`\mathrm{\Delta }M_K=2\mathrm{R}\mathrm{e}K^0|H_{\mathrm{eff}}^{\mathrm{\Delta }S=2}|\overline{K}^0,`$
$`ϵ={\displaystyle \frac{1}{\sqrt{2}\mathrm{\Delta }M_K}}\mathrm{Im}K^0|H_{\mathrm{eff}}^{\mathrm{\Delta }S=2}|\overline{K}^0.`$ (2)
Here $`H_{\mathrm{eff}}^{\mathrm{\Delta }S=2}`$ is the effective Hamiltonian for the $`\mathrm{\Delta }S=2`$ transition. It can be expressed via the Operator Product Expansion (OPE) as

$$H_{\mathrm{eff}}^{\mathrm{\Delta }S=2}=\underset{i}{}C_i(\mu )Q_i,$$
(3)
where $`C_i(\mu )`$ are the Wilson coefficients and $`Q_i`$ are the relevant local operators. The main uncertainty in this calculation arises from the matrix elements of $`Q_i`$, whereas the Wilson coefficients can be reliably calculated at high energies and evolved down to low energies via the Renormalization Group (RG) running.

In supersymmetric extensions of the SM, the dominant chargino contribution to $`H_{\mathrm{eff}}^{\mathrm{\Delta }S=2}`$ comes from the “super-box” diagram in Fig.1. We perform our calculations in the super CKM basis, i.e. the basis in which the gluino-quark-squark vertices are flavor-diagonal. In this basis, the chargino – left quark – left squark vertices involve the usual CKM matrix:

$$\mathrm{\Delta }=g\underset{k}{}\underset{a,b}{}V_{k1}K_{ba}^{}d_L^ai\sigma _2(\stackrel{~}{\chi }_{kL}^+)^{}\stackrel{~}{u}_L^b,$$
(4)
where $`K`$ is the CKM matrix, $`a,b`$ are the flavor indices, $`k=1,2`$ labels the chargino mass eigenstates, and $`V,U`$ are the chargino mixing matrices defined by

$`M_{\chi ^+}=\left(\begin{array}{cc}M_2& \sqrt{2}M_W\mathrm{sin}\beta \\ \sqrt{2}M_W\mathrm{cos}\beta & \mu \end{array}\right),`$
$`U^{}M_{\chi ^+}V^1=\mathrm{diag}(m_{\chi _1^+},m_{\chi _2^+}).`$ (5)
Only the gaugino components of the charginos lead to significant contributions to the $`K\overline{K}`$ mixing since the higgsino couplings are suppressed by the quark masses (except for the stop coupling) and are not important even at large $`\mathrm{tan}\beta `$. The stop loop contribution is suppressed by the CKM mixing at the vertices: each vertex involving the stop is suppressed by $`\lambda ^2`$ or $`\lambda ^3`$ with $`\lambda `$ being the Cabibbo mixing, whereas we will be working in $`𝒪(\lambda )`$ order. The super-box involving higgsino interactions with the stops depends on the left-right mass insertions and, as will be clear later, does not lead to useful constraints on the SUSY flavor structures.

Due to the gaugino dominance, chargino-squark loops will generate a significant contribution to only one operator

$$Q_1=\overline{s}_L^\alpha \gamma ^\mu d_L^\alpha \overline{s}_L^\beta \gamma _\mu d_L^\beta ,$$
(6)
similarly to the Standard Model ($`\alpha ,\beta `$ are the color indices). The corresponding Wilson coefficient is given by the sum of the Standard Model and the chargino contributions:

$$C_1(\mu )=C_1(\mu )^{SM}+C_1(\mu )^{\stackrel{~}{\chi }^+}.$$
(7)
Generally, there are additional contributions from gluinos and the Higgs sector, but they are not correlated with the chargino contributions and are unimportant for the present study. We calculate $`C_1(\mu )^{\stackrel{~}{\chi }^+}`$ using the mass insertion approximation. That is, we express the left-left squark propagator as

$$\stackrel{~}{u}_L^a\stackrel{~}{u}_L^b=\mathrm{i}(k^2\mathrm{𝟏}m^2\mathrm{𝟏}\delta m^2)_{ab}^1\frac{\mathrm{i}\delta _{ab}}{k^2m^2}+\frac{\mathrm{i}(\delta m^2)_{ab}}{(k^2m^2)^2},$$
(8)
where $`\mathrm{𝟏}`$ is the unit matrix and $`m`$ is the average up-squark mass. The SUSY contributions are parameterized in terms of the dimensionless parameters $`(\delta _{LL}^u)_{ab}(\delta m^2)_{ab}/m^2`$. The corresponding Wilson coefficient is calculated to be

$`C_1(M_W)^{\stackrel{~}{\chi }^+}={\displaystyle \frac{g^4}{768\pi ^2m^2}}\left({\displaystyle \underset{a,b}{}}K_{a2}^{}(\delta _{LL}^u)_{ab}K_{b1}\right)^2{\displaystyle \underset{i,j}{}}|V_{i1}|^2|V_{j1}|^2{\displaystyle \frac{x_ih(x_i)x_jh(x_j)}{x_ix_j}},`$ (9)
where $`x_im_{\stackrel{~}{\chi }_i^+}^2/m^2`$ and

$$h(x)=\frac{2+5xx^2}{(1x)^3}+\frac{6x\mathrm{ln}x}{(1x)^4}.$$
(10)
It is interesting to note that “flavor-conserving” mass insertions $`(\delta _{LL}^u)_{aa}`$ contribute to $`C_1(M_W)`$, unlike for the gluino case. Such mass insertions arise from non-degeneracy of the squark masses and are proportional to the difference of the average squark mass squared and the diagonal matrix elements of the squark mass matrix. If the diagonal elements are equal, the “flavor-conserving” mass insertions drop out of the sum due to the GIM cancellations.

The flavor structure appearing in Eq.9 can be expanded in powers of $`\lambda `$:

$$\underset{a,b}{}K_{a2}^{}(\delta _{LL}^u)_{ab}K_{b1}=(\delta _{LL}^u)_{21}+\lambda \left[(\delta _{LL}^u)_{11}(\delta _{LL}^u)_{22}\right]+𝒪(\lambda ^2).$$
(11)
Assuming the presence of one type of the mass insertions at a time in Eq.11 at each order in $`\lambda `$, one can derive constraints on $`(\delta _{LL}^u)_{21}`$ and $`\delta (\delta _{LL}^u)_{11}(\delta _{LL}^u)_{22}`$ imposed by $`\mathrm{\Delta }M_K`$ and $`\epsilon `$. A much weaker constraint on $`(\delta _{LL}^u)_{31}`$ can also be obtained if we are to keep $`𝒪(\lambda ^2)`$ terms in Eq.11.

To derive constraints on the mass insertions, one has to take into account the RG evolution of the Wilson coefficients. In our numerical numerical analysis, we use the NLO QCD result $`C_1(\mu )^{\stackrel{~}{\chi }^+}0.8C_1(M_W)^{\stackrel{~}{\chi }^+}`$ with $`\mu =2`$ GeV . The matrix element of $`Q_1`$ is computed via $`K^0|Q_1|\overline{K}^0=\frac{1}{3}M_Kf_K^2B_1(\mu )`$ with the lattice value $`B_1(\mu )=0.61`$ . In addition, the SM contribution should be taken into account. Its detailed discussion can be found in Ref.. In our analysis, we assume a zero CKM phase which corresponds to a conservative bound on the mass insertion. The Wolfenstein parameters are set to $`A=0.847`$ and $`\rho =0.4`$. The other relevant constants are $`M_K=0.498`$ GeV and $`f_K=0.16`$ GeV.

The resulting bounds on $`(\delta _{LL}^u)_{21}`$ and $`\delta `$ as functions of $`M_2`$ and the average squark mass $`m`$ are presented in Tables 1 and 2. We find that these bounds are largely insensitive to $`\mathrm{tan}\beta `$ in the range 3–40 and to $`\mu `$ in the range $`200500`$ GeV. This can be understood since these parameters do not significantly affect the gaugino components of the charginos and their couplings. Note that $`\delta `$ is real due to the hermiticity of the squark mass matrix and therefore does not contribute to $`\epsilon `$. The presented bounds on the real part of $`(\delta _{LL}^u)_{21}`$ are a bit stronger than those derived from the gluino contribution to the $`D\overline{D}`$ mixing , whereas the imaginary part of $`(\delta _{LL}^u)_{21}`$ is not constrained by any other FCNC processes.

In principle, it is possible to constrain the $`(\delta _{LL}^u)_{31}`$ mass insertion as well. At order $`\lambda ^2`$, there are two contributions in Eq.11: from $`(\delta _{LL}^u)_{31}`$ and $`(\delta _{LL}^u)_{12}+(\delta _{LL}^u)_{21}`$. Assuming no cancellations between these two terms, the constraints on $`(\delta _{LL}^u)_{31}`$ are obtained by multiplying the bounds in Tables 1 and 2 by $`(A\lambda ^2)^124`$. Clearly, this leaves the real part of $`\left[(\delta _{LL}^u)_{31}\right]^2`$ essentially unconstrained, while the bound on $`\sqrt{\mathrm{Im}\left[(\delta _{LL}^u)_{31}\right]^2}`$ is of order $`10^1`$. We note that a similar constraint on $`(\delta _{LR}^u)_{13}`$ can be derived from the higgsino-stop contribution, however such a constraint is typically satisfied automatically (especially if the squarks are heavy) since $`(\delta _{LR}^u)_{13}ϵm_t/m`$ with $`ϵ1`$ being the 1-3 mixing the left-right sector.

Next let us consider the chargino contribution to $`\epsilon ^{}`$ using the same approximations. The $`\epsilon ^{}`$ parameter is a measure of direct CP violation in $`K\pi \pi `$ decays given by

$$\frac{\epsilon ^{}}{\epsilon }=\frac{\omega }{\sqrt{2}|\epsilon |\mathrm{Re}A_0}\left(\mathrm{Im}A_0\frac{1}{\omega }\mathrm{Im}A_2\right),$$
(12)
where $`A_{0,2}`$ are the amplitudes for the $`\mathrm{\Delta }I=1/2,3/2`$ transitions and $`\omega \mathrm{Re}A_2/\mathrm{Re}A_01/22`$. Experimentally it has been found to be $`\mathrm{Re}(\epsilon ^{}/\epsilon )1.9\times 10^3`$ which provides firm evidence for the existence of direct CP violation. This value can be accommodated in the Standard Model although the theoretical prediction involves large uncertainties.

The effective Hamiltonian for the $`\mathrm{\Delta }S=1`$ transition is given by

$`H_{\mathrm{eff}}^{\mathrm{\Delta }S=1}={\displaystyle \underset{i}{}}\widehat{C}_i(\mu )\widehat{Q}_i.`$ (13)
Just as in the Standard Model, two operators, $`\widehat{Q}_6`$ and $`\widehat{Q}_8`$, play the dominant role. They originate from the gluon and electroweak pengiun diagrams (Fig.2) and are defined by

$`\widehat{Q}_6=(\overline{s}^\alpha d^\beta )_{VA}{\displaystyle \underset{q=u,d,s}{}}(\overline{q}^\beta q^\alpha )_{V+A},`$
$`\widehat{Q}_8={\displaystyle \frac{3}{2}}(\overline{s}^\alpha d^\beta )_{VA}{\displaystyle \underset{q=u,d,s}{}}e_q(\overline{q}^\beta q^\alpha )_{V+A},`$ (14)
with $`(\overline{f}f)_{VA}\overline{f}\gamma _\mu (1\gamma _5)f`$. Their matrix elements are enhanced by $`(m_K/m_s)^2`$ compared to those of the other operators:

$`(\pi \pi )_{I=0}|Q_6|K^0=4\sqrt{{\displaystyle \frac{3}{2}}}\left[{\displaystyle \frac{m_K}{m_s(\mu )+m_d(\mu )}}\right]^2m_K^2(f_Kf_\pi )B_6,`$
$`(\pi \pi )_{I=2}|Q_8|K^0=\sqrt{3}\left[{\displaystyle \frac{m_K}{m_s(\mu )+m_d(\mu )}}\right]^2m_K^2f_\pi B_8,`$ (15)
where $`B_{6,8}`$ are the bag parameters. In addition, the contributions of these operators are enhanced by the QCD corrections. Although the Wilson coefficient of $`\widehat{Q}_8`$ is suppressed by $`\alpha /\alpha _s`$ compared to that of $`\widehat{Q}_6`$, its contribution to $`\epsilon ^{}`$ is enhanced by $`1/\omega `$ and is significant. In fact, it provides the dominant contribution in our analysis.

The relevant QCD corrections in the context of the MSSM with minimal flavor violation have been studied in Ref. and later, in more detail, in Ref.. To account for a general flavor structure in the mass insertion approximation, only the loop functions of Ref. are to be modified. In our numerical analysis we use the parameterization of Ref. and express the chargino contribution to $`\epsilon ^{}/\epsilon `$ as

$$\left(\frac{\epsilon ^{}}{\epsilon }\right)^{\stackrel{~}{\chi }^+}=\mathrm{Im}\left(\underset{a,b}{}K_{a2}^{}(\delta _{LL}^u)_{ab}K_{b1}\right)F_\epsilon ^{},$$
(16)
where

$$F_\epsilon ^{}=(P_X+P_Y+P_Z)F_Z+\frac{1}{4}P_ZF_\gamma +P_EF_g.$$
(17)
Here we have omitted the box diagram contributions which are negligible . The parameters $`P_i`$ include the relevant matrix elements and NLO QCD corrections, and are given by $`P_X=0.58`$, $`P_Y=0.48`$, $`P_Z=7.67`$, and $`P_E=0.82`$. The quantities $`F_i`$ are functions of supersymmetric parameters resulting from the gluon, photon, and Z penguin diagrams (Fig.2) and are calculated in the mass insertion approximation. Explicitly,

$`F_g`$ $`=`$ $`2{\displaystyle \frac{m_W^2}{m^2}}{\displaystyle \underset{i}{}}|V_{i1}|^2f_g(x_i),`$
$`F_\gamma `$ $`=`$ $`2{\displaystyle \frac{m_W^2}{m^2}}{\displaystyle \underset{i}{}}|V_{i1}|^2f_\gamma (x_i),`$
$`F_Z`$ $`=`$ $`{\displaystyle \frac{1}{8}}2{\displaystyle \underset{i}{}}|V_{i1}|^2f_Z^{(1)}({\displaystyle \frac{1}{x_i}},{\displaystyle \frac{1}{x_i}})`$ (18)
$`+`$ $`2{\displaystyle \underset{i,j}{}}V_{j1}^{}V_{i1}\left[U_{i1}U_{j1}^{}f_Z^{(2)}(x_j,x_i)V_{j1}V_{i1}^{}f_Z^{(1)}(x_j,x_i)\right],`$
where $`x_im_{\stackrel{~}{\chi }_i^+}^2/m^2`$ and the loop functions are given by

$`f_g(x)={\displaystyle \frac{16x+18x^210x^33x^4+12x^3\mathrm{ln}x}{18(x1)^5}},`$
$`f_\gamma (x)={\displaystyle \frac{2260x+45x^24x^33x^4+3(39x^2+4x^3)\mathrm{ln}x}{27(x1)^5}},`$
$`f_Z^{(1)}(x,y)={\displaystyle \frac{(y1)\left[(x1)(x^2x^2y+xy^2y^2)+x^2(y1)\mathrm{ln}x\right](x1)^2y^2\mathrm{ln}y}{16(x1)^2(y1)^2(yx)}},`$
$`f_Z^{(2)}(x,y)=\sqrt{xy}{\displaystyle \frac{(y1)\left[(x1)(xy)+x(y1)\mathrm{ln}x\right](x1)^2y\mathrm{ln}y}{8(x1)^2(y1)^2(yx)}}.`$ (19)
As noted in ref., the dominant contribution typically comes from the Z-penguin diagram, especially if the SUSY particles are heavy. This can be seen as follows. Due to the gauge invariance, the $`g\overline{s}_Ld_L`$ and $`\gamma \overline{s}_Ld_L`$ vertices are proportional to the second power of the momentum transfer, i.e. $`(q_\mu q_\nu g_{\mu \nu }q^2)/m^2`$. This momentum dependence is cancelled by the gluon (photon) propagator which leads to the suppression factor $`1/m^2`$ in the final result. On the other hand, the $`Z\overline{s}_Ld_L`$ vertex exists at $`q^2=0`$ due to the weak current non-conservation and is momentum-independent to leading order. It is given by a dimensionless function of the ratios of the SUSY particles’ masses. The Z propagator then leads to the suppression factor $`1/M_Z^2`$ which is much milder than $`1/m^2`$ appearing in the gluon and photon contributions.

The resulting bounds on $`\mathrm{Im}(\delta _{LL}^u)_{21}`$ are presented in Table 3. Note that there is no SM contribution to $`\epsilon ^{}`$ since we assume a vanishing CKM phase. In the limit of heavy superpartners, these bounds become insensitive to the SUSY mass scale. This occurs due to the dominance of the Z penguin contribution. Indeed, the contributions of the photon and gluon penguins fall off as $`1/(\mathrm{SUSY}\mathrm{scale})^2`$ as can be seen from Eq.18. On the other hand, the Z-penguin contribution stays constant in the “decoupling” limit. This may seem to conflict with the intuitive expectation of the decoupling of heavy particles. However, the proper decoupling behaviour is obtained when the flavor violating parameters $`\delta m_{ab}^2`$ are kept constant (or if they grow slower than the masses of the superpartners). To put it in a slightly different way, the decoupling should be expected when the mass splittings among the squarks grow slower than the masses themselves.

It is noteworthy that (for universal GUT scale gaugino masses of around 200 GeV) the bounds on $`\mathrm{Im}(\delta _{LL}^u)_{21}`$ are slightly stronger than those on $`\mathrm{Im}(\delta _{LL}^d)_{21}`$ derived from the gluino contribution to $`\epsilon ^{}`$ . The suppression due to the weaker coupling is compensated by a larger loop function mainly due to the presence of the diagram on the right in Fig.2.

These results show that to have a chargino-induced $`\epsilon ^{}`$ would require a relatively large LL mass insertion ($`𝒪(10^1)`$) which typically violates the constraints from $`\mathrm{\Delta }M_K`$ and $`\epsilon `$. Yet, it is possible to saturate $`\epsilon `$ and $`\epsilon ^{}`$ with the chargino contributions in corners of the parameter space. For instance, taking $`M_2=450`$ GeV and $`m=300`$ GeV, $`\epsilon `$ requires

$$2\mathrm{Im}(\delta _{LL}^u)_{21}\mathrm{Re}(\delta _{LL}^u)_{21}2.3\times 10^4.$$
(20)
Then, assuming $`\mathrm{Im}(\delta _{LL}^u)_{21}0.12`$ to produce $`\epsilon ^{}`$, $`\mathrm{Re}(\delta _{LL}^u)_{21}`$ has to be $`9\times 10^4`$. These values are in marginal agreement with the $`\mathrm{\Delta }M_K`$ bound:

$$\sqrt{\left|\left[\mathrm{Re}(\delta _{LL}^u)_{21}\right]^2\left[\mathrm{Im}(\delta _{LL}^u)_{21}\right]^2\right|}0.12.$$
(21)
The main lesson, however, is that combining the chargino and the gluino contributions can provide fully supersymmetric $`\epsilon `$ and $`\epsilon ^{}`$ in considerable regions of the parameter space*** In principle it is possible to saturate both $`ϵ^{}`$ and $`ϵ`$ with gluino contributions , see also . It was noted in Ref. that the LR mass insertions of the required size may lead to the charge and color breaking minima. However, this is not true in general, i.e. when the Higgs, squark, and slepton mass parameters are unrelated, see Ref. . . For example, only a small ($`\mathrm{Im}(\delta _{LR}^d)_{21}10^5`$) mass insertion in the down-sector is required to generate the observed value of $`\epsilon ^{}`$ . Then $`\epsilon `$ can be entirely due to the mass insertions in the up-sector: $`\mathrm{Im}(\delta _{LL}^u)_{21}10^3`$ and $`\mathrm{Re}(\delta _{LL}^u)_{21}10^2`$. Generally, this does not require large SUSY CP-phases and may be accommodated in the framework of approximate CP symmetry , which is motivated by the strong EDM bounds (for a review see ). Alternatively, the CP-phases can be order one but enter only flavor-off-diagonal quantities which occurs in models with hermitian flavor structures . Clearly, the regions of the parameter space where CP violation can be fully supersymmetric are significantly enlarged if both the gluino and the chargino contributions are included. Of course, it remains a challenge to build a realistic well-motivated model with all of the required features.