text:
# ITFA-00-13 July 2000 Transversality of the logarithmic divergences in the Classical Finite Temperature 𝑆𝑈(𝑁) Self-Energy

## I Introduction

As it is well known, low-momentum excitations in a high temperature plasma behave classically. Nevertheless, to study these excitations one cannot simply replace the quantum thermal field theory that describes them by a classical thermal field theory, because of the appearance of ultraviolet Rayleigh-Jeans type divergences. One should think of a classical thermal field theory as an effective theory at large scales (larger than the typical interparticle distance $`\mathrm{}/T`$). This involves the introduction of a cut-off $`\mathrm{\Lambda }`$ in the momentum of the order of the temperature $`\mathrm{\Lambda }T/\mathrm{}`$. The resulting cut-off dependences reflect directly the divergences of the classical theory and indicate its different high-momentum behavior with respect to the quantum theory. One might hope that for low-momentum correlation functions this different behavior does not play a role so that the physics involved at the cut-off scale $`T/\mathrm{}`$ is unimportant. However, this is not so as the high-momentum modes affect, through interactions, also the low-momentum sector of the theory in an essential way .
A general strategy to improve the classical theory is to include the dominant quantum contributions from the high-momentum (hard) modes. In this context it is important to understand the classical divergences, since they correspond to the dominant hard-mode contributions in the quantum theory. For instance, the linear divergences in classical non-abelian gauge theories have a one-to-one correspondence to the well-known hard thermal loops (HTLs) . The fact that hard thermal loops have to be included in an effective theory for the soft (low-momentum) modes has been known since the work of Braaten and Pisarski . These hard thermal loops have the following remarkable properties, namely, they are gauge invariant, they satisfy abelian-like Ward identities  and they allow a kinetic formulation in terms of a Vlasov equation .
Besides linear divergences, classical theories contain also logarithmic divergences starting at two-loop. (The explanation of how loops arise in a classical theory can be found in .) Higher-loop diagrams are superficially finite, although they may contain linear or logarithmic divergences in the form of one- or two-loop subdiagrams, respectively. Thus, we have two essential kinds of classical divergences, linear and logarithmic.
The relevance and the physical significance of the linear divergences (the HTLs) has been extensively studied, in particular, an effective action which incorporates them has been developed (see  or  for example).
A further step towards an effective theory beyond hard thermal loops would then be to include the logarithmic divergences (“log-divergences”) into the effective action. At this point some natural questions arise, namely: Do the log-divergences have the same or similar properties as the linear divergences, namely the hard thermal loops? What are the quantum contributions corresponding to the log divergences?
In any case, an extra term in the effective action containing the physics of the logarithmic divergences would enter the equations of motion of the fields as a current, which must be conserved for the effective theory to be consistent. In this letter we show that the logarithmically divergent part of the classical $`SU(N)`$ gauge self-energy is transverse, therefore providing current conservation at the level of two-point functions in the equations of motion. For that we will use the Ward identities for the self-energy at finite temperature and the results of , where it was argued that the degree of divergence $`d`$ of a classical diagram with $`\mathrm{}`$ loops is $`d=2\mathrm{}`$ and also that there are no log divergences at one-loop.

## II Transversality of Logarithmic divergences

Gauge invariance provides us with identities among different Green functions, the so-called Ward identities. In $`SU(N)`$ gauge theories the Ward identitiy for the full retarded propagator $`D_{\mu \nu }^{\text{R, full}}(K)`$ in the covariant gauge ($`_\mu A^\mu =0`$) is

$$K^\mu K^\nu D_{\mu \nu }^{\text{R, full}}=\alpha $$
(1)
and in general, in linear gauges ($`F_\mu A^\mu =0`$) it is given by

$$F^\mu F^\nu D_{\mu \nu }^{\text{R, full}}=\alpha ,$$
(2)
where $`\alpha `$ is the gauge parameter in both cases. The identity is the same at $`T=0`$ and $`T0`$
The identities can be written for the self energy $`\mathrm{\Pi }_{\mu \nu }`$ using the relation $`(D_{\text{full}})^1=(D_0)^1\mathrm{\Pi }`$. Now, a difference arises for the two cases $`T=0`$ and $`T0`$. This difference can be better seen in the covariant gauge. At zero temperature, $`\mathrm{\Pi }_{\mu \nu }`$ must be a linear combination of the two available tensors $`g_{\mu \nu }`$ and $`K_\mu K_\nu `$, which immediately leads to the result that the self-energy is transverse, namely $`K^\mu \mathrm{\Pi }_{\mu \nu }=0`$. However, at finite temperature the difference arises due to the presence of a heat bath with four-velocity $`U^\mu `$. From (2) and the presence of $`U^\mu `$ the Lorentz structure of $`D_{\text{R, full}}^{\mu \nu }`$ allows different independent tensors combination of both $`K_\mu `$ and $`U_\mu `$, for instance $`g_{\mu \nu }`$, $`K_\mu K_\nu `$, $`U_\mu U_\nu `$ and $`K_\mu U_\nu +U_\mu K_\nu `$. More convenient are the dimensionless tensors $`A_{\mu \nu }`$, $`B_{\mu \nu }`$, $`C_{\mu \nu }`$ and $`D_{\mu \nu }`$ detailed in . Here we will use the fact that both $`A_{\mu \nu }`$ and $`B_{\mu \nu }`$ are transverse to the four-momentum $`K^\mu `$, i.e. $`K^\mu A_{\mu \nu }=0`$ and $`K^\mu B_{\mu \nu }=0`$, whereas $`C_{\mu \nu }`$ and $`D_{\mu \nu }`$ are not. Moreover, $`A_{\mu \nu }`$ and $`B_{\mu \nu }`$ are respectively transverse and longitudinal to the spatial momentum $`𝐤`$. With these tensors the self-energy can be written as

$$\mathrm{\Pi }_{\mu \nu }=\mathrm{\Pi }_TA_{\mu \nu }+\mathrm{\Pi }_LB_{\mu \nu }+\mathrm{\Pi }_CC_{\mu \nu }+\mathrm{\Pi }_DD_{\mu \nu },$$
(3)
where $`\mathrm{\Pi }_T`$ and $`\mathrm{\Pi }_L`$ denote the four-momentum transverse components (transverse and longitudinal to the spatial momentum respectively) and $`\mathrm{\Pi }_C`$ and $`\mathrm{\Pi }_D`$ are the non-transverse components. This decomposition of the self-energy in the above basis of tensors is not only valid for the covariant gauge, but also for the temporal and Coulomb gauges and, in general, for linear gauges that do not break rotational invariance .
From the Ward identity (2) and the decomposition (3) above one can derive a relation between the different transverse and non-transverse components of the self-energy :

$$[\mathrm{\Pi }_C(K)]^2=[K^2+\mathrm{\Pi }_L(K)]\mathrm{\Pi }_D(K).$$
(4)
We note that this result does not imply that the self-energy is transverse. Indeed, it is well known that already at one-loop the self-energy is not transverse . However, remarkably, the hard thermal loop part of the self-energy is, i.e. $`K^\mu \mathrm{\Pi }_{\mu \nu }^{\text{HTL}}=0`$. This is due to the fact that HTLs satisfy abelian-like ward identities .
The Ward identity (4) will be a starting point in our discussion on the transversality of the divergent parts of the classical self-energy.

### A Linear divergences

Let us now consider the classical approximation of $`SU(N)`$ gauge theory, which is obtained by taking the $`\mathrm{}0`$ limit of the quantum theory. The classical theory is expected to be a good approximation at low energies because the classical and low-energy limit of the Bose-Einstein distribution function $`n(\omega _𝐤)`$ yield the same result:

$$n(\omega _𝐤)=\frac{1}{e^{\mathrm{}\omega _𝐤/T}1}n_{\text{cl}}(\omega _𝐤)\frac{T}{\omega _𝐤\mathrm{}},$$
(5)
where $`\omega _𝐤=|𝐤|`$ is the frequency at wavenumber $`𝐤`$. As mentioned in the introduction, the classical theory as an effective theory requires the introduction of a cut-off $`\mathrm{\Lambda }`$ which appears in the calculation of diagrams, reflecting directly the divergences of the theory.

The linearly divergent terms in the classical theory correspond to the HTLs in the quantum theory . This is so because the HTLs are proportional to the $`\omega _{\text{pl}}^2`$, where $`\omega _{\text{pl}}gT\mathrm{}^{1/2}`$ is the plasmon frequency. Thus, the HTLs behave as $`1/\mathrm{}`$, which become the linear divergences in the classical theory as we take $`\mathrm{}0`$.

The fact that HTLs are transverse indicates that the linearly divergent terms should also be so. This can be checked by making use of (4). Consider the case $`K^2=K_\mu K^\mu 0`$ and let us start at one-loop. Since $`\mathrm{\Pi }_C=0`$ at tree level, it begins at order $`O(g^2)`$ (one-loop). Now, from (4) we notice that $`\mathrm{\Pi }_D`$ should start at $`O(g^4)`$ (two-loops). The two-loop contribution $`\mathrm{\Pi }_D^{[2]}`$ (the superscript denotes loop order) is superficially log divergent, containing at most a linear subdivergence. Hence by (4) we see that the one-loop contribution $`\mathrm{\Pi }_C`$ cannot have a linear divergence. Thus, at one-loop both $`\mathrm{\Pi }_C`$ and $`\mathrm{\Pi }_D`$ vanish, and therefore,

$$K^\mu \mathrm{\Pi }_{\mu \nu }^{\text{[1], lin}}=0,$$
(6)
as we expected from our considerations of hard thermal loops above. In fact, this is another way of showing that HTLs are transverse.

### B Logarithmic divergences

At one-loop there are no logarithmic divergences , so we consider the case of two-loop, i.e. $`O(g^4)`$. We first split the two-loop self-energy component $`\mathrm{\Pi }_D^{\text{[2]}}`$ in a logarithmically divergent part, a part that may contain a linear subdivergence and a finite part:

$$\mathrm{\Pi }_D^{\text{[2]}}=\mathrm{\Pi }_D^{\text{[2], log}}+\mathrm{\Pi }_D^{\text{[2], sublin}}+\mathrm{\Pi }_D^{\text{[2], fin}}.$$
(7)
We insert this expression into the Ward identity (4). The terms in $`\mathrm{\Pi }_C`$ at the right-hand side that match with the $`O(g^4)`$ at the left-hand side are those corresponding to one-loop, which does not have logarithmic divergences, and therefore

$$\mathrm{\Pi }_D^{\text{[2], log}}=0.$$
(8)
We saw already that $`\mathrm{\Pi }_C`$ does not contain a linear divergence, thus also

$$\mathrm{\Pi }_D^{\text{[2], sublin}}=0.$$
(9)
Next, we consider $`\mathrm{\Pi }_C`$. Analogously to (7), we split it in terms of the different types of divergences

$$\mathrm{\Pi }_C^{\text{[2]}}=\mathrm{\Pi }_C^{\text{[2], log}}+\mathrm{\Pi }_C^{\text{[2], sublin}}+\mathrm{\Pi }_C^{\text{[2], fin}}.$$
(10)
We use the Ward identity (4) at $`O(g^8)`$, which we may write as

$`\left(\mathrm{\Pi }_C^{\text{[2], log}}+\mathrm{\Pi }_C^{\text{[2], sublin}}+\mathrm{\Pi }_C^{\text{[2], fin}}\right)^2+2\mathrm{\Pi }_C^{\text{[1], fin}}\mathrm{\Pi }_C^{\text{[3]}}`$ (11)
$`=K^2\mathrm{\Pi }_D^{\text{[4]}}+\mathrm{\Pi }_L^{\text{[1]}}\mathrm{\Pi }_D^{\text{[3]}}+\mathrm{\Pi }_L^{\text{[2]}}\mathrm{\Pi }_D^{\text{[2]}}.`$ (12)
We now focus on the terms that could lead to a logarithmic divergence. We keep from (12) terms proportional to $`(\mathrm{log}\mathrm{\Lambda })^2`$. This results in

$$\left(\mathrm{\Pi }_C^{\text{[2], log}}\right)^2+2\mathrm{\Pi }_C^{\text{[1], fin}}\mathrm{\Pi }_C^{\text{[3]}}=K^2\mathrm{\Pi }_D^{\text{[4]}}+\mathrm{\Pi }_L^{\text{[1]}}\mathrm{\Pi }_D^{\text{[3]}}.$$
(13)
As a consequence of (8), $`\mathrm{\Pi }_D^{\text{[2]}}`$ does not contain a logdivergence, therefore the last term on the r.h.s. of (12) does not contribute to (13). Let us consider the products of one- and three-loop contributions. Since at one-loop there are no log-divergences, the three-loop diagrams must contain a double log-divergence for these products to contribute. Now, schematically, the expression for a three-loop diagram is

$$\mathrm{\Pi }^{\text{[3]}}(P)=g^6T^3\frac{d^4K}{(2\pi )^4}\frac{d^4K^{}}{(2\pi )^4}\frac{d^4K^{\prime \prime }}{(2\pi )^4}f^{\text{[3]}}(K,K^{},K^{\prime \prime },P),$$
(14)
where $`K`$, $`K^{}`$ and $`K^{\prime \prime }`$ are the internal momenta, $`P`$ is the external momentum, $`g`$ the coupling constant and $`T`$ the temperature. The result after the integration over any two arbitrary internal momenta can be regarded as either an expression for two disjunct one-loop diagrams or a two-loop diagram, with external lines depending on the other momenta. Consider for example the three-loop diagram in figure 1.

In the case that the integration over two arbitrary internal momenta ($`K^{}`$ and $`K^{\prime \prime }`$ in Fig. 1.) corresponds to a two-loop diagram (as in Fig. 1b.), it can at most give a single logarithmic divergence, which we denote as log$`\mathrm{\Lambda }`$. When it does, the integration over the remaining momentum ($`K`$ in Fig. 1) cannot give an extra log$`\mathrm{\Lambda }`$, since the superficial degree of divergence of the total diagram is $`d=2\mathrm{}=1`$. In the case that the integration over $`K^{}`$ and $`K^{\prime \prime }`$ does not give a logdivergence, the integration over $`K`$ could still lead to one log$`\mathrm{\Lambda }`$. Hence, a three-loop diagram can at most give a single log-divergence. Therefore, the product of one- and three-loop diagrams cannot contribute to (13).
The above argument can be repeated for the four-loop contribution to the self-energy. In this case there are four internal momenta. The result after integration over any three given internal momenta can be regarded as the expression for a disjunct two and one-loop diagram or three disjunct one-loop diagrams. Therefore it can at most give a single logarithmic divergence, and since a four-loop contribution to the self-energy is finite, the remaining integration cannot give an extra log-divergence and as in the case above, $`\mathrm{\Pi }^{\text{[4]}}`$ cannot contribute to (13). Thus, we find from (13) that

$$\mathrm{\Pi }_C^{\text{[2], log}}=0.$$
(15)
Note that we cannot say that the two-loop contribution to $`\mathrm{\Pi }_C`$ containing a linear divergence from one-loop subdiagrams equals zero, as we could for $`\mathrm{\Pi }_D`$.
Since both $`\mathrm{\Pi }_C^{\text{[2]}}`$ and $`\mathrm{\Pi }_D^{\text{[2]}}`$ vanish , then we conclude that the logarithmically divergent part of the two-loop classical self-energy is transverse

$$K^\mu \mathrm{\Pi }_{\mu \nu }^{\text{[2], log}}=0.$$
(16)
Our argument does not hold for three- or higher-loop log-divergences. However, since those diagrams are superficially finite, the divergences can only come through two-loop subdiagrams and in general, we do not expect the whole contribution to the self-energy to be transverse. The important point is that, as we mentioned in the introduction, the essential divergences that appear in the classical theory are linear (at one-loop) and logarithmic (at two-loop) and they are both transverse.

## III Conclusions and outlook

Here we showed that the logarithmically divergent parts of the classical finite temperature $`SU(N)`$ self-energy are transverse. This also holds for the hard thermal loops, which correspond to classical linear divergences. Therefore we see that all divergences appearing in the classical self-energy are transverse. We would like to comment on the importance of this result regarding the construction of a low-energy effective theory beyond the HTL approximation.

Consider for instance the effective action which results from integrating out the hard modes with momenta $`P>\mu `$, with $`\mu `$ an intermediate scale such that $`\omega _{\mathrm{pl}}<\mu <T/\mathrm{}`$. In a $`\mathrm{}`$ or high $`T`$ expansion the effective action for the soft modes would schematically be written as

$`\mathrm{\Gamma }_{\text{eff}}=g^2T\left({\displaystyle \frac{T}{\mathrm{}}}c_1\mu \right)\overline{\mathrm{\Gamma }}_{\text{HTL}}+\left(g^2T\right)^2\text{log}\left({\displaystyle \frac{c_2T}{\mathrm{}\mu }}\right)\overline{\mathrm{\Gamma }}_{\text{log}}`$ (17)
$`+S_{\text{cl}}+O(g^2\mathrm{},{\displaystyle \frac{\omega _{\text{pl}}}{\mu }},{\displaystyle \frac{\mathrm{}\mu }{T}}),`$ (18)
where $`c_1`$ and $`c_2`$ are constants that depend on the regularization scheme. The first term in the expansion, which corresponds to the HTLs, is proportional to $`\mathrm{}^1`$, being thus linearly divergent in the classical limit $`\mathrm{}0`$. The second term is proportional to log ($`T/\mathrm{}`$) and so corresponds to the logarithmic divergences in the classical theory. The third term is the classical action and the other terms in the expansion are unimportant contributions in either a high $`T`$ or classical regime. Thus we see that in a $`\mathrm{}`$ or high-temperature expansion the next-to-leading order terms are given by the classical log-divergences.

A consistent scheme to include hard-mode contributions beyond HTLs in the classical theory seems to require the inclusion of terms that diverge in the classical limit. An effective action of the form (18) gives rise to currents in the equations of motion for the classical $`SU(N)`$ gauge field

$$\delta _{A_\mu }S_{\text{cl}}=j_{\text{HTL}}^\mu +j_{\text{log}}^\mu ,$$
(19)
where $`j_{\text{HTL}}^\mu \delta _{A_\mu }\overline{\mathrm{\Gamma }}_{\text{HTL}}`$ and $`j_{\text{log}}^\mu \delta _{A_\mu }\overline{\mathrm{\Gamma }}_{\text{log}}`$. The current $`j_{\text{HTL}}^\mu `$ generated by the HTLs is conserved. For consistency, it is necessary that the current $`j_{\text{log}}^\mu `$ generated by the log-divergences is also conserved. Our result (16) shows that the logarithmic divergent part of the self-energy is transverse, which implies that $`j_{\text{log}}^\mu `$ is indeed conserved. We stress that this is a special property at finite temperature that should not generally be expected, and in fact, this result encourages the study of the classical logarithmic divergences towards the construction of a feasible low-energy effective theory beyond hard thermal loops.


format: text

text:
# Moving Wigner Glasses and Smectics: Dynamics of Disordered Wigner Crystals

\[

## Abstract

We examine the dynamics of driven classical Wigner solids interacting with quenched disorder from charged impurities. For strong disorder, the initial motion is plastic—in the form of crossing winding channels. For increasing drive, the disordered Wigner glass can reorder to a moving Wigner smectic–with the electrons moving in non-crossing 1D channels. These different dynamic phases can be related to the conduction noise and $`I(V)`$ curves. For strong disorder, we show criticality in the voltage onset just above depinning. We also obtain the dynamic phase diagram for driven Wigner solids and prove that there is a finite threshold for transverse sliding, recently found experimentally.

\]

Ordered arrays of charged particles have been studied in the context of colloidal suspensions, ion rings, atomic-ion Wigner crystals, quantum computers, astrophysics, biophysics, plasmas, electrons deposited on liquid Helium surfaces, semiconductor heterostructures, Wigner droplets, and arrays of metallic islands interconnected by tunnel junctions . A revival of general interest in charged arrays has been fueled by the observation, in 2D heterostructures, of nonlinear I-V curves exhibiting a threshold as a function of an externally applied electric field , indicating the presence of a Wigner solid (WS)  that has been pinned by the disorder in the sample. The depinning threshold can vary up to two orders of magnitude in different samples . Indeed, experiments  based on transport and photoluminescence provide indirect evidence (without observing the location of the charges) for the existence of the WS, and demonstrate the very important role disorder plays in the dynamics of the WS . It is the purpose of this work to study how disorder affects the transport properties of the driven WS.

The onset of broadband conduction noise has been interpreted as a signature of the sliding of a defected WS . If the electrons were to retain their order and slide collectively, narrow band noise resembling that seen in sliding charge-density wave systems should be observed instead. Simulations by Cha and Fertig on classical Wigner crystals interacting with charged defects indicate that a disorder-induced transition, from a clean to a defected WS, can occur as a function of increasing pinning strength . For strong pinning, the initial depinning is plastic and involves tearing of the electron crystal (b). However, many aspects of this transport have not been systematically characterized, including the current-voltage $`I(V)`$ curves, conduction noise, transverse meandering, transverse threshold for sliding, and the electron lattice structure and motion for varying applied drives.

Similar plastic depinning transitions have been observed in the related system of driven vortex lattices in disordered superconductors. The plastic depinning of the vortices is associated with flux motion through intricate river-like channels . Defects in the vortex lattice strongly affect the depinning thresholds and the voltage noise signatures produced by the system. For increasing drive the initially defected vortex lattice can reorder to a moving lattice or a moving smectic phase . In the fast moving lattice phase the vortex lattice regains order in both transverse and longitudinal directions with respect to the driving force, while in the smectic state only order transverse to the driving direction appears.

It is unclear a priori whether the same type of reordering transitions can occur in a system containing long-range pinning, such as in the case of the WS interacting with charged impurities. This is in contrast with the vortex system, where the pinning interacts with the flux lines only on a very short length scale. Additionally, in the strong pinning limit, critical behavior may occur near the depinning threshold, leading to velocity-force relations  of the form $`v(FF_T)^\xi `$. For instance, for transport in metallic dots, $`\xi =5/3`$ theoretically , and $`\xi =2.0`$ and $`1.58`$ experimentally .

In this work, we use numerical simulations to study dynamical transitions in a 2D electron system forming a classical WS in the presence of charged impurities. For strong pinning, the electron crystal is highly defected and depins plastically, with certain electrons flowing in well-defined channels or interconnecting rivers, while others remain immobile. For increasing applied driving force, the initially disordered electrons can partially reorder to a moving Wigner smectic state where the electrons flow in 1D non-crossing winding channels. The reordering transition is accompanied by a saturation in the $`dI/dV`$ curves and by a change from a broad-band to a narrow-band voltage noise signature. For weak disorder, the depinning is elastic and a narrow-band noise signal can be observed at all drives above depinning. For strong disorder, where the depinning is plastic, we find criticality in the velocity force curves in agreement with transport in metallic dots . We map out the dynamic phase diagram as a function of disorder strength and applied driving force. Also, we find unequivocal evidence for a transverse finite threshold for sliding conduction, in agreement with recent experiments . Our results can also be tested in other systems listed in the first sentence of this paper (e.g., with charged colloidal systems suspended in water in the presence of a disordered substrate).

We conduct overdamped molecular dynamics (MD) simulations using a model similar to that studied by Cha and Fertig . The energy from interactions is

$$U=\underset{ij}{}\frac{e^2}{|𝐫_i𝐫_j|}\underset{ij}{}\frac{e^2}{\sqrt{(|𝐫_i𝐫_j^{(p)}|^2+d^2)}}$$
(1)
The first term is the electron-electron (Coulomb) repulsion and the second term is the electron-impurity interaction, where the impurities are positively charged defects out of plane. Here, $`𝐫_i`$ is the location of electron $`i`$, and $`𝐫_j^{(p)}`$ is the in-plane location of a positive impurity which is located at an out-of-plane distance $`d`$ (measured here in units of $`a_0`$, the average lattice constant of the WS). The number of electrons $`N_i`$ equals the number of impurities $`N_p`$, and the disorder strength is varied by changing $`d`$. We have also considered many other cases (e.g., $`\pm `$ charged impurities, $`N_iN_p`$, etc.) with consistent results. The long-range Coulomb interactions are evaluated with a fast converging sum technique by Jensen  which is computationally more efficient than the Ewald-sum technique used in previous studies .

To obtain the initial electron positions, we perform simulated annealing in which we start from a high temperature above the electron lattice melting transition and slowly cool to a low temperature. Temperature is implemented via Langevin random kicks. Once the electron configuration has been initialized, the critical depinning force is determined by applying a very slowly increasing uniform driving force $`f_d`$ which would correspond to an applied electric field. After each drive increment, we wait $`10^4`$ MD time steps before taking data, which we average over the next $`10^4`$ time steps. For each drive increment, we measure the average electron velocity ($``$ current) in the direction of drive, $`V_x=(1/N_i)_i^{N_i}\widehat{𝐱}𝐯_i`$. The $`V_x`$ versus $`f_d`$ curve corresponds to an $`I(V)`$ experimental curve and we will thus use the notation $`V_x=I`$ and $`f_d=V`$. The depinning force $`f_d^c`$ is defined as the drive $`f_d`$ at which $`V_x`$ reaches a value of $`0.01`$ that of the ohmic response. We have studied system sizes $`N_i`$ from $`64`$ to $`800`$ and find similar behavior at all sizes. Most of the results presented here are for systems with $`N_i=256`$.

When driven through a sample containing strong pinning, the electron lattice undergoes a gradual reordering transition as the driving force is increased. We illustrate such a reordering transition in Fig. 1 by plotting the eastbound electron trajectories at different driving forces for a sample containing strong disorder of $`d=0.65`$. In Fig. 1(a) the onset of motion occurs through the opening of a single winding channel. Electrons outside of the channel remain pinned, and the overall electron lattice is disordered. At $`f_d/f_d^c=1.5`$, shown in Fig. 1(b), several channels have opened, some of which are interconnecting.

The original channel in Fig. 1(a) has grown in width, but regions of pinned electrons are still present. Electrons moving past a pinned electron perturb it, causing it to move (like a revolving door) in a circular orbit around the center of the potential minima in which it is trapped. Several of these electron turnstiles can be observed in Fig. 1(b). In Fig. 1(c) for $`f_d/f_d^c=2.25`$ there appear to be regions where electrons do not flow, but none of the electrons are permanently pinned. Some electrons become pinned for a period of time before moving again. If the trajectories are drawn for a sufficiently long time, the electron flow appears everywhere in the sample, although there are still preferred paths in which more electrons flow. In Fig. 1(d) for $`f_d/f_d^c=3.0`$ the electron flow occurs more uniformly across the sample. In Fig. 1(e) at $`f_d/f_d^c=4.0`$ the electrons begin to flow predominantly in certain non-crossing channels, although some electrons jump from channel to channel. In Fig. 1(f) for $`f_d/f_d^c=5.0`$ the electron flow occurs in well defined non-

crossing channels, which can contain different numbers of electrons. A similar channel motion exists for driven vortices in disordered superconductors .

For samples containing very weak disorder, the pinned WS has six-fold ordering and depins elastically, without generating defects in the lattice. In this regime, the electron crystal flows in 1D channels with each channel containing the same number of electrons. Here, the transverse wanderings of the electrons are considerably reduced compared to the case of strong pinning.

To better illustrate the change in the amount of disorder in the electron lattice, in Fig. 2(a) we show the Delaunay triangulation for the electrons for the same drive as in Fig. 1(f). Defects, in the form of 5-7 disclination pairs, appear with their Burgers vectors oriented perpendicular to the direction of the drive. In the structure factor there are only two prominent peaks for order in the direction transverse to the drive, consistent with a moving Wigner smectic state. In Fig. 2(b) we show the Delaunay triangulation for the moving state in a weakly pinned sample with $`d=1.76`$ where the initial depinning is elastic. Here the moving lattice is defect free. The structure factor in this case shows four longitudinal peaks in addition to the two transverse peaks, although the transverse peaks are more prominent. Much larger systems would be necessary to determine whether the system is in a smectic state or in a moving Bragg glass in this case of weak disorder.

In order to connect the reordering sequence with a measure that is readily accessible experimentally, we plot in Fig. 3(a) the $`I(V)`$ and $`dI/dV`$ curves. A peak occurs in $`dI/dV`$ at $`f_d/f_d^c3`$, when the electrons are undergoing very disordered plastic flow. In Fig. 3(b) we plot the fraction $`P_6`$ of six-fold coordinated electrons as a function of drive. A perfect triangular lattice would have $`P_6=1`$. For drives $`f_d/f_d^c<2`$ the lattice is highly defected. For $`f_d/f_d^c>3`$ (i.e., past the peak in $`dI/dV`$), the order in the lattice begins to increase. The value of $`P_6`$ saturates near $`f_d/f_d^c6`$ which also coincides with the saturation of the $`dI/dV`$ curve. Thus the experimentally observable $`I(V)`$ characteristic can be considered a good measure of

the degree of order and the nature of the flow in the system. To quantify the degree of plasticity of the electron flow, we plot in Fig. 3(c) the fraction of electrons $`D_{tr}`$ that wandered a distance of more than $`a_0/2`$ in the direction transverse to the drive during an interval of 8000 MD steps. The peak in $`D_{tr}`$ coincides with the peak in $`dI/dV`$. $`D_{tr}`$ then slowly declines until it saturates at $`f_d/f_d^c6`$, indicating the gradual formation of the non-crossing channels as seen in Fig. 1(f). The saturation in $`D_{tr}`$ coincides with the saturations in both $`P_6`$ and $`dI/dV`$.

The $`I(V)`$ curve in Fig. 3(a) corresponds to a highly irregular voltage signal as a function of time when the electrons are in the plastic flow regime ($`f_d/f_c=1.5`$). The corresponding voltage noise spectrum in the inset of Fig. 3(c) shows that only broad-band noise is present. In contrast, at $`f_d/f_d^c=4.0`$, in the moving smectic regime, a roughly regular signal is obtained, and narrow-band noise appears, as shown in the inset of Fig. 3(c). In a system with $`d=1.57`$ when the depinning is elastic only, an even more pronounced narrow band noise signal is observable. In experiments, broad band noise has been observed above depinning , but narrow band noise has not been seen. Our results suggest that plasticity may be playing an important role in most experiments.

In Fig. 3(d,e) we examine the critical behavior above the depinning threshold of $`V_x`$ versus $`f_d`$ for disorder strengths $`d=0.5`$ and $`0.65`$. In Fig. 3(d), a log-log plot of $`V_x[(f_df_d^c)/f_d^c]^\xi `$, indicates that the curves are fit well by a power law over one and half decades with $`\xi =1.61\pm 0.10`$ and $`1.71\pm 0.10`$, respectively. These values agree well with the predicted value of $`\xi =5/3`$  for electron flow through disordered arrays, and the experimentally observed values of $`\xi =1.58`$ and $`2.0`$ .

The theory of the moving ordered phase  predicts a barrier to transverse motion once channels similar to those in Fig. 1(e) form. As shown in the inset of Fig. 4, for a system in the reordered phase we observe a transverse depinning threshold that is about 1/6 the size of the longitudinal depinning threshold $`f_d^c`$. In recent experiments, Perruchot et al.  also find evidence for a transverse barrier that is about 1/10 the size of the longitudinal threshold. These thresholds are much larger than those observed in vortex matter interacting with short-ranged pinning, where ratios of 1/100 are seen .

In Fig. 4 we present the WS dynamic phase diagram as a function of disorder strength and driving force. For $`d<1.4`$ there are a considerable number of defects in the WS, and the initial depinning is plastic. We label this region the pinned Wigner glass. For increasing disorder, the pinned region grows while for strong enough drive the electrons reorder to a moving Wigner smectic state. For $`d>1.4`$ there are few or no defects in the pinned state and the initial depinning is elastic. The presence of a crossover from elastic to plastic depinning with decreasing disorder strength is in agreement with .

In summary, we have investigated the pinning and dynamics of an electron solid interacting with charged disorder. We find that for strong disorder the depinning transition is plastic with electrons flowing in a network of winding channels. For increasing drives, the electrons partially reorder and flow in non-crossing channels forming a moving Wigner smectic. We show that the onset of these different phases can be inferred from the transport characteristics. In the plastic flow regime, the noise has broad-band characteristics, while in the moving smectic or elastic flow phase, a narrow band noise signal is observable. We also show that the onset of the plastic flow phase shows critical behavior with critical exponents in agreement with predictions for transport in arrays of metallic dots . We map out the dynamic phase diagram as a function of disorder and applied driving force. We obtain a finite threshold for transverse sliding, in agreement with recent experiments .

We gratefully acknowledge J. Groth, Niels Grønbech-Jensen, and Alan H. MacDonald for useful discussions, and B. Janko for his kind hospitality. This work was partially supported by DOE Office of Science No. W-31-109-ENG-38, CLC and CULAR (LANL/UC), and NSF DMR-9985978.


text:
# Prediction of Long Term Stability by Extrapolation

## Chapter 1 Introduction

This paper studies the possibility of using the survival function, to predict long term stability by extrapolation . The survival function is a function of the initial coordinates and is the number of turns a particle will survive for a given set of initial coordinates. To determine the difficulties in extrapolating the survival function, tracking studies were done to compute the survival function. The survival function was found to have two properties that may cause difficulties in extrapolating the survival turns function. One is the existence of rapid oscillations, and the second is the existence of plateaus.It was found that it appears possible to extrapolate the survival function to estimate long term stability by taking the two difficulties into account.

A model is proposed which pictures the survival function to be a series of plateaus with rapid oscillations superimposed on the plateaus. The tracking studies give results for the widths of these plateaus and for the seperation between adjacent plateaus which can be used to extrapolate and estimate the location of plateaus that indicate survival for longer times than can be found by tracking.

## Chapter 2 The survival function, $`j_{trns}`$.

For a given set of initial coordinate, $`x_0,p_{x0},y_0,p_{y0}`$, one can find by tracking the survival time in turns, which is the number of turns the particle will survive before becoming unstable and which will be denoted by $`j_{trns}`$. This determines the function $`j_{trns}(x_0,p_{x0},y_0,p_{y0})`$ which will be called the survival function . If one limits the tracking to $`\mathrm{1.0\hspace{0.33em}10}^6`$ turns or less, one can find $`j_{trns}`$ for those $`x_0,p_{x0},y_0,p_{y0}`$ for which $`j_{trns}`$ is less than or equal to $`\mathrm{1.0\hspace{0.33em}10}^6`$.

A tracking study was done of particle motion with no rf present, using an older version of the RHIC lattice. Random and systematic field error multipoles are present up to order 10. The particle momemtum is $`dp/p=0`$. As the first case studied, the initial coordinates $`x_0,p_{x0},y_0,p_{y0}`$ are chosen along the direction in phase space given by $`p_{x0}=0,p_{y0}=0`$ and $`ϵ_{x0}=ϵ_{y0}`$, where $`ϵ_{x0},ϵ_{y0}`$ are the linear emittances in the absence of the error multipoles. Along this direction in phase space , $`j_{trns}`$ may be considered to be a function of $`x_0`$. For a given initial coordinate, $`x_0`$, one can find by tracking the survival time in turns, which is the number of turns the particle will survive before becoming unstable and will be denoted by $`j_{trns}`$. This determines the survival function  $`j_{trns}(x_0)`$. If one limits the tracking to $`\mathrm{1\hspace{0.33em}10}^6`$ turns or less, one can find $`j_{trns}`$ for those $`x_0`$ for which $`j_{trns}`$ is less than or equal to $`\mathrm{1\hspace{0.33em}10}^6`$. For $`dp/p=0`$, it is assumed that on the closed orbit $`x=0`$ and $`x_0`$ is also the initial betatron oscillation amplitude.

The tracking was first done using $`x_0`$ which are seperated by $`dx_0=.1`$ mm. The results for $`j_{trns}`$ versus $`x_0`$ are shown in Fig. 2.1. The results in this figure may be looked at as the results of a search in $`x_0`$ starting at large $`x_0`$ and decreasing $`x_0`$ in steps of $`dx_0=.1mm`$ . The figure shows an apparent stability limit for 1e6 turns of $`u_{sl}=`$15.2 mm.

### 2.1 Rapid oscillations in the survival function

Fig. 2.1 shows rapid oscillations in $`j_{trns}`$ with $`x_0`$. Large changes in $`j_{trns}`$ occur when $`x_0`$ changes by .1 mm. The oscillations extend over about .3 mm. Tracking studies show that the oscillations become more rapid when the search interval, $`dx_0`$, is decreased. This is indicated in Fig. 2.2 where $`dx_0`$ is decreased to .05 mm and in Fig. 2.3 where $`dx_0`$ is decreased to .025 mm. Results for different search intervals, $`dx_0`$, are shown in Table 2.1. $`dx_0`$ is decreased from .1mm to .0001 mm . The wavelength of the oscilations , $`\mathrm{\Delta }x_0`$, as measured from the stability limit for 1e6 turns, $`u_{sl}`$, to the location of the first peak in $`j_{trns}`$ decreases from about .3 mm to .0002 mm. Also listed for each value of $`dx_0`$ are the apparent stability limit for 1e6 turns, $`u_{sl}`$, and $`j_{trns}`$ at $`u_{sl}+dx_0`$.

Table 2.1 shows that the wavelength of the oscilations , $`\mathrm{\Delta }x_0`$ is roughly proportional to the size of the search interval, $`dx_0`$. The value of $`j_{trns}`$ at $`u_{sl}+dx_0`$ shows that near $`u_{sl}`$, $`j_{trns}`$ changes appreciably in the small change in $`x_0`$ given by $`dx_0`$. Tracking results show that this seems to hold even at extremely small $`dx_0`$.The computed results appear to indicate that $`j_{trns}(x_0)`$ is not a continuous function of $`x_0`$. For a continuous function of $`x_0`$, one can find a small enough interval in $`x_0`$ such that the differnce between the values of the function for any two $`x_0`$ in that interval is very small. This does not appear to be true for $`j_{trns}(x_0)`$.

The existence of the rapid oscillations in the survival function, $`j_{trns}(x_0)`$, would seem to make it difficult to extrapolate to find those $`x_0`$ that survive for more than 1e6 turns. However one could view the survival function shown in Fig. 2.1 as being made up of the rapid oscillations superimposed on a smoother, more slowly varying function which could be used for the extrapolation. This is discussed further in section 2.3.

### 2.2 Plateaus in the survival function

Looking at Fig. 2.1, one can make out plateaus in the survival function, $`j_{trns}(x_0)`$. The plateaus are regions where $`j_{trns}`$ oscillates rapidly around an almost constant value of $`j_{trns}`$. The plateaus can be seen somewhat more clearly if one reduces the search interval $`dx_0`$, as shown in Fig. 2.2 where $`dx_0`$ is decreased to $`dx_0=.05mm`$. One can make out 4 plateaus located at about $`j_{trns}`$=1.5e5, 2e4, 3500, 400 turns. Possible plateaus with $`j_{trns}`$ less than 100 turns are being ignored. It will be seen that the width of the plateaus do not depend on the search interval, $`dx_0`$. This is also true of the location of the plateaus in $`j_{trns}`$. This is shown in Fig. 2.3 where $`dx_0`$=.025 mm.

The existence of the plateaus in the survival function, $`j_{trns}(x_0)`$, would seem to make it difficult to extrapolate to find those $`x_0`$ that survive for more than 1e6 turns. If one does the extrapolation using points which are close to the stability limit for 1e6 turns, $`u_{sl}`$, one may get incorrect results as these points may lie on one of the plateaus.

### 2.3 Extrapolation of the survival function

The data given above leads to a model of the survival functioon, which pictures it as sequence of plateaus. Within the plateaus, $`j_{trns}`$ oscillates about some constant vaue of $`j_{trns}`$ which will be called the level of the plateau. The existence of the plateaus makes extrapolation of the survival function appear difficult. The last plateau that was measured has a level of about 1.5e5 turns and a width of about 1.2 mm. An interesting question is what is the level of the next plateau, at lower $`x_0`$, after the last plateau that was measured. In the following this plateau will be referred to as the ’next plateau’.

Some help with the extrapolation is also provided by plotting $`1/Log(j_{trns})`$ against $`x_0`$ as shown in Fig. 2.4, where the search interval is $`dx_0`$=.1mm and points with $`j_{trns}`$ less than 400 turns are omitted. It will be seen below that the seperation between adjacent plateaus does not vary greatly when measured as the change in $`1/Log(j_{trns})`$.

To help locate the ’next plateau’, long runs were done starting with $`x_0`$=15.2mm, and decreasing $`x_0`$ 1n steps of .1mm. In order to detect the beginning of the ’next plateau’, the runs have to be long enough not to be confused by the rapid oscillations in $`j_{trns}`$ that occur within each plateau. Runs of length 2e7 turns were chosen, and these runs take about 10 days for the RHIC lattice used.

The results are shown in Fig. 2.5 . One sees from this figure that the ’next plateau’ appears to start at $`x_0`$=15.3 mm and ends at about 14.0 mm, and has a level of about $`j_{trns}`$=15e6 turns and a width of 1.2mm . The previous plateau goes from 16.5mm to 15.2 mm and has a width of 1.3mm and a level of jtrns=1.5e5 turns. The width of the plateau is measured here from the end of one plateau in x0 to the end of the adjacent plateau and includes the transition region where the points move from one plateau to the next. The width of the ’next plateau’ is difficult to measure, as the adjacent plateau at still lower $`x_0`$ is estimated to have a level of 3e9 turns, and cannot be found by tracking. One can see that the width is larger than 1.2mm. The $`x_0`$ at 14.0 mm survived more than 89.9e6 turns. The data given above for the ’next plateau’ is somewhat in error but it seems better to use it than to throw away this information.

The results are also shown as a $`1/Log(j_{trns})`$ versus $`x_0`$ plot in Fig. 2.6. One can see three plateaus near the stability boundary at the levels of about jtrns= 15e6, 1.5e5,and 2e4 turns. The properties of the plateaus are summarized in Table 2.2. The seperation between the plateaus in the $`1/Log(j_{trns})`$ plot are given by .054 and .039, so that the seperation between the plateaus in the $`1/Log(j_{trns})`$ are not too different. The widths of the plateaus that were measured were 1.2+, 1.3 and 1.7 mm.

The plateau model will now be used to extrapolate and investigate long term stability in RHIC. In RHIC, $`j_{trns}`$=1e9 turns coresponds to a survival time of about 3.5 hours. The data given above will be used to extrapolate to find the plateau whose level is greater than or equal to $`j_{trns}`$=1e9 turns. This plateau will be called the 1e9 plateau. The existence of the plateaus indicates that there are limits on the accuracy one can hope to achieve by extrapolation. The two parameters one needs to extrapolate the survival function are the plateau width and the plateau level seperation. One cannot be certain what these parameters will be at $`j_{trns}`$ of the order of 1e9 turns. However, one can use the data found for these two parameters at $`j_{trns}`$ of the order of 1e6 turns , to make the best estimate for these parameters at $`j_{trns}`$ of the order of 1e9 turns.

The plateau width and the plateau level seperation were studied for three different cases coresponding to three different directions in phase space. The results are summarized in chapter 5. The widths of the plateaus were found to be roughly constant when measured in terms of $`X_0`$=$`[\beta _{x0}(ϵ_{x0}+ϵ_{y0})]^{.5}`$ with an average value $`\mathrm{\Delta }X_0`$=2.00 mm. The plateau level seperation measured in $`1/Log(j_{trns})`$ varied from .054 to .019 with an average value of .033. It is suggested that in the extrapolation, the average value of these two parameters be used. For the case being considered here, it is assumed in the extrapolation that the plateau widths in $`x_0`$ will be 1.4 mm corresponding to the average value of $`\mathrm{\Delta }X_0`$=2.00 mm, and the plateau level seperation is .033 in $`1/Log(j_{trns})`$. Note that in this case $`X_0`$=1.414 $`x_0`$

The results found using these assumptions are shown in Fig. 2.7, where two new plateaus are shown that were found by extrapolation. The plateau just below the next plateau in x0 has the level of $`j_{trns}`$=3.3e9 turns and is the 1e9 plateau in this example. The 1e9 plateau goes from $`x_0`$=13.9mm to $`x_0`$=12.6mm. The aperture for 1e9 turns may be taken as about 13.9mm, which may be compared to the 15.2 mm found using runs of 1e6 turns. This indicates a loss of about 1.3mm or about 9% due to the required survival time of 1e9 turns. One may note that the level of the 1e9 plateau being 3.3e9 turns, one may expect that some of the $`x_0`$ on this plateau will not survive 1e9 turns due to the oscillations in jtrns that will occur on this plateau.To be safer one could assume the aperture for 1e9 turns to be 12.6 mm, which is the beginning of the adjacent plateau with a level of 1.2e14 turns, giving a loss of 17%.

The procedure used indicates that in this case the result is not sensitive to changes in the two assumptions made regarding the width and level seperation of the plateaus. If the plateau model continues to hold for the extrapolated plateaus, and the width and level seperation remain roughly 1.3 mm and .034 respectively, then the result for the aperture for 1e9 turns will be about the same.

The plateau model developed above avoids certain problems that arise in tracking studies. In trying to find the aperture for the survival time of 1e9 turns, the more usual approach is to try find the $`x_0`$ such that all smaller $`x_0`$ will survive for 1e9 turns. If while doing a tracking search starting from large $`x_0`$, one finds a $`x_0`$ that survives for 1e9 turn, then one has to ask whether all smaller $`x_0`$ will survive for 1e9 turns. This is a difficult question to answer. It is even possible, that there are are always smaller $`x_0`$ that will not survive 1e9 turns, although these $`x_0`$ may become very scarce at smaller $`x_0`$. In the plateau model, in trying to find the aperture for the survival time of 1e9 turns, the approach is to try to find the plateau whose level is greater than 1e9 turns. This is a better defined target, and the expectation is that this plateau will indicate a region of $`x_0`$ where most $`x_0`$ will survive 1e9 turns.

## Chapter 3 Survival function along the $`ϵ_{y0}=0`$ direction

The tracking results presented above were all done along the $`ϵ_{x0}=ϵ_{y0}`$ direction. Results will now be given for the $`y_0=0,p_{y0}=0`$ or $`ϵ_{y0}=0`$ direction. Tracking studies of other cases will test the consistancy of the plateau model. The direction in the space of $`x_0,p_{x0},y_0,p_{y0}`$ is further defined by $`p_{x0}=0`$. The particle motion is 4 dimensional because of the presence of skew multipoles. $`j_{trns}`$ may then be considered to be a function of $`x_0`$. Fig. 3.1 shows $`j_{trns}`$ as a function of $`x_0`$ as found with tracking runs of 1e6 turns. The apparent stability limit using 1e6 turns and $`dx_0`$=.1mm is 22.0 mm. To make the plateaus more visible , one can reduce the search interval $`dx_0`$ to $`dx_0`$=.05 mm. These results are shown in Fig. 3.2. One sees that there are two plateaus in the region shown with $`j_{trns}`$ greater than or equal to 1e4, whose levels are located at $`j_{trns}`$= .9e6 , and 1e4 turns .The oscillations on the plateau near the stability boundary, appear to be smaller than those seen in the $`ϵ_{x0}=ϵ_{y0}`$ case.

Runs of about 2e7 turns were done to find the shape of the ’next plateau’ in the region where $`x_0`$ is less than or equal to 21.0 mm. Enough tracking runs were done to determine the beginning , and the level of the ’next plateau’. The results are shown in Fig. 3.3. Including these longer runs , one sees that the ’next plateau’ begins at x0=21.0mm and the level of the ’next plateau’ is about $`j_{trns}`$=8e6 turns. Here, the begining of the plateau is the edge at larger $`x_0`$ and the end is the edge at lower $`x_0`$. Th end of the ’next plateau’ is some what difficult to determine, as the adjacent plateau at lower $`x_0`$ has a high level of about 1e9 turns. The end of the ’ next plateau’ was taken to be at $`x_0`$=19.0 mm. The adjacent plateau at higher x0 is at the level of 1.5e6 turns and with the width of 1.8 mm. These two plateaus are seperated by .017 in $`1/Log(j_{trns})`$ which is smaller than the .054 found in the $`ϵ_{x0}`$=$`ϵ_{y0}`$ case. The results are also shown as $`1/Log(j_{trns})`$ versus $`x_0`$ plot in Fig. 3.4.

The data found for these two plateaus will be used to extrapolate and find the 1e9 plateau in this case. In the extrapolation, the results found in chapter 5 for the average plateau width and level seperation will be used. In this case, this leads to the assumptions that each plateau is 2.0mm wide and a plateau level seperation of .033 in $`1/Log(j_{trns})`$ will be used here.

The results are shown in Fig. 3.5, where two extrapolated plateaus are shown with the plateau levels 8.5e8, and 4.5e12 . The aperture for 1e9 turns was taken to be 17.1mm, the beginning of the plateau with a level of 4.5e12 turns. 17.1 mm is to be compared with the aperture of 21.9 mm found with runs of 1e6 turns. A loss of 4.8 mm or 22%.

## Chapter 4  Survival function along $`ϵ_{x0}=0`$ direction

Results will now be given for the $`x_0=0,p_{x0}=0`$ or $`ϵ_{x0}=0`$ direction.In order to be able to compare results with those of the above two cases, $`j_{trns}`$ will plotted against $`X_0=(\beta _{x0}/\beta _{y0})^{.5}y_0`$. Fig. 4.1 shows $`j_{trns}`$ as a function of $`X_0`$ as found with tracking runs of about 2e7 turns.The apparent stability limit using 1e6 turns and $`dx_0`$=.1mm is $`u_{sl}`$=24.9 mm. In Fig. 4.1 one can make out two plateaus with levels larger than 1e4 turns. The levels of these two plateaus are at 1.7e6 turns and 1.5e7 turns. The end of the 1.5e7 plateau was taken as $`X_0`$=21.4 mm. The results are also shown as a $`1/Log(j_{trns})`$ plot in Fig. 4.2.

The data found for these two plateaus will be used to extrapolate and find the 1e9 plateau in this case. In the extrapolation, the results found in chapter 5 for the average plateau width and level seperation will be used. In this case, this leads to the assumptions that each plateau is 2.0 mm wide and a plateau level seperation of .033 in $`1/Log(j_{trns})`$ will be used here.

The results are shown in Fig. 4.3, where two extrapolated plateaus are shown with the plateau levels 3.3e9, and 7.7e13 turns . Most of the particles on the 3.3e9 plateau will be assumed to survive 1e9 turns, and the aperture for 1e9 turns is then 21.4mm. 21.4mm is to be compared with the aperture of 24.9 mm found with runs of 1e6 turns. A loss of 3.5 mm or 17%.

## Chapter 5  Extrapolation parameters for the plateau model

The extrapolation depends on two parameters, the width of the plateaus and the seperation between consecutive plateau levels as measured in terms of $`1/Log(j_{trns})`$. The behaviour of these two parameters was studied by doing tracking runs for 3 different cases, corresponding to the three different directions in phase space, $`ϵ_{y0}`$=0, $`ϵ_{x0}`$=$`ϵ_{y0}`$ and $`ϵ_{x0}`$=0. The two parameters were measured for these 3 cases using runs of about 2e7 turns. The results are summarized in Table 5.1. Altogether, 7 plateaus were found and the two parameters for these 7 plateaus were measured. In Table 5.1 one sees that the plateau width, as measured as $`\mathrm{\Delta }X_0`$, $`X_0=(\beta _{x0}(ϵ_{x0}+ϵ_{y0}))^{.5}`$, is relatively constant with an average value of $`\mathrm{\Delta }X_0`$=2.00 mm. The plateau level seperation when measured in $`1/Log(j_{trns})`$ varies considerably with an average value of .033.

Based on the above results, it is proposed that in extrapolating the survival function, it is assumed that the extrapolated plateau widths are given by $`\mathrm{\Delta }X_0`$=2.00 mm, and the plateau level seperations in $`1/Log(j_{trns})`$ are .033.

## Chapter 6  Conclusions

Tracking studies lead to a model of the survival function, which pictures it as sequence of plateaus. Within the plateaus, the survival time in turns, $`j_{trns}`$, oscillates about some constant vaue of $`j_{trns}`$ which will be called the level of the plateau. Studying the survival function along different directions in phase space, using an older version of the RHIC lattice , one finds that the width of the plateaus, $`\mathrm{\Delta }X_0`$, $`X_0=(\beta _{x0}(ϵ_{x0}+ϵ_{y0}))^{.5}`$, remains roughly constant at about 2.00 mm. The seperation between the levels of adjacent plateaus has the same order of magnitude when measured in terms of the change in $`1/Log(j_{trns})`$, and has an average vaue of .033. Using, these results for the width of the plateaus and the seperation between plateau levels, one can extrapolate to estimate the location of the plateaus that correspond to longer survival times than can be found by tracking. For the case treated, it was found that a required survival time of 1e9 turns reduced the aperture by about 15% as compared to the aperture found by tracking using 1e6 turns.

The plateau model also leads to new criteria to be used in tracking studies to find the aperture for particles to survive a given number of turns. In the plateau model, one finds the first plateau whose level is higher than the given number of turns, in order to find the aperture for the given number of turns. This is to be compared with often used method, where one does a search starting at large amplitudes until one finds an amplitude that survives the given number of turns. In the latter method one cannot be sure that a finer search would not find unstable runs at smaller amplitudes or how frequently these unstable runs will occur. In the plateau model, while there may be unstable runs at smaller amplitudes, there is the assumption that they will not occur frequently.


text:
# 1 Introduction

## 1 Introduction

Non-thermal radiation is observed in several clusters of galaxies at frequencies varying from the radio \[see (Ensslin 1999) for a recent review\] to the UV to the soft and hard X-rays. These observations are hard to reconcile with our classical view that would interpret the radio radiation as synchrotron emission and the hard X-rays and UV radiation as ICS of the same electrons off the photons of the cosmic microwave background. This interpretation in general requires a large rate of injection of cosmic rays at recent times and intracluster magnetic fields smaller than those obtained by Faraday rotation measurements (Eilek 1999; Clarke et al. 1999) by a factor $`50100`$. These can be considered at least circumstantial evidence for some problems in the conventional interpretation, although additional observational tests are required in order to have more solid proofs.

The only events that seem to be able to provide the required energetics are mergers of clusters of galaxies. In this paper we study in detail all the channels that may contribute to the production of non-thermal radiation during a cluster merger, including the contribution of both primary and secondary electrons. Since in some of the clusters observed in radio and X-rays there does not seem to be evidence for an ongoing merger, we concentrate our attention on mergers occurred in the past. Although the calculations are carried out in the most general form, we apply them specifically to the case of the Coma cluster, for which a complete set of multiwavelength observations have been carried out. A more general discussion of our results and more technical details can be found in (Blasi 2000).

For the case of Coma, fitting the non-thermal multiwavelength observations has the following implications: 1) Magnetic fields $`50`$ times smaller than the ones measured by Eilek (1999) and Clarke et al. (1999) are necessary; 2) the gamma ray flux at $`100`$ MeV exceeds the EGRET upper limits; 3) the merger must have been just ended.

Since the main reason for these problems resides in the synchrotron plus ICS model, some attempts have been made to look for alternative interpretations of the hard X-ray excess, in particular invoking the bremsstrahlung emission from a non-thermal tail in the electron distribution, produced by stochastic acceleration (Blasi 2000a; Ensslin, Lieu and Biermann 1999; Dogiel 1999). The X-ray spectra can be fitted by this model without requiring small magnetic fields, although large injection rates of MHD waves are needed (compatible with the expectations in a merger of two clusters of galaxies).

The paper is structured as follows: in section 2 we describe mergers as particle accelerators; in section 3 we describe our calculations; in section 4 we discuss an alternative model of the non-thermal radiation produced during cluster mergers and its consequences. We conclude in section 5.

## 2 Cluster mergers as particle accelerators

During the merger between two clusters of galaxies, a typical gravitational energy of $`E_{merger}1.4\times 10^{64}`$ erg is available (we assumed here a typical total mass of a cluster to be $`5\times 10^{14}M_{}`$ and a distance of $`1.5`$ Mpc between the two clusters). During the approach, this energy is mainly converted into kinetic energy of the dark matter component, which is weakly interacting. However, at some point a strong shock is formed in the baryonic component and energy is transferred from dark matter to baryons (and electrons). This process is thought to produce the heating of the intracluster gas. A small fraction of this energy however can be converted into kinetic energy of particles out of equilibrium (non-thermal) through first order Fermi acceleration at the shock. Both electrons and protons are expected to be accelerated although the process should be more efficient for protons. A detailed description of this problem can be found in (Levinson 1994; McClements et al. 1997), but a brief description can be useful here: protons with low energy do interact with Alfvén waves, while electrons do not. Therefore protons can be more easily injected into the accelerator while electrons need to interact with some other kind of waves (for instance whistlers). Moreover, in order for a particle to be accelerated at the shock, its Larmor radius must be larger than the width of the shock, comparable with the Larmor radius of thermal protons. It is easy to see that this condition is realized only for electrons with energy larger than $`510`$ MeV (much larger than the typical cluster temperature). Electrons need therefore a preacceleration, otherwise the fraction of electrons injected in the shock region remains insignificantly small. We introduce here the coefficient $`\xi <1`$ as the ratio of the spectra of electrons and protons at injection at fixed energy. For the injection spectrum of electrons and protons we use here a power law in momentum with index $`\gamma =2.32`$, needed to fit the radio spectrum of Coma (this value is compatible with the compression ratios at the merger’s shocks observed in the simulations of Takizawa and Naito (2000)). These spectra are changed by propagation and losses effects, which we calculate by solving the full transport equations for primary electrons, protons and secondary electrons, generated by the decay of the charged pions due to $`pp`$ inelastic scattering in the intracluster medium \[see (Sarazin 1999) for a detailed description of the energy losses and radiative properties of primary electrons\].

The cosmic ray confinement of the proton component over cosmological scales (Berezinsky, Blasi and Ptuskin 1997;Volk, Aharonian and Breitschwerdt 1996; Colafrancesco and Blasi 1998) and the fast energy losses of electrons make the results very weakly dependent on the choice of the diffusion coefficient, with the exception of the maximum energies achievable in the acceleration. We use a diffusion coefficient derived by applying the quasi linear theory to a Kolmogorov spectrum of magnetic fluctuations (Blasi and Colafrancesco 1999) with the size of the largest eddy $`200`$ kpc (typical size between two galaxies). Higher maximum energies can be attained adopting a Bohm diffusion coefficient, but the fluxes far away from the tails do not change appreciably, as stressed above.

## 3 The calculations

Our calculations are fully time-dependent, therefore we can evaluate the non-thermal spectra at each time during or after the merger. Some general findings are the following: i) sufficiently high energy electrons, required to generate radio radiation and hard X-rays, exist in the cluster only if the merger has just ended or is currently ongoing. In particular the higher frequencies in the radio spectrum can be produced only if the merger was not over more than $`20`$ million years ago. ii) The secondary electron component is time independent: in fact the protons injected during the merger are confined in the intracluster medium and serve as a continuous source of new electrons, even if the proton injection ended long ago. iii) The gamma ray fluxes due to secondary electrons or due to pion decays are time independent. They cannot be washed out by the time evolution.

On this basis, we assume here that the merger just ended in Coma and calculate the fluxes of non-thermal radiation. The results for the radio emission of Coma are plotted in Fig. 1, where a magnetic field of $`0.15\mu G`$ was required. The thick curve that fits the observations (Feretti et al. 1995) refers to synchrotron emission of primary electrons. The cutoff is due to the maximum energy of the electrons and therefore might be in a different place for different choices of the diffusion coefficient. The real presence of a cutoff at $`>2`$ GHz was questioned anyway by Deiss et al. (1997). The two thin lines represent the synchrotron contribution of the secondary electrons for $`\xi =0.1`$ (lower curve) and $`\xi =0.01`$ (upper curve).

The results of our calculations for the hard X-rays and UV-soft X-ray radiation are plotted in Fig. 2, together with the expectation due to thermal bremsstrahlung of a gas of electrons at the temperature of Coma (thick solid curve). The data points are from BeppoSAX (Fusco-Femiano et al. 1999). The thick band is an estimate of the UV excess (Lieu et al. 1999). The thin solid line is the ICS contribution of primary electrons while the two dashed lines represent the ICS radiation of secondary electrons for $`\xi =0.1`$ (lower curve) and $`\xi =0.01`$ (upper curve).

With the parameters used to fit the radio and X-ray fluxes, we also calculate the gamma ray emission, due to several channels. The results are plotted in Fig. 3. The thick solid curve is the bremsstrahlung contribution of primary electrons, the two solid thin lines represent the fluxes of gamma rays from pion decay, the dashed lines are ICS fluxes from secondary electrons and the dashed dotted lines are the bremsstrahlung fluxes from secondary electrons. The fluxes produced by secondaries are always plotted for $`\xi =0.1`$ (lower curves) and $`\xi =0.01`$ (upper curves). The EGRET upper limit (Sreekumar et al. 1996 - arrow in the figure) is exceeded by a factor 3-4 . If a Bohm diffusion coefficient is used, then the maximum energy of primary electrons becomes large enough to generate an appreciable flux of gamma rays due to ICS. In this case the excess is a factor $`15`$.

Adopting magnetic fields even slightly larger than $`0.15\mu G`$ causes the radio radiation to exceed observations if the X-ray flux (e.g. the cosmic ray normalization) is fixed. Viceversa, if the magnetic field is increased and we fit the radio observations, the X-ray fluxes become too small.

## 4 Cluster mergers: alternative views

In the previous section we emphasized that the conventional interpretation of the non-thermal phenomena observed in clusters of galaxies leads to implications that seem to be in contrast with observations: a) small magnetic fields compared with Faraday Rotation measurements; b) large gamma ray fluxes.

There is an alternative way of explaining observations (Blasi 2000a; Dogiel 1999) without these problems: the idea consists of the following key points (Blasi, 2000a): 1) magnetic fields in the intracluster medium are at $`\mu G`$ level, as indicated by Faraday rotation measurements; 2) radio radiation is generated by synchrotron emission of relativistic electrons, with a very small energy content; 3) during the merger the electron thermal distribution is changed due to resonant interactions with perturbations in the magnetic field and acquire a non-Maxwellian tail; 4) hard X-rays are the result of bremsstrahlung of this modified electron distribution.

Detailed calculations of the development of this tail were carried out by Blasi (2000a): all processes responsible for the thermalization, resonant interaction with waves, and energy losses were included in the form of terms in a time-dependent non-linear Fokker-Planck (FP) equation, which was then solved numerically. This is equivalent to study the process of thermalization of a plasma in the presence of a perturbed magnetic field. We find that the electron distribution is not a Maxwell-Boltzmann distribution, which is a rather general result. There are two aspects to keep in mind: first, the thermalization time of electrons at low energies is extremely short, therefore the only region where possible distorsions from a thermal distribution can occur is on its tail; second, the FP equation is such that even if the energy transfer from the waves to electrons occurs only on the tail (lower energy electrons do not resonate), also the rest of the distribution is affected.

As a consequences of the points just stressed, the resonant interaction of electrons with waves in the intracluster medium mainly results into two effects: 1) the bulk of the electrons is heated up (the effective temperature increases); 2) a non-Maxwellian tail develops.

The results of Blasi (2000a) are plotted in Fig. 4a for the electron spectra and 4b for the X-ray spectra obtained as bremsstrahlung emission of the modified electron distribution. The data points are from BeppoSAX while the upper limits are from OSSE (Rephaeli et al. 1994). The different curves refer to different times.

This alternative model can be tested: since the electron spectrum is changed, the spectrum of the upscattered photons of the microwave background radiation (the so-called Sunyaev-Zeldovich effect) is also affected. Blasi, Olinto and Stebbins (2000) calculated in detail this modified SZ effect and found that: i) accurate measurements of the SZ spectrum should show distorsions distinguishable from the ordinary SZ effect; ii) these changes in the SZ spectra would affect the estimate of the Hubble constant and other cosmological parameters by $`20\%`$; iii) the masses estimated from X-rays should be a factor $`2`$ smaller than those derived from gravitational lensing, sensitive to the total energy (mass plus thermal energy plus non-thermal energy) in the cluster.

## 5 Conclusions

Radio radiation, hard and soft X-ray excesses and UV radiation have been observed from several clusters of galaxies, as a clear demonstration of the presence of cosmic rays and magnetic fields in the intracluster medium. The energetics involved, in terms of rate of injection of cosmic rays in clusters, seems to be compatible with the ones expected during mergers of clusters of galaxies. For this reason, in this paper we investigated the propagation and the radiative processes of both primary and secondary electrons and of protons, and we calculated the expected fluxes of radio, X and gamma radiation, applying our results to the specific case of Coma, for which a complete set of multiwavelength observations is available.

The general conclusion is that primary electrons can provide a good fit to the observations, provided the merger just ended (or is still ongoing). The rate of injection of primary electrons needed to fit the X-ray data is $`23\%`$ of the total energy available during the merger, and a magnetic field of $`0.15\mu G`$ is required to fit the radio data. The flux of UV radiation is automatically obtained in the right order of magnitude. All these fluxes are rapidly fading away with time, due to energy losses of high energy electrons. The contribution of the secondary electrons is slightly smaller but time independent, due to the confinement of cosmic ray protons on cosmological time scales.

The main consequence of the model is that the gamma ray fluxes above $`100`$ MeV exceed the EGRET upper limit by a factor between 3-4 to $`15`$ depending on the diffusion coefficient, which determines the maximum energy of the primary electrons. Besides this, the use of a magnetic field of $`0.15\mu G`$ seems to be inconsistent with the results of recent Faraday rotation measurements (Eilek 1999; Clarke et al. 1999).

There are two things to do in order to clarify the situation: 1) carry out a new measurement of the gamma ray fluxes in the GeV range, in order to check that EGRET did not just miss the detection; 2) carry out gamma ray measurements in the range $`>100`$ GeV, as proposed by Blasi (1999). These measurements can give precious information about the non-thermal content of clusters and they can be carried out even with current detectors like STACEE (Ong 1998).

A more theoretical approach consists in looking for alternative interpretations of the non-thermal radiation. We discussed this possibility in section 3, where we illustrated the calculations of Blasi (2000a), relative to the spectrum of a population of electrons that thermalize in the presence of a perturbed magnetic field. By solving the FP equation, Blasi (2000a) found that the electron spectrum presents an high energy tail, and that the bremsstrahlung emission from these electrons can explain the whole spectrum of X-rays, from the thermal region to the hard X-ray region. Magnetic fields compatible with the ones obtained by Faraday rotation measurements and low gamma ray fluxes can be accomodated naturally within this approach.

Moreover it is possible to prove (or disprove) this alternative model by looking at precision measurements of the Sunyaev-Zeldovich (SZ) effect: the electron distribution, modified by the resonant interactions with waves in the intracluster magnetic field produces a modified SZ effect that can be distinguished by the ordinary one, as discussed by Blasi, Olinto and Stebbins (2000).

Aknowledgments This work was supported by the DOE and the NASA grant NAG 5-7092 at Fermilab.

## 6 References

Berezinsky, V.S., Blasi, P., Ptuskin, V.S.: 1997 Astrophys. J. 487, p. 529.

Blasi, P.: 2000 submitted to Astropart. Phys.

Blasi, P.: 2000a Astrophys. J. Lett. 532, p. L9.

Blasi, P.: 1999 Astrophys. J. 525, p. 603.

Blasi, P., Colafrancesco, S.: 1999 Astropart. Phys. 12, p. 169.

Blasi, P., Olinto, A.V., Stebbins, A.: 2000 Astrophys. J. Lett. in press.

Clarke, T.E., Kronberg, P.P., Böringer, H.: 1999 in Ringberg Workshop on “Diffuse Thermal and Relativistic Plasma in Galaxy Clusters”, Eds: H. Böringer, L. Feretti, P. Schuecker, MPE Report No. 271.

Colafrancesco, S., Blasi, P.: 1998 Astropart. Phys. 9, p. 227.

Deiss, R.M., Reich, W., Lesch, H., Wielebinski, R.: 1997 A&A 321, p. 55.

Dogiel, V.A.: 1999 in Ringberg Workshop on “Diffuse Thermal and Relativistic Plasma in Galaxy Clusters”, Eds: H. Böringer, L. Feretti, P. Schuecker, MPE Report No. 271, 1999.

Eilek, J.A.: 1999 preprint astro-ph/9906485, to appear in ’Diffuse Thermal and Relativistic Plasma in Galaxy Clusters’, 1999, Ringberg Workshop, Germany, MPE-Report.

Ensslin, T.A.: 1999 Proceedings of The Universe at Low Radio Frequencies, ASP Conference Series, (preprint astro-ph/0001433).

Ensslin, T.A., Lieu, R., Biermann, P.: 1999 A&A 344, p. 409.

Feretti, L., Dallacasa, D., Giovannini, G., Taglianai, A.: 1995 A&A 302, p. 680.

Fusco-Femiano, R., Dal Fiume, D., Feretti, L., Giovannini, G., Grandi, P., Matt, G., Molendi, S., Santangelo, A.: 1999 Astrophys. J. Lett. 513, p. L21.

Levinson, A.: 1994 Astrophys. J. 426, p. 327.

Lieu, R., Ip, W.-H., Axford, W.I., Bonamente, M.: 1999 Astrophys. J. 510, p. 25.

McClements, K.G., Dendy, R.O., Bingham, R., Kirk, J.G., Drury, L. O’C.: 1997 MNRAS 291, p. 241.

Rephaeli, Y., Ulmer, M., Gruber, D.: 1994 Astrophys. J. 429, p. 554.

Sarazin, C.L.: 1999 Astrophys. J. 520, p. 529.

Ong, R.A.: 1998 Phys. Rep. 305, p. 93.

Sreekumar, P. et al.: 1996 Astrophys. J. /bf 464, p. 628.

Takizawa, M., Naito, T.: 2000 preprint astro-ph/0001046 (accepted for publication in Astrophys. J.).

Volk, H.J., Aharonian, F.A., Breitschwerdt, D.: 1996 Space Sci. Rev. 75, p 279.